水稻微生物组时间序列分析4-随机森林回归

写在前面
哪些菌可以作为生育时间的biomarkers?
回归分析
- 读取文件
- 随机森林回归
- 交叉验证
- feature重要性
- 美化feature贡献度柱状图
- 图4.2. 绘制时间序列热图
猜你喜欢
写在后面

写在前面

之前分享了3月底发表的的
《水稻微生物组时间序列分析》的文章，大家对其中图绘制过程比较感兴趣。一上午收到了超30条留言，累计收到41个小伙伴的留言求精讲。

我们将花时间把此文的原始代码整理并精讲，祝有需要的小伙伴能有收获。

本系列按原文4幅组图，共分为4节。本文是第4节，随机森林回归。

之前我们用了三篇文章，对随机森林的应用、分类、回归进行讲解和实战如下：
- 一文读懂随机森林在微生态中的应用
- 随机森林randomForest 分类Classification
- 随机森林randomForest 回归Regression

今天以图3中的一个子图，来实践一下冲击图的绘制。

前情提要

水稻微生物组时间序列分析
1模式图与PCoA
2a-相关分析
2b-散点图拟合
3-冲击图

先回顾一下图4的内容。

哪些菌可以作为生育时间的biomarkers?

图4. 水稻生育期相关的微生物标记物(biomarkers)。

A. 采用随机森林方法在两地点的两品种样本中鉴定了23个纲与生育时间相关。其中按贡献度由大到小排序。其中的子图为交叉验证评估的结果。

B. 热图展示23个年龄相关的biomarkers相对丰度。

方法简介：本图A采用R语言的RandomForest包进行分析，结果采用ggplot2的柱状图进行可视化，biomarkers按贡献度由大到小排序，并进行交叉验证模型的准确度和biomarkers数量的选择依据。图B采用pheatmap展示每个时间点biomarkers的相对丰度均值，其中biomarkers按出现最高丰度的时间排序。

回归分析

统计分析，主要基于两个表：OTU表和实验设计表，对于想进一步讨论分类级，别需要OTUs的物种注释文件。

这样基于这3个文件，可以制作出千变万化的统计分析图片，来作为论据支持你的文章(Story)。

时间序列做回归，主要是想建模来预测其它样品的生育时间。主要分为两部分，训练集建模，测试集验证。

我们主要有两个品种，种植在两个地点。这里先以A50建模，IR24验证的方案来演示。本实验较复杂，具体的方法会有多种组合。

读取文件

# 读取实验设计、和物种分类文件
tc_map =read.table("../data/design.txt",header = T, row.names = 1)
# 物种分类文件，由qiime summarize_taxa.py生成，详见扩增子分析流程系列
# 本研究以纲水平进行训练，其实各层面都可以，具体那个层面最优，需要逐个测试寻找。推荐纲、科，不建议用OTU，差异过大
otu_table =read.table("../data/otu_table_tax_L3.txt",header = T, row.names = 1)
# 筛选品种作为训练集
sub_map = tc_map[tc_map$genotype %in% c("A50"),] # ,"IR24"
# 筛选OTU
idx = rownames(sub_map) %in% colnames(otu_table)
sub_map = sub_map[idx,]
sub_otu = otu_table[, rownames(sub_map)]

随机森林回归

library(randomForest)
set.seed(315)
rf = randomForest(t(sub_otu), sub_map$day, importance=TRUE, proximity=TRUE, ntree = 1000)
print(rf)# 模型结果如下
Call:randomForest(x = t(sub_otu), y = sub_map$day, ntree = 1000, importance = TRUE,      proximity = TRUE) Type of random forest: regressionNumber of trees: 1000
No. of variables tried at each split: 61Mean of squared residuals: 97.60524% Var explained: 92.97

交叉验证

## 交叉验证选择Features
set.seed(315) # 随机数据保证结果可重复，必须
# rfcv是随机森林交叉验证函数：Random Forest Cross Validation
result = rfcv(t(sub_otu), sub_map$day, cv.fold=10)
result$error.cv

查看错误率表，23时错误率最低，为最佳模型

  184        92        46        23        12         6         3         1
