Microbial Metagenomic Trait Statistical Analysis Workflow

Script

require(vegan)
require(ape)
require(MASS)
require(ggplot2)

mg = eval(parse(text=traitsIn));
mg.phylo = eval(parse(text=phyloIn));
mg.phylo.hel = eval(parse(text=phyloHelIn));
mg.func.hel = eval(parse(text=funcHelIn));


#####################################################################################
#16S rRNA PHYLOGENETIC BETA-DIVERSITY
#####################################################################################
#NMDS
mg.phylo.metanmds<-metaMDS(mg.phylo.hel) #Stress 10.30507
#pdf(phylo_nmds)
#plot(gos.phylo.metanmds, type="t")
#dev.off()
mg.phylo.nmds<-data.frame(mg.phylo.metanmds$points)
mg.phylo.nmds<-cbind(mg.phylo.nmds, mg[,1], rowSums(mg.phylo))
colnames(mg.phylo.nmds)[3]<-"Habitat"
colnames(mg.phylo.nmds)[4]<-"Num16SrRNA"
#cols<-c("Coastal" = "cyan1","Open Ocean" = "blue","Fresh Water" = "yellow", "Estuary" = "chartreuse", "Mangrove"="darkorange1", "Hypersaline"="red", "Warm Seep"="darkblue", "Reef"="cornflowerblue", "Harbor"="darkgoldenrod4")
cat(capture.output(mg.phylo.nmds), file=phyloNmds, sep="\n")
cat(capture.output(mg.phylo.metanmds), file=phyloNmdsLog, sep="\n")
png(phylo_nmds)
print(qplot(MDS1, MDS2, data=mg.phylo.nmds, size=Num16SrRNA, colour=Habitat, xlab="Axis 1", ylab="Axis 2")
      +scale_colour_brewer(type="qual", palette = "Paired")+theme_bw()+ opts(legend.key = theme_blank()))
dev.off()

#####################################################################################
#SUBSYSTEMS FUNCTIONAL BETA-DIVERSITY
#####################################################################################
#NMDS
mg.func.metanmds<-metaMDS(mg.func.hel)  #Stress 11.98248
#pdf("gos.func.nmds.names.pdf")
#plot(gos.func.metanmds, type="t")
#dev.off()
mg.func.nmds<-data.frame(mg.func.metanmds$points)
mg.func.nmds<-cbind(mg.func.nmds, mg[,1], mg$TotalMb)
colnames(mg.func.nmds)[3]<-"Habitat"
colnames(mg.func.nmds)[4]<-"TotalMb"
cat(capture.output(mg.func.nmds), file=funcNmds, sep="\n")
cat(capture.output(mg.func.metanmds), file=funcNmdsLog, sep="\n")
png(func_nmds)
#cols <- c("Coastal" = "cyan1","Open Ocean" = "blue","Fresh Water" = "yellow", "Estuary" = "chartreuse", "Mangrove"="darkorange1", "Hypersaline"="red", "Warm Seep"="darkblue", "Reef"="cornflowerblue", "Harbor"="darkgoldenrod4")
print(qplot(MDS1, MDS2, data=mg.func.nmds, size=TotalMb, colour=Habitat, xlab="Axis 1", ylab="Axis 2")
      + scale_colour_brewer(type="qual", palette = "Paired")+theme_bw()+ opts(legend.key = theme_blank()))
dev.off()

R Server

Script

require(vegan)
require(ape)
require(MASS)
require(ggplot2)

##########
rarefaction.shannon<-function(x, random)

{	
#Arguments
community<-x
number_randomizations<-random

#Parameters and transformations
community_sorted<-apply(community,1,sort)
number_communities<-dim(community_sorted)[2]
a<-vector(length=number_randomizations)
number_sp<-vector(length=number_communities)


#Outputs
result<-matrix(NA,nrow=2,ncol=number_communities) 
rownames(result)=c("Mean_Shannon","Sd_Shannon")
colnames(result)=colnames(community_sorted)

#Loop for communities
for (i in 1:number_communities)
{
commX<-community_sorted[,i]
y<-commX[commX>0]
number_sp[i]<-length(y)
}
subsample_number_sp<-min(number_sp)

#Another loop for communities
for (i in 1:number_communities)
{
commX<-community_sorted[,i]
y<-commX[commX>0]
	#Loop for randomizations
	for (l in 1:number_randomizations)
	{
		#Calculation of Shannon
		w<-sample(y, subsample_number_sp, replace=FALSE)
		total<-sum(w)
    	x<-w/total
    	z<- -x * log(x)
    	H<-sum(z)
    	a[l]<-H
	}
#Final Calculations
result[1,i]=mean(a)
result[2,i]=sd(a)
}
#Final Output
return (result)
}

############

