的文章，大家對(duì)其中圖繪制過(guò)程比較感興趣。一上午收到了超30條留言，累計(jì)收到41個(gè)小伙伴的留言求精講。
我們將花時(shí)間把此文的原始代碼整理并精講，祝有需要的小伙伴能有收獲。
本系列按原文4幅組圖，共分為4節(jié)。本文是第4節(jié)，隨機(jī)森林回歸。
之前我們用了三篇文章，對(duì)隨機(jī)森林的應(yīng)用、分類、回歸進(jìn)行講解和實(shí)戰(zhàn)如下：
今天以圖3中的一個(gè)子圖，來(lái)實(shí)踐一下沖擊圖的繪制。
前情提要
先回顧一下圖4的內(nèi)容。
哪些菌可以作為生育時(shí)間的?
圖4. 水稻生育期相關(guān)的微生物標(biāo)記物()。
A. 采用隨機(jī)森林方法在兩地點(diǎn)的兩品種樣本中鑒定了23個(gè)綱與生育時(shí)間相關(guān)。其中按貢獻(xiàn)度由大到小排序。其中的子圖為交叉驗(yàn)證評(píng)估的結(jié)果。
B. 熱圖展示23個(gè)年齡相關(guān)的相對(duì)豐度。
方法簡(jiǎn)介：本圖A采用R語(yǔ)言的包進(jìn)行分析，結(jié)果采用的柱狀圖進(jìn)行可視化，按貢獻(xiàn)度由大到小排序時(shí)間序列分析實(shí)例，并進(jìn)行交叉驗(yàn)證模型的準(zhǔn)確度和數(shù)量的選擇依據(jù)。圖B采用展示每個(gè)時(shí)間點(diǎn)的相對(duì)豐度均值，其中按出現(xiàn)最高豐度的時(shí)間排序。
回歸分析
統(tǒng)計(jì)分析，主要基于兩個(gè)表：OTU表和實(shí)驗(yàn)設(shè)計(jì)表時(shí)間序列分析實(shí)例，對(duì)于想進(jìn)一步討論分類級(jí)，別需要OTUs的物種注釋文件。
這樣基于這3個(gè)文件，可以制作出千變?nèi)f化的統(tǒng)計(jì)分析圖片，來(lái)作為論據(jù)支持你的文章(Story)。
時(shí)間序列做回歸，主要是想建模來(lái)預(yù)測(cè)其它樣品的生育時(shí)間。主要分為兩部分，訓(xùn)練集建模，測(cè)試集驗(yàn)證。
我們主要有兩個(gè)品種，種植在兩個(gè)地點(diǎn)。這里先以A50建模，IR24驗(yàn)證的方案來(lái)演示。本實(shí)驗(yàn)較復(fù)雜 ，具體的方法會(huì)有多種組合。
讀取文件
# 讀取實(shí)驗(yàn)設(shè)計(jì)、和物種分類文件
tc_map =read.table("../data/design.txt",header = T, row.names = 1)
# 物種分類文件，由qiime summarize_taxa.py生成，詳見擴(kuò)增子分析流程系列
# 本研究以綱水平進(jìn)行訓(xùn)練，其實(shí)各層面都可以，具體那個(gè)層面最優(yōu)，需要逐個(gè)測(cè)試尋找。推薦綱、科，不建議用OTU，差異過(guò)大
otu_table =read.table("../data/otu_table_tax_L3.txt",header = T, row.names = 1)
# 篩選品種作為訓(xùn)練集
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)] ?隨機(jī)森林回歸
library(randomForest)
set.seed(315)
rf = randomForest(t(sub_otu), sub_map$day, importance=TRUE, proximity=TRUE, ntree = 1000)
print(rf)
# 模型結(jié)果如下
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交叉驗(yàn)證
## 交叉驗(yàn)證選擇Features
set.seed(315) # 隨機(jī)數(shù)據(jù)保證結(jié)果可重復(fù)，必須
# rfcv是隨機(jī)森林交叉驗(yàn)證函數(shù)：Random Forest Cross Validation
result = rfcv(t(sub_otu), sub_map$day, cv.fold=10)
result$error.cv查看錯(cuò)誤率表，23時(shí)錯(cuò)誤率最低，為最佳模型
?184 ? ? ? ?92 ? ? ? ?46 ? ? ? ?23 ? ? ? ?12 ? ? ? ? 6 ? ? ? ? 3 ? ? ? ? 1
116.56613 109.18251 100.78435 ?99.51937 110.35282 171.82919 286.35414 679.31080 ? ?
繪制驗(yàn)證結(jié)果
with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))重要性
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)
# 簡(jiǎn)單可視化
varImpPlot(rf, main = "Top 23 - Feature OTU importance",n.var = 25, bg = par("bg"),
? ? ? ? ? color = par("fg"), gcolor = par("fg"), lcolor = "gray" )想繪制圖中樣式的圖，可以使用imp的值，和名稱中對(duì)應(yīng)的物種注釋進(jìn)行數(shù)據(jù)篩選即可。這屬于美化工作，下面開始，個(gè)人風(fēng)格，僅供參考。
美化貢獻(xiàn)度柱狀圖
軟件內(nèi)部的可以快速可視化貢獻(xiàn)度，簡(jiǎn)單全面，但發(fā)表還是要美美噠，美是需要代碼的，就是花時(shí)間
基本思路同繪制Top 23 柱狀圖，按門著色，簡(jiǎn)化綱水平名字
# 讀取所有feature貢獻(xiàn)度
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保持隊(duì)形
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. 繪制時(shí)間序列熱圖
# 加載熱圖繪制包
library(pheatmap)
# 數(shù)據(jù)篩選23個(gè)feature展示
sub_abu = sub_otu[rownames(imp),]
# 直接自動(dòng)聚類出圖
pheatmap(sub_abu, scale = "row")這是樣品特征自動(dòng)聚類熱圖，已經(jīng)有時(shí)間序列的感覺(jué)，但樣本太多(n=249)，需要進(jìn)一步按組合并。
# 按時(shí)間為組合并均值
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)原文中又進(jìn)一步按各物種出現(xiàn)時(shí)間進(jìn)行排序，屬于個(gè)性化美化，這里不再贅述。
本分析的全部文件和代碼，會(huì)在 上持續(xù)更新，也可以后臺(tái)回復(fù)“時(shí)間序列”獲得百度云下載鏈接。
如果本文分享的技術(shù)幫助了你的科研，歡迎引用下文，支持國(guó)產(chǎn)SCI越來(lái)越好。
:
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 shift in rice
with time in the field and stage. Sci China Life Sci 61,