116.56613 109.18251 100.78435  99.51937 110.35282 171.82919 286.35414 679.31080

绘制验证结果

with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))

feature重要性

imp= as.data.frame(rf$importance)
imp = imp[order(imp[,1],decreasing = T),]
head(imp)
write.table(imp,file = "importance_class.txt",quote = F,sep = '\t', row.names = T, col.names = T)
# 简单可视化
varImpPlot(rf, main = "Top 23 - Feature OTU importance",n.var = 25, bg = par("bg"),color = par("fg"), gcolor = par("fg"), lcolor = "gray" )

想绘制图中样式的图，可以使用imp的值，和名称中对应的物种注释进行数据筛选即可。这属于美化工作，下面开始，个人风格，仅供参考。

美化feature贡献度柱状图

软件内部的varImpPlot可以快速可视化贡献度，简单全面，但发表还是要美美哒，美是需要代码的，就是花时间

基本思路同绘制Top 23 feature柱状图，按门着色，简化纲水平名字

# 读取所有feature贡献度
imp = read.table("importance_class.txt", header=T, row.names= 1, sep="\t")
# 分析选择top23分组效果最好
imp = head(imp, n=23)
# 反向排序X轴，让柱状图从上往下画
imp=imp[order(1:23,decreasing = T),]# imp物种名分解
# 去除公共部分
imp$temp = gsub("k__Bacteria;p__","",rownames(imp),perl=TRUE)
# 提取门名称
imp$phylum = gsub(";[\\w-\\[\\]_]+","",imp$temp,perl=TRUE) # rowname unallowed same name
imp$phylum = gsub("[\\[\\]]+","",imp$phylum,perl=TRUE)
# 提取纲名称
imp$class = gsub("[\\w\\[\\];_]+;c__","",imp$temp,perl=TRUE)
imp$class = gsub("[\\[\\]]+","",imp$class,perl=TRUE)
# 添加纲level保持队形
imp$class=factor(imp$class,levels = imp$class)# 图1. 绘制物种类型种重要性柱状图p=ggplot(data = imp, mapping = aes(x=class,y=X.IncMSE,fill=phylum)) + geom_bar(stat="identity")+coord_flip()+main_theme
p
ggsave(paste("rf_imp_feature",".pdf", sep=""), p, width = 4, height =4)

图4.2. 绘制时间序列热图

# 加载热图绘制包
library(pheatmap)# 数据筛选23个feature展示
sub_abu = sub_otu[rownames(imp),]# 直接自动聚类出图
pheatmap(sub_abu, scale = "row")

这是样品特征自动聚类热图，已经有时间序列的感觉，但样本太多(n=249)，需要进一步按组合并。

# 按时间为组合并均值
sampFile = as.data.frame(sub_map$day2,row.names = row.names(sub_map))
colnames(sampFile)[1] = "group"
mat_t = t(sub_abu)
mat_t2 = merge(sampFile, mat_t, by="row.names")
mat_t2 = mat_t2[,-1]
mat_mean = aggregate(mat_t2[,-1], by=mat_t2[1], FUN=mean) # mean
otu_norm_group = do.call(rbind, mat_mean)[-1,]
colnames(otu_norm_group) = mat_mean$group
pheatmap(otu_norm_group,scale="row",cluster_cols = F, cluster_rows = T)
pheatmap(otu_norm_group, scale = "row", filename = "heatmap_groups.pdf", width = 5, height = 5)

原文中又进一步按各物种出现时间进行排序，属于个性化美化，这里不再赘述。

本分析的全部文件和代码，会在 https://github.com/YongxinLiu/RiceTimeCourse 上持续更新，也可以后台回复“时间序列”获得百度云下载链接。

如果本文分享的技术帮助了你的科研，欢迎引用下文，支持国产SCI越来越好。

Citition:
Zhang, J., Zhang, N., Liu, Y.X., Zhang, X., Hu, B., Qin, Y., Xu, H., Wang, H., Guo, X., Qian, J., et al. (2018). Root microbiota shift in rice correlates
with resident time in the field and developmental stage. Sci China Life Sci 61, https://doi.org/10.1007/s11427-018-9284-4

点我下载PDF

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1700+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA