基于模型的约束排序，并探究OTUs对pH的响应特征——单峰or线性？_python

排序是微生物生态学中分析微生物格局常用方法，包括非约束排序和约束排序。非约束排序可以很直观地观察各样点或处理间微生物群落结构的变化情况，常用的有基于各种生态距离（beta多样性，最常用的Bray-Curtis距离、Jaccard距离等）的NMDS、PCoA和t-SNE降维，基于卡方距离的对应分析（CA，correspondence analysis ）和基于欧氏距离的PCA分析。约束排序可以探索环境因子对微生物群落结构的影响常用的有：RDA、CCA或基于距离的RDA（db-RDA）。以上这些方法我们称之为基于距离的排序，这些方法历史悠久，都诞生在上个世纪，在生态学中得到了长时间的检验，具体的内容可以参考Legendre的《R在数量生态学中的应用》一书。这里我们主要介绍一种新的排序方法——基于模型的排序。

1. 基于模型的排序与基于距离的排序相比有哪些优缺点？

两者的优缺点可以笼统地概括为下面的表格：

方法	优点	缺点
基于距离的排序	运行速度快、历史悠久	生态学意义提炼有限
基于模型的排序	可获得更多生态学意义解释	运行速度慢

下面具体介绍基于模型的方法的优点，主要表现在以下5个方面：

对mean-variance关系进行了考虑，从而区分trends in location (the natural variation)和changes in dispersion (the associated mean–variance relationships)；
检验模型与我们假设的关系，并比较不同的模型；
确定系数的不确定性；
用于预测，可用该模型生成新的数据；
考虑空间和时间造成的自相关性；

看到这五点，你可能看蒙了吧，确实不好理解，建议参考英文原文中的introduction部分：

gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r
Fast model-based ordination with copulas
Distance-based multivariate analyses confound
location and dispersion effects

2. 基于模型的方法目前有哪些实现方法？

这里推荐三个比较新的或者比较快的或者比较全面的方法，并在接下来的部分用gllvm的方法进行约束排序，看看其中的优势。

gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r
Fast model-based ordination with copulas
A new joint species distribution model for faster and more accurate inference of species associations from big community data
How to make more out of community data? A conceptual framework and its implementation as models and software

3. 如何利用基于模型的方法进行约束排序？

基于模型的约束排序的难点并不在于R代码的编写，而是如何将自己的假设和想要验证的问题用模型的形式表达出来，并进行解释。这里以OTUs对环境因子pH的响应作为例子，pH在很多生态系统中都被证明是影响群落结构的主要影响因子，OTUs对环境因子如pH的相应有两种模式——线性相应和单峰响应。线性响应是指OTU随研究范围内pH的增加单调地增加或减少，单峰响应是指OTU随研究范围内pH的增加呈现先增加后减少的趋势，那么如何区别这两者呢？基于模型的约束排序可以实现这一目的。下面是R代码：

nifhotu_norep <- read.delim("~/Desktop/htl_data_analysis/R code/Effect of biological replicates on RDA/nifhotu_no_reps.txt", row.names=1)
env_norep <- read.delim("~/Desktop/htl_data_analysis/R code/Effect of biological replicates on RDA/env_no_reps.txt", row.names=1)

library(gllvm)
Y <- nifhotu_norep
Y <- Y[,colSums(Y>0)/nrow(Y) > 0.8]
summary(rowSums(Y>0))
X <- env_norep 
X <- scale(X)
#设置了两个约束潜变量：num.lv.c = 2
CGLLVM <- gllvm(Y,X=X,family="negative.binomial", num.lv.c = 2)
summary(CGLLVM)
ordiplot(CGLLVM, biplot=FALSE)
#设置了两个约束潜变量：num.lv.c = 2
QCGLLVM <- gllvm(Y,X=X,family="negative.binomial", num.lv.c = 2,quadratic = TRUE)
summary(QCGLLVM)
ordiplot(QCGLLVM, biplot=FALSE)

此处我们可以看到考虑到二次效应的模型QGLLVM比仅考虑一次效应模型的AIC值更低，因此QGLLVM更优。然后我们再看LV1与LV2的标准偏差，我们可以看到LV1的变异较大为0.4088，LV2仅为0.0501。因此与LV1相比，我们可以忽略LV2。

且与其他变量相比pH对LV1的影响最大：

我们提取OTUs对潜变量LV1的二次响应系数与标准差，判断OTUs对潜变量LV1的响应特征——线性or单峰？

theta <- as.data.frame(QCGLLVM$params$theta)
sigma <- as.data.frame(QCGLLVM$sd$theta)
#theta$`CLV1^2` + 2*sigma$`CLV1^2`>0
#说明二次方的系数为0在两倍的标准差区间内，
#因此接受“二次方的系数为0”这一零假设。
#这样的OTUs为线性响应OTUs，这里称为lr_otu(linear responsing OTUs);
#否则为单峰响应OTUs，称为um_otu(unimodal OTUs)
lr_otu <- rownames(theta)[theta$`CLV1^2` + 2*sigma$`CLV1^2`>0]
um_otu <- rownames(theta)[theta$`CLV1^2` + 2*sigma$`CLV1^2`<0]

通过上面的分析我们发现pH是影响群落结构变异的主要环境因子（这一点并不新鲜），同时我们根据OTUs对pH响应特征进行了划分，得到了单峰响应OTUs，um_otu和线性响应OTUs，lr_otu。然后再根据这样的划分进行进一步的分析，例如每个样地不同响应特征OTUs的比例、其遗传发育特征等等。

难道就这？基于模型的方法还可以用于计算residue co-occurrence network，将在下一期的博客中介绍哈。

测试数据及R代码可参考如下资源：