https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-standard-operating-procedure-v2

Note that this workflow is continually being updated. If you want to use the below commands be sure to keep track of them locally.

Last updated: 9 Oct 2018 (see revisions above for earlier minor versions and Metagenomics SOPs on the SOP for earlier major versions)

Important points to keep in mind:

这是那篇原始的文章：Arumugam, M., Raes, J., et al. (2011) Enterotypes of the human gut microbiome, Nature,doi://10.1038/nature09944 在谷歌上一搜，作者竟然做了个分析肠型的教程在这，学习一下：http://enterotyping.embl.de/enterotypes.html 这是2018年大佬们的共识文章：这是国人翻译的这篇文章，http://blog.sciencenet.cn/blog-3334560-1096828.html 当然，如果你只需要获得自己的结果或者自己课题的结果，不需要跑代码的，有最新的网页版分型，更好用，网址也放在这，同样也是上面翻译的那篇文章里提到的网址：http://enterotypes.org/ 只需要把菌属的含量比例文件上就能很快得到结果。

下面我就边学习边做来尝试着来个分析，并把代码放在这里备忘。其实作者已经整理好了代码，我学习一下，争取实现对手上的数据进行分析。

我表示对R语言还只是一知半解的状态，所以，先跑下，然后能用上自己的数据， 当个工具用就暂知足啦。我是黑苹果10.11的系统，运行这个软件提示少了Xquartz，于是装了个，windows和linux应该不需要。原代码中还提示『没有”s.class”这个函数』，百度了一下发现有个老兄的新浪博客说了是这个包，于是加了句library(ade4)就ok了。 Xquartz的下载地址Mac 10.6+：https://dl.bintray.com/xquartz/downloads/XQuartz-2.7.11.dmg

I get a list of the null count for each column:

What I’m trying to do is create a new dataframe which has the column header alongside the null count, e.g.

This is the complete workflow used to generate a random forest model using output data from an hmmsearch of your protein coding genes against eggNOG gamma proteobacterial protein HMMs. Before running this notebook, run the parse_hmmsearch.pl script to get a tab-delimited file containing bitscores for all isolates.

