lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

R包randomForest的随机森林方法鉴定组间代表物种并进行样本分类

已有 7929 次阅读 2019-1-21 00:46 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _机器学习, _分类

RrandomForest的随机森林方法鉴定组间代表物种并进行样本分类

 

此处结合微生物群落研究中的16S扩增子分析数据,给大家分享怎样在R中使用随机森林分析数据。

R中,可用于进行随机森林分析的R包有很多可供选择,此处作为示例介绍其中的“randomForest”包,该R包的官方介绍详见该链接:

https://www.stat.berkeley.edu/~breiman/RandomForests/

 

首先介绍示例数据。我们此处共有13516S测序样本,这135个样本均来自土壤。其中,45个为对照组(c组),其余的95个为在土壤中添加了某种化学物质后搜集的。对于这95个添加了某种化学物质后的土壤样本,有45个样本中添加了高浓度的化学物质(h组);剩余的45个样本分为5组,每组9个,梯度添加了这种化学物质(m组,根据浓度梯度细分为m1m2m3m4m5)。

1.png

我们已经得知,对照组土壤的细菌群落与高浓度化学物质处理下的细菌群落具有显著差异。因此此处希望借助随机森林:

1)找出能够稳定区分对照组(c组)与高浓度处理组h组)OTU

2)构建模型,根据代表OTU的丰度,指示所添加化学物质的浓度。

 

