Step 7.4 Data visualization and analyses in R - shenjean/diversity GitHub Wiki

The commands below are run in the R environment. Export OTU and taxonomy tables from QIIME2 accordingly.

Import metadata, OTU and taxonomy tables

OTUtable=read.delim("biom/rplB.OTU.tsv",header=T,row.names=1,sep="\t")
TAXtable=read.delim("biom/rplB.TAX.tsv",header=T,row.names=1,sep="\t")

ampvis2

ampvis2 is a R package to visualize and analyze ordination and relative abundance data. For the metadata file, the sample IDs do not have to be converted to row names.

ampMETA=read.delim("singlem.metadata.tsv",header=T,sep="\t")

# Convert to ampvis2 object
d = amp_load(otutable=OTUtable,metadata=ampMETA,taxonomy=TAXtable)

Example commands for ordination analyses in ampvis2

PCA=amp_ordinate(d,type="PCA",transform="hellinger",species_plot=TRUE,
             sample_shape_by="depth_cat",sample_color_by="project_short",
             envfit_factor=c("combine"), envfit_show=FALSE,
             filter_species=0,detailed_output=TRUE)

CCA=amp_ordinate(d,type="CCA",transform="hellinger",species_plot=TRUE,
                 sample_color_by="depth_cat",sample_shape_by="project_short",
                 filter_species=0,detailed_output=TRUE,constrain=c("combine"))

phyloseq

phyloseq is a R package to visualize and analyze alpha diversity, ordination, and relative abundance data. Unlike ampvis2, the sample IDs should be converted to row names for the metadata file.

meta=read.delim("singlem.metadata.tsv",header=T,row.names=1,sep="\t")
# Convert to phyloseq object
sampledata=sample_data(meta)
rOTU = otu_table(OTUtable, taxa_are_rows = TRUE)
rTAX = tax_table(as.matrix(TAXtable))
physeq = phyloseq(rOTU, rTAX)
physeqr = merge_phyloseq(physeq, sampledata)

Check OTUs and samples with too many missing values. This is done using the WGCNA package

library(WGCNA)
gsg=goodSamplesGenes(t(otu_table(physeqr)))
gsg$allOK

Bray-Curtis ordination with cailliez/lingoes correction

rdist=ordinate(physeqr,distance="bray",method="PCoA",correction="lingoes")

plot_ordination(physeqr, rdist, "samples",color="project_short")+theme_bw()+
  geom_point(size=4)+ggtitle("PCoA of % recruitment of gammaproteobacterial MAGs based on Bray-Curtis distances")

# Check for negative eigenvalues 
hdist$values$Eigenvalues

pairwiseAdonis

pairwiseAdonis is a wrapper function for multilevel pairwise comparison using adonis2 (~Permanova) from package 'vegan'. The function returns adjusted p-values using p.adjust().

pairwise.adonis command does not accept interaction between factors neither strata

# Extract metadata from phyloseq data
metaR <- as(sample_data(physeqr), "data.frame")
# Run pairwise.adonis 
Radonis=pairwise.adonis(distance(physeqr,method="bray"),metaR$category,p.adjust.m="fdr")
write.table(Radonis,"pairwiseadonis.categories.tsv",sep="\t")

pairwise.adonis2 accepts interaction between factors and strata.

strata=pairwise.adonis2(distance(physeqr,method="bray")~depth_cat,data=metaR,strata='project_short',p.adjust.m="fdr")
strata