?
本文首先发布于“宏基因组”公众号原创。
作者：舟行天下
编辑：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 dataotu <- t(otu.tab)#data normolization (Legendre and Gallagher,2001)##by logenv.data.log <- log1p(env.data)####delete NAenv <- na.omit(env.data.log)###hellinger transformotu.hell <- decostand(otu, “hellinger”)#DCA analysis sel <- decorana(otu.hell)selotu.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)#删除掉共线性的环境因子，删掉最大的变量，直到所有的变量都小于10otu.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 againotu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env.data.log)#方差膨胀因子分析,目前所有变量都已经小于10vif.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 ?
?
51224936