mg<-read.table(traits, sep="\t", header=T, row.names=1)
mg.dinuc<-read.table(dinuc, sep="\t", header=T, row.names=1)
mg.codon<-read.table(codon, sep="\t", header=T, row.names=1)
mg.aa<-read.table(aa, sep="\t", header=T, row.names=1)
mg.func<-read.table(func, sep="\t", header=T, row.names=1)
mg.phylo<-read.table(phylo, sep="\t", header=T, row.names=1)


mg.func.hel<-decostand(mg.func, method="hellinger")
mg.phylo.hel<-decostand(mg.phylo, method="hellinger")

#Calculate new traits from PCA of additional tables
#Dinucleotides
mg.dinuc.pca<-prcomp(mg.dinuc, scale=F)

#Codons
mg.codon.pca<-prcomp(mg.codon, scale=F)

#Aa
mg.aa.pca<-prcomp(mg.aa, scale=F)

#Subsystems
mg.func.pca<-prcomp(mg.func.hel, scale=F)

#Phylogenetic
mg.phylo.pca<-prcomp(mg.phylo.hel, scale=F)


cat("Dinucleotide PCA summary\n", capture.output(summary(mg.dinuc.pca)), "\nCodon PCA summary", capture.output(summary(mg.codon.pca)), "\nAminoacid PCA summary", capture.output(summary(mg.aa.pca)),"\nPhylogenetic PCA summary", capture.output(summary(mg.phylo.pca)), "\nSubsystem PCA summary", capture.output(summary(mg.func.pca)), file=pcaLog, sep="\n") #79.1% PCA1

#Calculate new traits using Shannon (H) diversity index
mg.func.H<-t(rarefaction.shannon(mg.func, random=1000))
mg.phylo.H<-t(rarefaction.shannon(mg.phylo, random=1000))

mg<-cbind(mg, mg.dinuc.pca$x[,1])
colnames(mg)[ncol(mg)]<-"Dinucleotides" 
mg<-cbind(mg, mg.codon.pca$x[,1])
colnames(mg)[ncol(mg)]<-"Codon"
mg<-cbind(mg, mg.aa.pca$x[,1])
colnames(mg)[ncol(mg)]<-"Aa"
mg<-cbind(mg, mg.func.pca$x[,1])
colnames(mg)[ncol(mg)]<-"Functional" 
mg<-cbind(mg, mg.func.H[,1])
colnames(mg)[ncol(mg)]<-"Functional.H"
mg<-cbind(mg, mg.phylo.pca$x[,1])
colnames(mg)[ncol(mg)]<-"Phylogeny"
mg<-cbind(mg, mg.phylo.H[,1])
colnames(mg)[ncol(mg)]<-"Phylogeny.H"

write.table(mg, file=summary_traits, quote=F, sep="\t")


traitsOut = deparse(mg);
phyloHelOut = deparse(mg.phylo.hel)
funcHelOut = deparse(mg.func.hel)
dinucDfOut = deparse(mg.dinuc)
codonDfOut = deparse(mg.codon)
aaDfOut = deparse(mg.aa)
dinucDfOut = deparse(mg.dinuc)
phyloOut = deparse(mg.phylo)

R Server

Script

require(vegan)
require(ape)
require(MASS)
require(ggplot2)

mg = eval(parse(text=traitsIn));
mg.dinuc = eval(parse(text=dinucDfIn));
mg.codon = eval(parse(text=codonDfIn));
mg.aa = eval(parse(text=aaDfIn));
mg.phylo.hel = eval(parse(text=phyloHelIn));
mg.func.hel = eval(parse(text=funcHelIn));

#Create distance matrices
mg.traits.dist<- vegdist(decostand(mg[,4:ncol(mg)], method="standardize"), method="euclidean")
mg.dinuc.dist<-vegdist(mg.dinuc, method="euclidean")
mg.codon.dist<-vegdist(mg.codon, method="euclidean")
mg.aa.dist<-vegdist(mg.aa, method="euclidean")
mg.func.dist<-vegdist(mg.func.hel, method="bray")
mg.phylo.dist<-vegdist(mg.phylo.hel, method="bray")

#####################################################################################
#PHYLOGENETIC AND FUNCTIONAL RESPONSES
#####################################################################################
#Relationship between functional traits and phylogenetic/functional matrix
#Create a distance matrix for each individual trait
for (i in 4:ncol(mg)) {
  x<-vegdist(mg[,i], method="euclidean")
  assign(paste(colnames(mg)[i],"distance",sep="_"), x)
}

