用RDA進行微生物環境因子分析
本文首先發佈於“巨集基因組”公眾號原創。
作者：舟行天下
編輯：metagenome
前言
在進行微生物多樣性分析時，大家一定會做α，Β多樣性分析。α多樣品通俗來講就是樣本內的物種多樣性。Β多樣性是指在地區尺度上，物種組成沿著某個梯度方向從一個群落到另一個群落的變化率。即沿著某一環境梯度，物種替代的速率、物種週轉率等。
排序的過程是將樣品或微生物物種排列在一定的空間， 使得排序軸能夠反映一定的生態梯度 這些排序方法又可以分成間接梯度排序(indirect gradient analysis)和直接梯度排序(direct gradient analysis)。間接梯度排序又叫非約束性排序；尋求潛在的或在間接的環境梯度來解釋物種資料的變化包括PCA，PCoA，NDMS，直接排序又叫約束性排序；它是指在特定的梯度上（環境軸） 上探討物種的變化情況；方法包括 RDA, CCA, db-RDA。排序分析（Ordination analysis）。排序(ordination)的過程就是在一個視覺化的低維空間或平面重新排列這些樣本,使得樣本之間的距離最大程度地反映出平面散點圖內樣本之間的關係資訊。
RDA 介紹
distance-based redundancy analysis (db-RDA) 是目前在微生物領域應用的最為廣泛的環境因子分析，該分析方法內建在R中的vegan包中。相信大家一定都知道vegan包，該R包是進行生態學（包括微生物多樣性分析）研究的必備神器！vegan包中提供了所以基本排序分析的方法，可以說是一包在手搞定所有！ 關於vegan包的詳細介紹，請大家檢視vegan包的官方文件
微生物環境因子分析
要進行微生物環境因子分析，我們需要兩個檔案，一個是微生物多樣性的OTU 表格，另一個就是你所有樣品的環境因子資料。 比如，你進行土壤微生物研究，這時候你就需要知道你所測土壤的C,N,P,K等化學元素含量以及不同樣地的氣候資訊等等，總之，在分析之前可以多準備些環境因子資料，後期我們還可以對這些環境因子進行共線性，以及環境因子與資料擬合優良性判斷。
資料均一化
首先看看我們準備的OTU表格以及環境因子資料結構（圖1，圖2），讀取完資料之後，我們要把OTU的橫軸和縱軸調換位置，然後把OTU表格也要進行hellinger轉化，使資料均一性更好。並把環境因子進行log轉化，以減少同一種環境因子之間本身數值大小造成的影響。
RDA和CCA模型篩選
資料都進行均一化之後，我們要進行RDA和CCA的模型篩選。先用species-sample資料做DCA分析看分析結果中Lengths of gradient的第一軸的大小，如果大於4.0，就應該選CCA，如果在3.0-4.0之間，選RDA和CCA均可，如果小於3.0，RDA的結果要好於CCA。（圖3）
方差膨脹因子分析
在篩選完RDA和CCA分析後，我們需要利用方差膨脹因子分析，對所有環境因子進行共線性分析。我們要依次刪掉最大的變數，也就是刪除掉共線性的環境因子，直到所有的變數都小於10。
檢測最低AIC值
最後我們要用step模型檢測最低AIC值，在這一步中該模型會自動篩選出最優的環境因子。當“none”位於最頂端時意味著改模型篩選結束，位於none值上方的環境因子即為與OTU擬合最好的環境因子。
ANOVA 顯著性分析並出圖
在進行完以上的資料篩選之後，我們可以用篩選的結果重新進行一次環境因子與OTU的線性迴歸分析，這樣我們就拿到了最終的計算結果，並且用ANOVA進行顯著性檢驗，並且通過該分析我們還可以看到所篩選的環境因子的整體貢獻率，以及每個環境因子的單獨貢獻率。
本例中我們使用了內建ggplot2的vegan—ggvegan進行的分析。ggvegan的出圖結果可以用內建的ggplot2進行優化，使你的圖更為美觀，其具體用法與ggplot2的圖層疊加方式類似。詳情大家可以參考ggvegan的官網
# 首先要安裝devtools包，僅需安裝一次
install.packages("devtools")
# 載入devtools包
library(devtools)
# 下載ggvegan包
devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
otu.tab <- read.csv("otutab.txt", row.names = 1, header=T, sep="\t")
env.data <- read.csv("new_meta.txt", row.names = 1, fill = T, header=T, sep="\t")
#transform data
otu <- t(otu.tab)
#data normolization (Legendre and Gallagher,2001)
##by log
env.data.log <- log1p(env.data)##
##delete NA
env <- na.omit(env.data.log)
###hellinger transform
otu.hell <- decostand(otu, "hellinger")
#DCA analysis
sel <- decorana(otu.hell)
sel
otu.tab.0 <- rda(otu.hell ~ 1, env) #no variables
#Axis 第一項大於四應該用CCA分析
otu.tab.1<- rda(otu.hell ~ ., env)
#我們在篩選完RDA和CCA分析後，我們需要對所有環境因子進行共線性分析，利用方差膨脹因子分析
vif.cca(otu.tab.1)
#刪除掉共線性的環境因子，刪掉最大的變數，直到所有的變數都小於10
otu.tab.1 <- rda(otu.hell ~ N+P+K+Ca+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)
vif.cca(otu.tab.1)
#進一步篩選
otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)
vif.cca(otu.tab.1)
#test again
otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env.data.log)
#方差膨脹因子分析,目前所有變數都已經小於10
vif.cca(otu.tab.1)
##用step模型檢測最低AIC值
mod.u <- step(otu.tab.0, scope = formula(otu.tab.1), test = "perm")# "perm"增加P值等引數
mod.d <- step(otu.tab.0, scope = (list(lower = formula(otu.tab.0), upper = formula(otu.tab.1))))
mod.d
##本處篩選的結果，找到一個Mg環境因子適合模型構建，為了下一步畫圖，我們
#保留所有非共線性的環境因子
#choose variables for best model and rda analysis again#
(otu.rda.f <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env))
anova(otu.rda.f)
anova(otu.rda.f, by = "term")
anova(otu.rda.f, by = "axis")
#計算db-rda
## 用ggvegan繪圖
p<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = c("point", "text"), layers = c( "species","sites", "biplot", "centroids"), legend.position = "right", title = "db-RDA")
## 新增圖層
p + theme_bw()+theme(panel.grid=element_blank())
Reference
- Rognes, T., Flouri, T., Nichols, B., Quince, C., & Mahé, F. (2016). VSEARCH: a versatile open source tool for metagenomics. PeerJ, 4, e2584.
- Edgar, R.C. (2013) UPARSE: Highly accurate OTU sequences from microbial amplicon reads, Nature Methods [Pubmed:23955772, dx.doi.org/10.1038/nmeth.2604].
- UNOISE2: Improved error-correction for Illumina 16S and ITS amplicon read. bioRxiv, 2016