“`{r, read in data} # [...]]]> phenotype gproNOG.annotations bitscores

This is the complete workflow used to generate a random forest model using output data from an hmmsearch of your protein coding genes against eggNOG gamma proteobacterial protein HMMs. Before running this notebook, run the parse_hmmsearch.pl script to get a tab-delimited file containing bitscores for all isolates.

“`{r, read in data}
# library(gplots)
library(caret)
library(randomForest)
set.seed(1)
# set the directory you are working from
directory <- “/home/shenzy/UNISED/CII_ECOR_update_final_164genes”
# Reading in the eggNOG model scores, checking to make sure that top eggNOG model hit for each protein in the orthogroup is the same
traindata <- read.delim(paste(directory, “/bitscores.tsv”, sep=””))
traindata <- t(traindata)
phenotype <- read.delim(paste(directory, “/phenotype.tsv”, sep=””), header=F)
phenotype[,1] <- make.names(phenotype[,1])
traindata <- cbind.data.frame(traindata, phenotype=phenotype[match(row.names(traindata), phenotype[,1]),2])
traindata[is.na(traindata)] <- 0
# traindata <- na.roughfix(traindata)
traindata <- traindata[,-nearZeroVar(traindata)]
names(traindata) <- make.names(names(traindata))
“`

The following section is an optional step for picking the best values of mtry (number of gnees sampled per node) and ntree (number of trees in your random forest) for building your model. Instead of running all of the code at once, proceed through each step or model building and examine the figures produced. These will give you an indication of the point at which the performance of your model starts to level off.

In general, greater values of ntree and mtry will give you better stability in the top genes that are identified by the model. You can alternatively skip this step and proceed immediately to the next one, where values of 10,000 trees and p/10 genes per node (where p is total number of genes in the training data) have been chosen as a good starting point.

“`{r, train model}
# this section is for picking out the best parameters for building your model
set.seed(1)
# varying ntree
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=traindata, ntree=i)
error <- c(error, model$err.rate[length(model$err.rate)])
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(traindata)-1))
}
# varying mtry
error2 <- vector()
sparsity2 <- vector()
param <- ncol(traindata)-1
for(i in c(1, round(param/10), round(param/5), round(param/3), round(param/2), param)) {
model <- randomForest(phenotype ~ ., data=traindata, ntree=10000, mtry=i)
error2 <- c(error2, model$err.rate[length(model$err.rate)])
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(traindata)-1))
}
model <- randomForest(phenotype ~ ., data=traindata, ntree=10000, mtry=param/10, na.action=na.roughfix)
png(paste(directory, “/model_training/m1_error_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=error, xlab=”Number of trees”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m1_sparsity_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=sparsity, xlab=”Number of trees”, ylab=”% genes uninformative”, pch=16)
dev.off()
png(paste(directory, “/model_training/m1_error_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=error2, xlab=”Number of genes sampled per tree”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m1_sparsity_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=sparsity2, xlab=”Number of genes sampled per tree”, ylab=”% genes uninformative”, pch=16)
dev.off()
train2 <- traindata[,match(names(model$importance[model$importance[,1]>0,]), colnames(traindata))]
train2 <- cbind(train2, phenotype=traindata$phenotype)
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=train2, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(train2)-1))
}
error2 <- vector()
sparsity2 <- vector()
param <- ncol(train2)-1
for(i in c(1, round(param/10), round(param/5), round(param/3), round(param/2), param)) {
model <- randomForest(phenotype ~ ., data=train2, ntree=10000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(train2)-1))
}
model <- randomForest(phenotype ~ ., data=train2, ntree=10000, mtry=param/10, na.action=na.roughfix)
png(paste(directory, “/model_training/m2_error_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=error, xlab=”Number of trees”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m2_sparsity_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=sparsity, xlab=”Number of trees”, ylab=”% genes uninformative”, pch=16)
dev.off()
png(paste(directory, “/model_training/m2_error_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=error2, xlab=”Number of genes sampled per tree”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m2_sparsity_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=sparsity2, xlab=”Number of genes sampled per tree”, ylab=”% genes uninformative”, pch=16)
dev.off()
train3 <- train2[,match(names(model$importance[model$importance[,1]>quantile(model$importance[,1], 0.5),]), colnames(train2))]
train3 <- cbind(train3, phenotype=train2$phenotype)
error <- vector()
sparsity <- vector()
param <- ncol(train3)-1
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=train3, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(train3)-1))
}
error2 <- vector()
sparsity2 <- vector()
for(i in c(1, round(param/10), round(param/5), round(param/3), round(param/2), param)) {
model <- randomForest(phenotype ~ ., data=train3, ntree=10000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(train3)-1))
}
model <- randomForest(phenotype ~ ., data=train3, ntree=10000, mtry=param/10, na.action=na.roughfix)
png(paste(directory, “/model_training/m3_error_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=error, xlab=”Number of trees”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m3_sparsity_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=sparsity, xlab=”Number of trees”, ylab=”% genes uninformative”, pch=16)
dev.off()
png(paste(directory, “/model_training/m3_error_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=error2, xlab=”Number of genes sampled per tree”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m3_sparsity_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=sparsity2, xlab=”Number of genes sampled per tree”, ylab=”% genes uninformative”, pch=16)
dev.off()
train4 <- train3[,match(names(model$importance[model$importance[,1]>quantile(model$importance[,1], 0.5),]), colnames(train3))]
train4 <- cbind(train4, phenotype=train3$phenotype)
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=train4, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(train4)-1))
}
error2 <- vector()
sparsity2 <- vector()
param <- ncol(train4)-1
for(i in c(1, round(param/10), round(param/5), round(param/3), round(param/2), param)) {
model <- randomForest(phenotype ~ ., data=train4, ntree=10000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(train4)-1))
}
model <- randomForest(phenotype ~ ., data=train4, ntree=10000, mtry=param/10, na.action=na.roughfix)
png(paste(directory, “/model_training/m4_error_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=error, xlab=”Number of trees”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m4_sparsity_vs_ntree.png”, sep=””), width=350, height = 350)
plot(x=c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000), y=sparsity, xlab=”Number of trees”, ylab=”% genes uninformative”, pch=16)
dev.off()
png(paste(directory, “/model_training/m4_error_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=error2, xlab=”Number of genes sampled per tree”, ylab=”OOB error rate”, pch=16)
dev.off()
png(paste(directory, “/model_training/m4_sparsity_vs_mtry.png”, sep=””), width=350, height = 350)
plot(x=c(1, round(param/10), round(param/5), round(param/3), round(param/2), param), y=sparsity2, xlab=”Number of genes sampled per tree”, ylab=”% genes uninformative”, pch=16)
dev.off()
model$predicted
names(model$importance[order(model$importance, decreasing=T),][1:10])
save(model, train2, train3, train4, file=paste(directory, “/traindatanew.Rdata”, sep=””))
“`

This section allows you to build a model through iterative feature selection, using parameters that we feel are sensible. You can substitute in your own parameters chosen from the process above if you prefer.

“`{r, quicker model building}
set.seed(1)
param <- ncol(traindata)-1
model1 <- randomForest(phenotype ~ ., data=traindata, ntree=10000, mtry=param/10, na.action=na.roughfix)
model1
png(paste(directory, “/VI_fullnew.png”, sep=””), width=400, height=350)
plot(1:param, model1$importance[order(model1$importance, decreasing=T)], xlim=c(1,1000), ylab=”Variable importance”, xlab=”Top genes”)
dev.off()
pdf(paste(directory, “VI_fullnew.pdf”, sep=””), width=5, height=5)
plot(1:param, model1$importance[order(model1$importance, decreasing=T)], xlim=c(1,1000), ylab=”Variable importance”, xlab=”Top genes”)
dev.off()
train2 <- traindata[,match(names(model1$importance[model1$importance[,1]>0,]), colnames(traindata))]
train2 <- cbind(train2, phenotype=traindata$phenotype)
param <- ncol(train2)-1
model2 <- randomForest(phenotype ~ ., data=train2, ntree=10000, mtry=param/10, na.action=na.roughfix)
model2
train3 <- train2[,match(names(model2$importance[model2$importance[,1]>quantile(model2$importance[,1], 0.5),]), colnames(train2))]
train3 <- cbind(train3, phenotype=train2$phenotype)
param <- ncol(train3)-1
model3 <- randomForest(phenotype ~ ., data=train3, ntree=10000, mtry=param/10, na.action=na.roughfix)
model3
train4 <- train3[,match(names(model3$importance[model3$importance[,1]>quantile(model3$importance[,1], 0.5),]), colnames(train3))]
train4 <- cbind(train4, phenotype=train3$phenotype)
param <- ncol(train4)-1
model4 <- randomForest(phenotype ~ ., data=train4, ntree=10000, mtry=param/10, na.action=na.roughfix)
model4
train5 <- train4[,match(names(model4$importance[model4$importance[,1]>quantile(model4$importance[,1], 0.5),]), colnames(train4))]
train5 <- cbind(train5, phenotype=train4$phenotype)
param <- ncol(train5)-1
model5 <- randomForest(phenotype ~ ., data=train5, ntree=10000, mtry=param/10, na.action=na.roughfix, proximity=T)
model5
train6 <- train5[,match(names(model5$importance[model5$importance[,1]>quantile(model5$importance[,1], 0.5),]), colnames(train5))]
train6 <- cbind(train6, phenotype=train5$phenotype)
param <- ncol(train6)-1
model6 <- randomForest(phenotype ~ ., data=train6, ntree=10000, mtry=param/10, na.action=na.roughfix, proximity=T)
model6
train7 <- train6[,match(names(model6$importance[model6$importance[,1]>quantile(model6$importance[,1], 0.5),]), colnames(train6))]
train7 <- cbind(train7, phenotype=train6$phenotype)
param <- ncol(train7)-1
model7 <- randomForest(phenotype ~ ., data=train7, ntree=10000, mtry=param/10, na.action=na.roughfix, proximity=T)
model7
model7$predicted
names(model7$importance[order(model7$importance[,1], decreasing=T),])[1:10]
png(paste(directory, “final_model_VInew.png”, sep=””), width=400, height=350)

plot(1:param, model7$importance[order(model7$importance, decreasing=T),], xlab=””, ylab=”Variable importance”)
dev.off()
save(model1, model2, model3, model4, model5,model6, model7, traindata, train2, train3, train4, train5, train6, train7, file=paste(directory, “finalmodelnew.Rdata”, sep=””))
“`

The following section will show you how the performance of your model has improved as you iterated through cycles of feature selection. It will give you an idea of whether you have performed enough cycles, or whether you need to carry on.

You can look at both the out-of-bag votes, which are just votes cast by trees on data they were not trained on. This gives you a good idea of how the model would score strains that are distantly related to your training data. The second plot shows you votes cast by all of the trees on all of the serovars, and gives you a better idea of how your model would scores similar strains.

“`{r, looking at votes}
votedata <- rbind.data.frame(cbind.data.frame(model=rep(“1″, 9), Serovar=row.names(model1$votes), Invasive=model1$votes[,2]), cbind.data.frame(model=rep(“2″, 9), Serovar=row.names(model1$votes), Invasive=model2$votes[,2]), cbind.data.frame(model=rep(“3″, 9), Serovar=row.names(model1$votes), Invasive=model3$votes[,2]), cbind.data.frame(model=rep(“4″, 9), Serovar=row.names(model1$votes), Invasive=model4$votes[,2]), cbind.data.frame(model=rep(“5″, 9), Serovar=row.names(model1$votes), Invasive=model5$votes[,2]), cbind.data.frame(model=rep(“6″, 9), Serovar=row.names(model1$votes), Invasive=model6$votes[,2]), cbind.data.frame(model=rep(“7″, 9), Serovar=row.names(model1$votes), Invasive=model7$votes[,2]))
votedata$Phenotype <- phenotype[match(votedata$Serovar, phenotype[,1]),2]
ggplot(votedata, aes(x=model, y=Invasive, col=Phenotype)) + geom_jitter(width=0.1) + theme_classic(8) + ylab(“Proportion of votes for seletected strains phenotype”) + xlab(“Model iteration”) + geom_hline(yintercept=0.5, lty=2, col=”grey”) + theme(legend.key.size = unit(0.3, “cm”))
ggsave(“votesnew.png”, width=7, height=2.5)
ggsave(“votesnew.pdf”, width=7, height=2.5)
votedata$Phenotype <- phenotype[match(votedata$Serovar, phenotype[,1]),2]
ggplot(votedata, aes(x=model, y=Invasive, col=Serovar)) + geom_jitter(width=0.1) + theme_classic(8) + ylab(“Proportion of votes for seletected strains phenotype”) + xlab(“Model iteration”) + geom_hline(yintercept=0.5, lty=2, col=”grey”) + theme(legend.key.size = unit(0.3, “cm”))
ggsave(“votes.png”, width=7, height=2.5)
ggsave(“votes.pdf”, width=7, height=2.5)
fullvotes <- rbind.data.frame(cbind.data.frame(model=rep(“1″, 9), Serovar=row.names(model1$votes), Invasive=predict(model1, traindata, type=”vote”)[,2]), cbind.data.frame(model=rep(“2″, 9), Serovar=row.names(model1$votes), Invasive=predict(model2, train2, type=”vote”)[,2]), cbind.data.frame(model=rep(“3″, 9), Serovar=row.names(model1$votes), Invasive=predict(model3, train3, type=”vote”)[,2]), cbind.data.frame(model=rep(“4″, 9), Serovar=row.names(model1$votes), Invasive=predict(model4, train4, type=”vote”)[,2]), cbind.data.frame(model=rep(“5″, 9), Serovar=row.names(model1$votes), Invasive=predict(model5, train5, type=”vote”)[,2]), cbind.data.frame(model=rep(“6″, 9), Serovar=row.names(model1$votes), Invasive=predict(model6, train6, type=”vote”)[,2]), cbind.data.frame(model=rep(“7″, 9), Serovar=row.names(model1$votes), Invasive=predict(model7, train7, type=”vote”)[,2]))
fullvotes$Phenotype <- phenotype[match(fullvotes$Serovar, phenotype[,1]),2]
ggplot(fullvotes, aes(x=model, y=Invasive, col=Phenotype)) + geom_jitter(width=0.1) + theme_classic(8) + ylab(“Proportion of votes for invasive phenotype”) + xlab(“Model iteration”) + geom_hline(yintercept=0.5, lty=2, col=”grey”) + theme(legend.key.size = unit(0.3, “cm”))
ggsave(“full_votesnew.png”, width=7, height=2.5)
ggsave(“full_votesnew.pdf”, width=7, height=2.5)
“`

This section shows you how frequently particular genes are identified as top predictors across different iterations of model building, and how much importance they are assigned. Because model building involved a lot of stochasticity, you may find that your top predictors are quite subject to change across iterations.

“`{r, assessing stability of predictors}
set.seed(1)
usefulgenes <- data.frame()
topgenes <- data.frame()
for(i in 1:10) {
model <- randomForest(phenotype ~ ., data=traindata, ntree=10000, mtry=param/10, na.action=na.roughfix)
usefulgenes <- rbind(usefulgenes, cbind(model=i, gene=names(model$importance[model$importance>0,]), model$importance[model$importance>0]))
topgenes <- rbind(topgenes, cbind(model=i, gene=names(model$importance[order(model$importance, decreasing=T),][1:20]), model$importance[order(model$importance, decreasing=T),][1:20]))
}
png(paste(directory, “/gene_usefulnessnew.png”, sep=””))
hist(table(usefulgenes$gene), col=”grey”, main=””, xlab=”Number of times each gene was useful (VI >) in a model”)
dev.off()
sum(table(usefulgenes$gene)==10)
sum(table(usefulgenes$gene)<10)
topgenes$V3 <- as.numeric(as.character(topgenes$V3))
ggplot(topgenes, aes(x=model, y=gene, fill=as.numeric(V3))) + geom_tile() + ggtitle(“Imporance values for top genes across model iterations”) + scale_fill_continuous(“Importance”)
table(topgenes$gene)
png(paste(directory, “/topgenesnew.png”, sep=””))
hist(table(topgenes$gene), col=”grey”, main=””, xlab=”Number of times each gene appeared in the\ntop 20 predictors”)
dev.off()
save(usefulgenes, topgenes, file=paste(directory, “/allgenemodels.Rdata”, sep=””))
“`

This section allows you to look at the COG categories of your top predictor genes, compared to the COG categories in your original training set, to see if any are over-represented.

“`{r, COG analysis}
annotations <- read.delim(paste(directory, “/gproNOG.annotations.tsv”, sep=””), header=F) # matcing annotation file for your chosen model set can be downloaded from http://eggnogdb.embl.de/#/app/downloads
nogs <- read.delim(paste(directory, “/models_used_20ref.tsv”, sep=””), header=F)
nogs[,2] <- sub(“.*NOG\\.”, “”, nogs[,2])
nogs[,2] <- sub(“.meta_raw”, “”, nogs[,2])
info <- annotations[match(nogs[,2], annotations[,2]),]
info$V5 <- as.character(info$V5)
# proportion of each COG catergory that comes up in the top indicators
background <- nogs[match(colnames(traindata), make.names(nogs[,1])),2]
background_COGs <- annotations[match(background, annotations[,2]),5]
library(plyr)
bg_COGs2 <- aaply(as.character(background_COGs), 1, function(i){
if(!is.na(i)) {
if(nchar(i)>1) {
char <- substr(i,1,1)
for(n in 2:nchar(i)) {
char <- paste(char,substr(i,n,n),sep=”.”)
}
return(char)
} else {
return(i)
}
} else{
return(NA)
}
})

background_COGs2 <- unlist(strsplit(bg_COGs2, “[.]”))
predictors <- nogs[match(colnames(train6), make.names(nogs[,1])),2]
predictor_COGs <- annotations[match(predictors, annotations[,2]),5]
p_COGs2 <- aaply(as.character(predictor_COGs), 1, function(i){
if(!is.na(i)) {
if(nchar(i)>1) {
char <- substr(i,1,1)
for(n in 2:nchar(i)) {
char <- paste(char,substr(i,n,n),sep=”.”)
}
return(char)
} else {
return(i)
}
} else{
return(NA)
}
})
predictor_COGs2 <- unlist(strsplit(p_COGs2, “[.]”))
barplot(rbind(table(background_COGs2), table(predictor_COGs2)[match(names(table(background_COGs2)), names(table(predictor_COGs2)))]))
“`

This section allows you to compare the performance of your model predicting true phenotype and a random phenotype.

“`{r, control}
set.seed(2)
control <- traindata
control$phenotype <- sample(traindata$phenotype)
param <- ncol(control)-1
cmodel1 <- randomForest(phenotype ~ ., data=control, ntree=10000, mtry=param/10, na.action=na.roughfix)
cmodel1
png(paste(directory, “/VI_full_control.png”, sep=””), width=400, height=350)
plot(1:param, cmodel1$importance[order(cmodel1$importance, decreasing=T)], xlim=c(1,1000), ylab=”Variable importance”, xlab=”Top genes”)
dev.off()
pdf(paste(directory, “VI_full_control.pdf”, sep=””), width=5, height=5)
plot(1:param, cmodel1$importance[order(cmodel1$importance, decreasing=T)], xlim=c(1,1000), ylab=”Variable importance”, xlab=”Top genes”)
dev.off()
ctrain2 <- control[,match(names(cmodel1$importance[cmodel1$importance[,1]>0,]), colnames(control))]
ctrain2 <- cbind(ctrain2, phenotype=control$phenotype)
param <- ncol(ctrain2)-1
cmodel2 <- randomForest(phenotype ~ ., data=ctrain2, ntree=10000, mtry=param/10, na.action=na.roughfix)
cmodel2
ctrain3 <- ctrain2[,match(names(cmodel2$importance[cmodel2$importance[,1]>quantile(cmodel2$importance[,1], 0.5),]), colnames(ctrain2))]
ctrain3 <- cbind(ctrain3, phenotype=ctrain2$phenotype)
param <- ncol(ctrain3)-1
cmodel3 <- randomForest(phenotype ~ ., data=ctrain3, ntree=10000, mtry=param/10, na.action=na.roughfix)
cmodel3
ctrain4 <- ctrain3[,match(names(cmodel3$importance[cmodel3$importance[,1]>quantile(cmodel3$importance[,1], 0.5),]), colnames(ctrain3))]
ctrain4 <- cbind(ctrain4, phenotype=ctrain3$phenotype)
param <- ncol(ctrain4)-1
cmodel4 <- randomForest(phenotype ~ ., data=ctrain4, ntree=10000, mtry=param/10, na.action=na.roughfix)
cmodel4
ctrain5 <- ctrain4[,match(names(cmodel4$importance[cmodel4$importance[,1]>quantile(cmodel4$importance[,1], 0.5),]), colnames(ctrain4))]
ctrain5 <- cbind(ctrain5, phenotype=ctrain4$phenotype)
param <- ncol(ctrain5)-1
cmodel5 <- randomForest(phenotype ~ ., data=ctrain5, ntree=10000, mtry=param/10, na.action=na.roughfix, proximity=T)
cmodel5
ctrain6 <- ctrain5[,match(names(cmodel5$importance[cmodel5$importance[,1]>quantile(cmodel5$importance[,1], 0.5),]), colnames(ctrain5))]
ctrain6 <- cbind(ctrain6, phenotype=ctrain5$phenotype)
param <- ncol(ctrain6)-1
cmodel6 <- randomForest(phenotype ~ ., data=ctrain6, ntree=10000, mtry=param/10, na.action=na.roughfix, proximity=T)
cmodel6

cmodel6$predicted
names(cmodel6$importance[order(cmodel6$importance[,1], decreasing=T),])[1:10]
# compare votes and phenotypes for the control and real datasets
cbind(cmodel6$votes, control$phenotype, model6$votes, train6$phenotype)

“`

This approach repeats the full model-building process 5 times to see how similar the top predictor genes are.

“`{r, testing the robustness of the top result}
topgenes <- vector()
# picking out the top predictors from the model
for(i in 1:5) {
set.seed(i)
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=traindata, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(traindata)-1))
}
error2 <- vector()
sparsity2 <- vector()
param <- ncol(control)-1
for(i in c(1, param/10, param/5, param/3, param/2, param)) {
model <- randomForest(phenotype ~ ., data=traindata, ntree=1000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(traindata)-1))
}
train2 <- traindata[,match(names(model$importance[model$importance[,1]>0,]), colnames(traindata))]
train2 <- cbind(train2, phenotype=traindata$phenotype)
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=train2, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(train2)-1))
}
error2 <- vector()
sparsity2 <- vector()
param <- ncol(control)-1
for(i in c(1, param/10, param/5, param/3, param/2, param)) {
model <- randomForest(phenotype ~ ., data=train2, ntree=20000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(train2)-1))
}
train3 <- train2[,match(names(model$importance[model$importance[,1]>quantile(model$importance[,1], 0.5),]), colnames(train2))]
train3 <- cbind(train3, phenotype=train2$phenotype)
error <- vector()
sparsity <- vector()
for(i in c(1, 10, 50, 250, 500, 1000, 1500, 2000, 5000, 10000)) {
model <- randomForest(phenotype ~ ., data=train3, ntree=i, na.action=na.roughfix)
error <- c(error, median(model$err.rate))
sparsity <- c(sparsity, (sum(model$importance[,1]<=0))/(ncol(train3)-1))
}
error2 <- vector()
sparsity2 <- vector()
param <- ncol(control)-1
for(i in c(1, param/10, param/5, param/3, param/2, param)) {
model <- randomForest(phenotype ~ ., data=train3, ntree=2000, mtry=i, na.action=na.roughfix)
error2 <- c(error2, median(model$err.rate))
sparsity2 <- c(sparsity2, (sum(model$importance[,1]<=0))/(ncol(train3)-1))
}
train4 <- train3[,match(names(model$importance[model$importance[,1]>quantile(model$importance[,1], 0.5),]), colnames(train3))]
train4 <- cbind(train4, phenotype=train3$phenotype)
topgenes <- c(topgenes, colnames(train4))
}
table(topgenes)