#Mantel test against phylogenetic distance matrix (controlling for functional response)
gc.mantel<-mantel.partial(mg.phylo.dist, GC_distance, mg.func.dist, method="spearman")
vargc.mantel<-mantel.partial(mg.phylo.dist, varGC_distance, mg.func.dist, method="spearman")
rRNA.mantel<-mantel.partial(mg.phylo.dist, rRNA_distance, mg.func.dist, method="spearman")
genes.mantel<-mantel.partial(mg.phylo.dist, Genes_distance, mg.func.dist, method="spearman")
ab.mantel<-mantel.partial(mg.phylo.dist, ABratio_distance, mg.func.dist, method="spearman")
tf.mantel<-mantel.partial(mg.phylo.dist, TF_distance, mg.func.dist, method="spearman")
class.mantel<-mantel.partial(mg.phylo.dist, Classified_distance, mg.func.dist, method="spearman")
dinuc.mantel<-mantel.partial(mg.phylo.dist, mg.dinuc.dist, mg.func.dist, method="spearman")
codon.mantel<-mantel.partial(mg.phylo.dist, mg.codon.dist, mg.func.dist, method="spearman")
aa.mantel<-mantel.partial(mg.phylo.dist, mg.aa.dist, mg.func.dist, method="spearman")
funcH.mantel<-mantel.partial(mg.phylo.dist, Functional.H_distance, mg.func.dist, method="spearman")
phylo.mantel<-mantel.partial(mg.phylo.dist, mg.phylo.dist, mg.func.dist, method="spearman")
phyloH.mantel<-mantel.partial(mg.phylo.dist, Phylogeny.H_distance, mg.func.dist, method="spearman")
traits.mantel<-mantel.partial(mg.traits.dist, mg.phylo.dist, mg.func.dist, method="spearman")

phylo.mantel.df<-cbind('GC content', gc.mantel$statistic, gc.mantel$signif)
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Variance of GC content', vargc.mantel$statistic, vargc.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Dinucleotides', dinuc.mantel$statistic, dinuc.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Number of rRNA/genome', rRNA.mantel$statistic, rRNA.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Number of genes/genome', genes.mantel$statistic, genes.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Codons', codon.mantel$statistic, codon.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Amino acids', aa.mantel$statistic, aa.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Acidic to basic amino acids ratio', ab.mantel$statistic, ab.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('% of Transcriptional factors', tf.mantel$statistic, tf.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('% of classified reads', class.mantel$statistic, class.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Functional diversity', funcH.mantel$statistic, funcH.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Taxonomic content', phylo.mantel$statistic, phylo.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Taxonomic diversity', phyloH.mantel$statistic, phyloH.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('All community traits', traits.mantel$statistic, traits.mantel$signif))

phylo.mantel.df<-as.data.frame(phylo.mantel.df)
#row.names(phylo.mantel.df)<-phylo.mantel.df$V1
#phylo.mantel.df$V1<-NULL

colnames(phylo.mantel.df)[1]<-'traits'
colnames(phylo.mantel.df)[2]<-'statistic'
colnames(phylo.mantel.df)[3]<-'signif'

write.table(phylo.mantel.df, file=phyloMantel, quote=F, sep="\t", row.names=F)

#Mantel test against functional (Subsystems) distance matrix (controlling for phylogenetic response)
gc.mantel<-mantel.partial(mg.func.dist, GC_distance, mg.phylo.dist, method="spearman")
vargc.mantel<-mantel.partial(mg.func.dist, varGC_distance, mg.phylo.dist, method="spearman")
rRNA.mantel<-mantel.partial(mg.func.dist, rRNA_distance, mg.phylo.dist, method="spearman")
genes.mantel<-mantel.partial(mg.func.dist, Genes_distance, mg.phylo.dist, method="spearman")
ab.mantel<-mantel.partial(mg.func.dist, ABratio_distance, mg.phylo.dist, method="spearman")
tf.mantel<-mantel.partial(mg.func.dist, TF_distance, mg.phylo.dist, method="spearman")
class.mantel<-mantel.partial(mg.func.dist, Classified_distance, mg.phylo.dist, method="spearman")
dinuc.mantel<-mantel.partial(mg.func.dist, mg.dinuc.dist, mg.phylo.dist, method="spearman")
codon.mantel<-mantel.partial(mg.func.dist, mg.codon.dist, mg.phylo.dist, method="spearman")
aa.mantel<-mantel.partial(mg.func.dist, mg.aa.dist, mg.phylo.dist, method="spearman")
func.mantel<-mantel.partial(mg.func.dist, mg.func.dist, mg.phylo.dist, method="spearman")
funcH.mantel<-mantel.partial(mg.func.dist, Functional.H_distance, mg.phylo.dist, method="spearman")
phyloH.mantel<-mantel.partial(mg.func.dist, Phylogeny.H_distance, mg.phylo.dist, method="spearman")
traits.mantel<-mantel.partial(mg.func.dist, mg.func.dist, mg.phylo.dist, method="spearman")

phylo.mantel.df<-cbind('GC content', gc.mantel$statistic, gc.mantel$signif)
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Variance of GC content', vargc.mantel$statistic, vargc.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Dinucleotides', dinuc.mantel$statistic, dinuc.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Number of rRNA/genome', rRNA.mantel$statistic, rRNA.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Number of genes/genome', genes.mantel$statistic, genes.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Codons', codon.mantel$statistic, codon.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Amino acids', aa.mantel$statistic, aa.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Acidic to basic amino acids ratio', ab.mantel$statistic, ab.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('% of Transcriptional factors', tf.mantel$statistic, tf.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('% of classified reads', class.mantel$statistic, class.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Functional content', func.mantel$statistic, func.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Functional diversity', funcH.mantel$statistic, funcH.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('Taxonomic diversity', phyloH.mantel$statistic, phyloH.mantel$signif))
phylo.mantel.df<-rbind(phylo.mantel.df, cbind('All community traits', traits.mantel$statistic, traits.mantel$signif))

phylo.mantel.df<-as.data.frame(phylo.mantel.df)
#row.names(phylo.mantel.df)<-phylo.mantel.df$V1
#phylo.mantel.df$V1<-NULL

colnames(phylo.mantel.df)[1]<-'traits'
colnames(phylo.mantel.df)[2]<-'statistic'
colnames(phylo.mantel.df)[3]<-'signif'

write.table(phylo.mantel.df, file=funcMantel, quote=F, sep="\t", row.names=F)

R Server

Script

require(vegan)
require(ape)
require(MASS)
require(ggplot2)

mg = eval(parse(text=traitsIn));

#Principal Component Analysis (PCA)
mg.pca <-prcomp(mg[,4:ncol(mg)], scale=T)
cat("Traits PCA summary\n", capture.output(summary(mg.pca)),file=traitsPcaLog, sep="\n") #79.1% PCA1
mg.trait.pca<-data.frame(mg.pca$x[,1:2])
mg.trait.pca<-cbind(mg.trait.pca, mg$Description, mg$TotalMb)
colnames(mg.trait.pca)[3]<-"Habitat"

#cols<-c("Coastal" = "cyan1","Open Ocean" = "blue","Fresh Water" = "yellow", "Estuary" = "chartreuse", "Mangrove"="darkorange1", "Hypersaline"="red", "Warm Seep"="darkblue", "Reef"="cornflowerblue", "Harbor"="darkgoldenrod4")
png(traitsPcaPlot)
print(qplot(PC1, PC2, data=mg.trait.pca, size=I(5), colour=Habitat, xlab="PCA1", ylab="PCA2")
      + scale_colour_brewer(type="qual", palette = "Paired")+theme_bw()+ opts(legend.key = theme_blank()))
dev.off()


png(traitsBiplot)
biplot(mg.pca,cex=0.5, col=c("white", "red"))
dev.off()

cor.prob <- function(X, dfr = nrow(X) - 2) {
  R <- cor(X, method="spearman")
  above <- row(R) < col(R)
  r2 <- R[above]^2
  Fstat <- r2 * dfr / (1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr)
  R
}
cor.prob (mg[,4: ncol(mg)])
#Plot (all against all)
png(bivariatePlot, height=2048, width=2048)
print(plotmatrix(mg[,4: ncol(mg)]) + geom_smooth(method="lm"))
dev.off()

R Server

Script

aaOut = null;
codonOut = null;
dinucOut = null;
funcOut = null;
phyloOut = null;
traitsOut = null;

errorOut = null;

String[] iports = {aa, codon, dinuc, func, phylo, traits};
String[] iportsDef = {"Aminoacid", "Codon", "Dinucleotides", "Functional diversity", "Phylogenetic diversity", "Traits"};

// Test if stream is empty
for (var o = 0; o < iports.length; o++){

try{


if (iports[o].length() == 0)
	throw new Exception("ERROR:" + iportsDef[o] + " file is empty.") ;	

}catch(Exception err){
errorOut = err;
return;
}
}






//analysis needs more than 3 metagenomes and 1 column



	
for (var o = 0; o < iports.length; o++){
	
try{
var rows = iports[o].split("\n");
if (rows.length < 5)	
	throw new Exception("ERROR:" + iportsDef[o] + " file needs more than 3 metagenomes.") ;
for (var i = 0; i < rows.length; i++) {
	field = rows[i].split("\t");
	if (field.length <= 1)
		throw new Exception("ERROR:" + iportsDef[o] + " file needs more than 1 column.") ;	
} 
	
}
catch (Exception err){
errOut = err; 
return;	
}

}




if (errorOut == null){
errorOut = "Data validation passed.";
aaOut = aa;
codonOut = codon;
dinucOut = dinuc;
funcOut = func;
phyloOut = phylo;
traitsOut = traits;	
}