作图示例文件、R脚本等,已上传至百度盘,无提取码(若失效请在下方留言

https://pan.baidu.com/s/1M91cNNN0h20pHbO1Fpne3A

同时,对于正文内容若有存在错误的地方,还请在下方留言告知,谢谢。



示例文件简要


示例文件“otu_table.txt”为OTU丰度表格,其内容展示如下。

每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度。

2.png

 

文件“group.txt”为样本分组信息,其内容展示如下。

第一列(names)为各样本名称,第二列(group)为各样本所属的分组信息,第三列为样本描述。对应上述样本分组文件,示例数据中共包含135例样本。其中,45个样本属于c分组,45个样本属于m分组(根据所添加化学物质的浓度梯度,具体分为m1m2m3m4m5),45个样本属于h分组。

3.png

 


数据预处理


首先导入数据。

#读取数据
otu <- read.table('otu_table.txt', sep = '\t', row.names = 1, header = TRUE, fill = TRUE)
group <- read.table('group.txt', sep = '\t', header = TRUE, fill = TRUE)[-3]
 
#作为演示,过滤低丰度(相对丰度总和 < 0.0005 丰度)群体,提高运算速率和准确度
#此处的 0.0005 丰度阈值为 979
otu <- otu[which(rowSums(otu) > 979), ]

 

然后我们将过滤后的OTU表以及样本分组文件进行合并,得到方便被randomForest识别计算的格式。

#合并
otu <- data.frame(t(otu))
otu$names <- rownames(otu)
otu_group <- merge(otu, group, by = 'names', all.x = TRUE)
rownames(otu_group) <- otu_group$names

此时我们查看合并后的数据框“otu_group”,在原OTU表的基础上进行了转置,这样每一列为一种OTU,每一行为一个样本,并将样本分组信息添加在了最后一列。

4.png

 

之后我们将所有数据中的m分组单独提取出来,得到一个新数据框“otu_train”;将剩余的95个样本(c组和h组)合并为另一新数据框“otu_test”。

#选择数据
otu_train <- subset(otu_group, group != 'm')[-1]
otu_test <- subset(otu_group, group == 'm')[-c(1, ncol(otu_group))]
otu_train$group <- factor(otu_train$group)

例如此处我们可选择查看其中一个数据框,不再多说了。

5.png

 

我们将一开始的数据分为了两部分,这样做的目的在于:

由于我们想要“寻找能够稳定区分对照组与高浓度处理组的OTU”,因此使用对照组(c组)数据集以及高浓度处理组(h组)数据集进行构建随机森林模型,这两个数据集在这里即为“训练集”,占总数据量的2/3;而低浓度梯度处理组(m组,占总数据量的1/3)可作为“测试集”来使用,以评估所构建的模型是否合理,同时实现另一目的“构建模型,根据代表OTU的丰度,指示所添加化学物质的浓度”。

因为m组中存在浓度梯度,所以理论上来讲,若我们所构建的随机森林模型合理,则通过模型对测试集进行分类时,可以将m组中的样本正确归类:m组中的较低浓度梯度处理(如浓度12)的样本更倾向被分类为“对照c组”,而较高浓度梯度处理(如浓度45)的样本更倾向被分类为“高浓度h组”。



构建随机森林模型,并初步确立代表OTU 


然后加载R包,构建随机森林模型。此处使用到randomForest包中的randomForest()函数,构建起来就一条命令(但是其本身的算法很复杂)……真实的应用中,需要仔细调整模型参数。

library(randomForest)
 
#随机森林计算(默认 500 树),randomForest()
#详细参数使用 ?randomForest 参看
otu_train.forest <- randomForest(group ~ ., data = otu_train, importance = TRUE)
 
#我们可简要查看结果,例如
otu_train.forest
#或
summary(otu_train.forest)
#或
varImpPlot(otu_train.forest, main = 'Top 30 - variable importance')

 

直接输入“otu_train.forest”后,打印在屏幕的简要输出结果中,我们看到了本次构建随机森林所使用的命令,构建的决策树的数量。

OOB estimate of error rate 表明了分类器的错误率,此处为为0%。其实在我们输入的OTU表中不难看出,在OTU(变量)种类和丰度组成上,c组和h组的差异非常明显,因此此处所构建的模型也比较精准。若是我们所期望分类的两组之间的变量差距不是很明显的话,越模糊错误率就会越高。

Confusion matrix表明了每个分类的详细的分类情况。例如,对于c 这个分组而言,基于随机森林算法的分类器,将45个样本分类到了c 这个分组,而对于被分类至c分组的这45个样本来讲,它们真实的分组恰好就是c分组,所以此处分类器的错误率为0h组同理。

此处由于我给的示例数据的组间差异实在太明显了,使得本次的随机森林模型准确度很高。其实在这里若结合有一定错误率的模型进行结果解读,可以更详细地帮助同学们对该模型结果的理解。想想还是懒得再换数据了,就继续往下来吧。

6.png

我们使用varImpPlot()命令简要展示了30个关键的OTU,将它们划分为“关键OTU”的依据为模型中的两个重要指标(两个指标下各自包含30OTU,默认由高往低排)。

此处“Mean Decrease Accuracy”和“Mean Decrease Gini”为随机森林模型中的两个重要指标。其中,“mean decrease accuracy”表示随机森林预测准确性的降低程度,该值越大表示该变量的重要性越大;“mean decrease gini”计算每个变量对分类树每个节点上观测值的异质性的影响,从而比较变量的重要性。该值越大表示该变量的重要性越大。

 

我们可以使用importance()命令,提取模型变量(每个 OTU)的重要性。

#提取模型变量(每个 OTU)的重要性,importance()
importance_otu <- data.frame(importance(otu_train.forest))
#可以考虑排个序,例如根据“Mean Decrease Accuracy”指标
importance_otu <- importance_otu[order(importance_otu$MeanDecreaseAccuracy, decreasing = T), ]

查看提取后的结果。根据排序后的结果,我们可以从中选择关键OTU

7.png

到了这一步,我们已经初步构建了一个准确的随机森林模型,并得到了一些关键OTU作为代表物种。

 


模型评估及测试


我们可以将构建好的模型进行自身验证,也可使用之前挑选的测试集进行测试,以更好地评估模型质量。

首先进行自身验证。

#自身测试
pred1 <- predict(otu_train.forest, data = otu_train)
freq1 <- table(pred1, otu_train$group)
sum(diag(freq1)/sum(freq1))
plot(margin(otu_train.forest, otu_test$group), main = '观测值被判断正确的概率图')

不难看出,对于绝大多数样本来讲,具有非常高的正确率。下图中,红点和蓝点分别代表c组和h组的样本。当然还是如先前所提及,我给的示例数据的组间差异实在太明显了,使得本次的随机森林模型准确度很高。若是样本间差异模糊,则在概率图中会观察到有相当一部分样本出现较高的分类错误。

8.png

 

然后使用测试集测试。

#测试数据测试
test_predict <- predict(otu_train.forest, otu_test)

如同我们所预期的那样,此处模型将m分组中的各样本进行分类。对于绝大多数样本,都能够被“正确分类”,即m组中的较低浓度梯度处理(如浓度12,对应样本名称m1.*m2.*)的样本更倾向被分类为“对照c组”,而较高浓度梯度处理(如浓度45,对应样本名称m4.*m5.*)的样本更倾向被分类为“高浓度h组”。

到了这一步,我们初步实现了构建模型的另一目的:“构建模型,根据代表OTU的丰度,指示所添加化学物质的浓度”。

9.png

 

我们可结合生态统计中的非约束排序分析,更直观地查看分类结果。

以下借助NMDS分析(非度量多维尺度分析)展示结果,关于NMDS分析此处不再进行更多的说明,不了解的同学们可以百度下。

##随便绘个排序图看下......
library(vegan)         #用于NMDS排序分析
library(ggplot2)      #借助ggplot2绘制排序图 
library(plyr)            #用于在绘图时绘制“分组边界”,得到geom_polygon()做图数据,使用其中的ddply()命令
 
#OTU 就用未过滤之前的
otu <- t(read.table('otu_table.txt', sep = '\t', row.names = 1, header = TRUE, fill = TRUE))
 
#替换一下 group 分组信息
group_train <- subset(group, group != 'm')
group_new <- data.frame(test_predict)
colnames(group_new) <- 'group'
group_new <- rbind(cbind(names = rownames(group_new), group_new, group2 = rep('test', ncol(group_new))), cbind(group_train, group2 = group_train$group))
 
#计算 Bray-Curtis 距离(beta 多样性分析常用距离),提取坐标,绘图展示等......
nmds <- metaMDS(otu, distance = 'bray') 
result <- nmds$points
result <- as.data.frame(cbind(result, rownames(result)))
colnames(result)[1:3] <- c('NMDS1', 'NMDS2', 'names')
result <- merge(result, group_new, by = 'names', all.x = TRUE)  
 
result$NMDS1 <- as.numeric(as.character(result$NMDS1)) 
result$NMDS2 <- as.numeric(as.character(result$NMDS2))
result$names <- as.character(result$names)
group2 <- ddply(result[-4], 'group2', function(df) df[chull(df[[2]], df[[3]]), ])
colnames(group2)[4] <- 'group'
 
ggplot(result[-5], aes(NMDS1, NMDS2, color = group)) +  
geom_polygon(data = group2, aes(fill = group), alpha = 0.2) +
geom_point()


NMDS排序图中,绿色的区块对应对照组(c组),蓝色的区块对应高浓度的化学物质处理后的分组(h组),红色的区块对应梯度处理后的分组(m组)。由于m组中的样本用作测试,因此它们被随机森林模型分类为“c”或“h”。

可以看到,m组中接近c组的点被划分为绿色,表示这些样本中的细菌群落组成与对照组的样本中的细菌群落组成更为接近;m组中接近h组的点被划分为蓝色,表示这些样本中的细菌群落组成与h组的样本中的细菌群落组成更为接近。

10.png

 


基于十折交叉验证,确立最终代表OTU


尽管我们在初始构建了一个较为准确的随机森林模型,不过我们仍然能够看出,还是有极少部分样本被“错误分类”。如“m1.2”,被错误划分至“h组”,即模型根据“m1.2”样本中的细菌OTU组成,将“m1.2”错误判定为“较高浓度梯度处理的样本”。因此,模型仍不完善,还需加以修正。

9.png

 

此处影响模型准确度的因素之一,就是某些OTU的存在带来了一定的错误率。由于这些OTU在不同分组(此处为c组和h组)中的丰度差异并不明显,即本身不具备组间代表性,导致在考虑了这些OTU去进行分类时,会带来一定的模糊识别。最好的方式是去除这些会带来误差的OTU

前面我们提到了一个重要的指标“Mean Decrease Accuracy”,可反映OTU本身的准确度。若该值过低,我们可以考虑将这些OTU去除,只保留重要的OTU(这些OTU可用作最终的组间代表OTU,对应我们在开头提到的第一个目的“找出能够稳定区分对照组与高浓度处理组的OTU”)。那么在这里,我们保留多少OTU合适呢?

首先我们根据计算得到的各OTUMean Decrease Accuracy”的值,将OTU由高往低排序,之后考虑执行重复5次的十折交叉验证(10-fold Cross Validation),根据交叉验证曲线对OTU进行取舍。交叉验证法的作用就是尝试利用不同的训练集/验证集划分来对模型做多组不同的训练/验证,来应对单独测试结果过于片面以及训练数据不足的问题。此处使用训练集本身进行交叉验证。

注:数据量越大,验证时运行速率越慢,需等待一会儿。此外,交叉验证每次所得结果可能不一致,因此可考虑根据多次结果取均值,或进行拟合。

library(doBy)  #方便基于分组取均值,使用其中的summaryBy()命令
library(reshape2)    #方便排列表格,使用其中的melt()命令
 
#5 次重复十折交叉验证
otu_train.5_10 <- replicate(5, rfcv(otu_train[-ncol(otu_train)], otu_train$group, cv.fold = 10,step = 1.5), simplify = FALSE)
 
#提取验证结果
otu_train.5_10 <- data.frame(sapply(otu_train.5_10, '[[', 'error.cv'))
otu_train.5_10$names <- rownames(otu_train.5_10)
otu_train.5_10 <- melt(otu_train.5_10, id = 'names')
otu_train.5_10$names<-as.numeric(as.character(otu_train.5_10$names))
otu_train.5_10 <- summaryBy(value~names, otu_train.5_10, FUN = mean)
 
#绘图展示交叉验证曲线
ggplot(otu_train.5_10, aes(names, value.mean)) +
geom_line(color = 'blue') +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +  
labs(title = '',x = 'Number of OTUs', y = 'Cross-validation error') +
geom_vline(xintercept = 20)

11.png

在这里,我们根据交叉验证曲线,可大致判定选取十几个OTU(根据“Mean Decrease Accuracy”值将OTU由高往低排)就具有足够的代表性了。

 

其实挺尴尬的,还是由于我给的示例数据的组间差异实在太明显了,使得交叉验证曲线很快下降至平缓,并未再出现小幅度上升。虽然在模型中,我们得到这样的结果其实是好事,说明模型可靠;但是考虑到在这里是作为分享教学,结果出现小波动才方便更好地解读……那这里就直接说吧,若是构建模型的训练集中,组间差异不是很明显的话,交叉验证曲线在下降至底部后,还会出现微小幅度的上升。这时,我们可认为导致曲线在后面出现微小幅度的上升的那些OTU会在一定程度上导致分类判别产生错误,无需添加至模型中。这时我们就可以根据曲线的特征,选取“谷值”对应的位置,将出现在此之前的OTU训练模型。

我的示例中,因为在一开始交叉验证曲线就直接降到底了,选几个OTU作为代表OTU好像也说不过去,那在这里就挑选20OTU吧。这些OTU即为最终得到的代表OTU

#之前我们已经根据“Mean Decrease Accuracy”指标对OTU排序过了
#从中选择 top20
importance_otu <- importance_otu[1:20, ]
otu_select <- rownames(importance_otu)

12.png

 


构建最终模型


到了上一步,“找出能够稳定区分对照组与高浓度处理组的OTU”这一目的最终完成了。

接下来,我们对原始OTU表进行筛选,只保留这20OTU,区域的全部过滤掉。之后仅基于这20OTU构建随机森林模型,并进行验证。

#分离 20 个重要 OTU,得到新的 OTU 丰度表,验证随机森林,方法同上
 
#选择数据
otu_train_20 <- otu_train[c(otu_select, 'group')]
otu_test_20 <- otu_test[otu_select]
 
#随机森林计算(默认 500 树) >>> randomForest()
otu_train.forest_20 <- randomForest(group ~ ., data = otu_train_20, importance = TRUE)
 
#测试数据测试
test_predict_20 <- predict(otu_train.forest_20, otu_test_20)
 
#比较两次(所有 OTU 数据集 & 20 OTU 数据集)测试结果
predict_two <- cbind(data.frame(test_predict), data.frame(test_predict_20))
for (i in 1:nrow(predict_two)) ifelse(predict_two[i,1] == predict_two[i,2], predict_two[i,3] <- 'True', predict_two[i,3] <- 'False')
#可见二者结果几乎一致,但第二次的结果质量明显提升
#然后懒得再画图展示了

 

以下是测试结果。其中,test_predict是使用所有OTU(过滤相对丰度总和 < 0.0005 的低丰度群体后)所得模型的测试结果;test_predict_20是使用20个挑选出来的代表OTU所得模型的测试结果;对于V3列,若test_predict=test_predict_20则为True,反之为False,方便快速比较两个结果的异同。可以看到,第二次的质量明显提升,例如,“m1.2”在第二次被正确归类。

13.png

 


参考文献


Robert I. Kabacoff. R语言实战(第二版)

Edwards J A , Christian M. Santos-Medellín, Liechty Z S , et al. Compositional shifts in root-associated bacterial and archaeal microbiota track the plant life cycle in field-grown rice. Plos Biology, 2018, 16(2):e2003862.

Root microbiota shift in rice correlates with resident time in the field and developmental stage. Science China(Life Sciences), 2018, v.61(06):3-11.




http://blog.sciencenet.cn/blog-3406804-1158196.html

上一篇:几种在NCBI中查询获取目的基因序列的方法
下一篇:R包vegan进行微生物群落主坐标分析(PCoA)及ggplot2作图示例

2 王敬敬 谢懿楠

该博文允许注册用户评论 请点击登录 评论 (3 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 18:07

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部