随机森林是一种基于决策树（Decisiontree）的高效的机器学习算法，可以用于对样本进行分类（Classification），也可以用于回归分析（Regression）。

它属于非线性分类器，因此可以挖掘变量之间复杂的非线性的相互依赖关系。通过随机森林分析，可以找出能够区分两组样本间差异关键OTU。

Feature Importance Scores表格-来源于随机森林结果

记录了各OTU对组间差异的贡献值大小。

注：一般地，选取Mean_decrease_in_accuracy值大于0.05的OTU，作进一步分析；对于组间差异较小的样本，该值可能会降至0.03。

2. 交叉验证分析

交叉验证（Crossvalidation)，是一种统计学上将数据样本切割成较小子集的实用方法。先在一个子集上做分析，而其它子集则用来做后续对此分析的确认及验证。一开始的子集被称为训练集。而其它的子集则被称为验证集或测试集。

其中最常见的为k-foldercross-validation，它指的是将所有数据分成k个子集，每个子集均做一次测试集，其余的作为训练集。交叉验证重复k次，每次选择一个子集作为测试集，并将k次的平均交叉验证识别正确率作为结果。

所有的样本都被作为了训练集和测试集，每个样本都被验证一次。

对随机森林方法筛选出的关键OTU的组合进行遍历，以期用最少的OTU数目组合构建一个错误率最低高效分类器。

