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

写在前面

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

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

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

之前我们用了三篇文章，对随机森林的应用、分类、回归进行讲解和实战如下：

• 随机森林randomForest 分类Classification

• 随机森林randomForest 回归Regression

今天以图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: regression                     Number of trees: 1000 No. of variables tried at each split: 61      Mean 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

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：生信宝典 学术图表 高分文章 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

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

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

点击阅读原文，跳转最新文章目录阅读