一般地，对随机森林分析筛选出的关键OTU，按照不同组合进行10倍交叉验证分析，找出能够最准确区分组间差异的最少的OTU组合，再做进一步的分析，如ROC分析等。

注：图中横坐标表示不同数量的OTU组合，纵坐标表示该数量OTU组合下分类的错误率。OTU组合数越少，且错误率越低，则该OTU组合被认为是能够区分组间差异的最少的OTU组合。

3. ROC曲线

接收者操作特征曲线（Receiveroperating characteristic curve，ROC 曲线）也是一种有效的有监督学习方法。ROC分析属于二元分类算法，用来处理只有两种分类的问题，可以用于选择最佳的判别模型,选择最佳的诊断界限值。

可依据专业知识，对疾病组和参照组测定结果进行分析，确定测定值的上下限、组距以及截断点(cut-offpoint)，按选择的组距间隔列出累积频数分布表，分别计算出所有截断点的敏感性(Sensetivity)、特异性和假阳性率(1-特异性:Specificity)。以敏感性为纵坐标代表真阳性率，(1-特异性)为横坐标代表假阳性率，作图绘成ROC曲线。ROC曲线越靠近左上角，诊断的准确性就越高。亦可通过分别计算各个试验的ROC曲线下的面积(AUC)进行比较，哪一种试验的AUC最大，则哪一种试验的诊断价值最佳。

注：图中横坐标为假阳性率false positive rate（FPR）：Specificity，纵坐标为真阳性率true positive rate（TPR）：Sensetivity。最靠近左上角的ROC曲线的点是错误最少的最好阈值，其假阳性和假阴性的总数最少。ROC曲线下的面积值在1.0和0.5之间。在AUC>0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在 0.5~0.7时有较低准确性，AUC在0.7~0.9时有一定准确性，AUC在0.9以上时有较高准确性。AUC=0.5时，说明诊断方法完全不起作用，无诊断价值。AUC<0.5不符合真实情况，在实际中极少出现。

4. Wilcoxon秩和检验分析

Wilcoxonrank-sum test，也叫曼-惠特尼U检验（Mann–WhitneyU test），是两组独立样本非参数检验的一种方法。其原假设为两组独立样本来自的两总体分布无显著差异，通过对两组样本平均秩的研究来实现判断两总体的分布是否存在差异，该分析可以对两组样品的物种进行显著性差异分析，并对p值计算假发现率（FDR）q值。

注：mean分别为两组样品物种的平均相对丰度，sd分别是两组样本物种相对丰度的标准差。P值为对两组检验原假设为真的概率值，p<0.05表示存在差异，p<0.01表示差异显著，q值为假发现率。

5. 差异菌群Heatmap分析

以10倍交叉验证（10-foldcross-validation）估计泛化误差（Generalizationerror）的大小，其余参数使用默认设置。建模结果同时包含“基线”误差（Baselineerror）的期望值，即数据集中属于最优势分类的样本全部被错误分类的概率。每个OTU根据其被移除后模型预报错误率增加的大小确定其重要度数值，重要度越高，该OTU对模型预报准确率的贡献越大。

根据挑选出来的差异OTU，根据其在每个样品中的丰度信息，对物种进行聚类，绘制成热图，便于观察哪些物种在哪些样品中聚集较多或含量较低。

注：图中越接近蓝色表示物种丰度越低，越接近橙红色表示丰度越高。左边的聚类树是根据各物种间的spearman相关性距离进行聚类；上边的聚类树是采用样本间距离算法中最常用的Bray-Curtis算法进行聚类。

6. 两组样本Welch’s t-test分析

两组不同方差的样本可使用Welch’st-test进行差异比较分析，通过此分析可获得在两组中有显著性差异的物种[或差异基因丰度—适用于元（宏）基因组]。

注：上图所示为不同基因丰度（或物种）在两组样品中的丰度比例，中间所示为95%置信度区间内，物种丰度的差异比例，最右边的值为p值，p值＜0.05，表示差异显著。

7. Shannon多样性指数比较盒状图

将不同分类或环境的多组样本的Shannon多样性指数进行四分位计算，比较不同样本组的组间Shannon指数差异。同时进行非参数Mann-Whitney判断样本组间的显著性差异

默认创建mkcp + 随机一种伪装头配置文件：

自定义v2ray配置文件:

查看v2ray配置:

warning: 如果用centos，需要先关闭防火墙

Alpha多样性（Alpha Diversity）是对某个样品中物种多样性的分析，包含样品中的物种类别的多样性——丰富度（Richness）和物种组成多少的整体分布——均匀度（Evenness）两个因素，通常用Richness,Chao1，Shannon，Simpson，Dominance和Equitability等指数来评估样本的物种多样性。

丰富度指数

Richness, Chao1，Shannon三个指数是常用的评估丰富度的指标，数值越高表明样品包含的物种丰富度就越高。

均匀度指数

Simpson，Dominance和Equitability三个指数是常用的评估均匀度的指标。

汇总表：

Next-generation sequencing (NGS) has become routine in oncology and opens up new avenues of therapies, particularly in personalized oncology setting. An increasing number of cases also implies a need for a more robust, automated, and reproducible processing of long lists of [...]]]> https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btz924/5674039

ClinVAP: A reporting strategy from variants to therapeutic options

Abstract