检测OTU序列遗传发育信号的R实现

  群落格局及其背后的构建机制是群落生态学研究的重点内容，Stegen et al. 在2011年使用了beta-NTI这一指标来推断群落构建，此后得到了微生物生态学家的广泛使用。使用beta-NTI推断群落构建有一个前提条件，就是该基因必须在较小的遗传距离上具有遗传发育信号(Q1:为什么这是必须的呢？)，即在较小的遗传距离下，遗传发育距离越大，其生态位距离越大(Q2:为什么是较小的遗传距离下有遗传距离发育信号即可，不需要整个的遗传距离呢？)。本文将介绍三种方法检测近距离下的遗传发育信号，并用R语言实现。

如果没有环境因子数据可以使用方法一和方法三，如果有环境因子数据，三种方法都可使用。

1. 栖息地偏好

1.1 栖息地偏好介绍

  假如有两个OTU在各个样品(不要包含重复)中的分布如下：

OTU1	OTU2
3	4
2	5
6	3
6	2

  那么这两个OTU在样品空间的距离就是栖息地距离，这里我们可以用Bray-Curtis距离作为衡量指标，也可以使用其它距离测度（可参考vegdist()函数的说明文档）。

1.2 R语言计算OTU间栖息地距离并与遗传发育距离比较

#用函数打包便于代码复用
phylosignal <- function(otu_niche,otu_tree,nclass){
  #加载所需的R包
  library(picante)
  library(ape)
  #去除tree中不在OTU表的OTU分枝
  prune_tree<-prune.sample(otu_niche,otu_tree)
  #去除OTU表中不在prune_tree中的OTU
  tip<-prune_tree$tip.label
  coln<-colnames(otu_niche)
  m<-NULL
  for(i in 1:length(coln)){
    if(!coln[i]%in%tip){
      print(coln[i])
      m<-cbind(m,coln[i])
    }
  }
  m<-as.vector(m)
  otu_niche<-otu_niche[,!colnames(otu_niche)%in%m]
  #计算OTU间的遗传发育距离
  otu_phydist <- cophenetic(prune_tree)
  #otu_bray<-vegdist(log1p(otu_niche), method="bray")
  #计算OTU间的栖息地距离
  otu_bray<-vegdist(t(otu_niche), method="bray")
  otu_bray<-as.matrix(otu_bray)
  length(otu_bray) == length(otu_phydist)
  #mntd(otu_niche,cophenetic(prune_tree))
  #确保两个距离矩阵的OTU顺序相同
  otu_bray <- otu_bray[colnames(otu_phydist),colnames(otu_phydist)]
  #整理数据并绘图
  data1 <- data.frame(cor = c(unlist(otu_bray)),phydist=c(unlist(otu_phydist)))
  colnames(data1) <- c("sorted_otu_bray_3","sorted_otu_phydist_3")
  data1<-data1[order(data1$sorted_otu_phydist_3),]
  fac1<-cut(data1$sorted_otu_phydist_3,nclass)
  data2<-cbind.data.frame(fac1,data1)
  library(plyr)
  data3<-ddply(data2,.(fac1),colwise(median))
  plot(data3$sorted_otu_phydist_3,data3$sorted_otu_bray_3)
  physin <- data.frame(phydist=data3$sorted_otu_phydist_3,habpre=data3$sorted_otu_bray_3)
  physin2 <- list(physin,otu_bray,otu_phydist)
  names(physin2) <- c("physin","otu_bray","otu_phydist")
  return(physin2)
}

physin <- phylosignal(otu_niche,otu_tree,600)
#去除上三角和对角上的数值减少计算量，并消除对角值的影响
physin$otu_bray[upper.tri(physin$otu_bray, diag=TRUE)] = NA
physin$otu_phydist[upper.tri(physin$otu_phydist, diag=TRUE)] = NA
#该步骤耗时较长，建议输出变量保留结果。
man2 <- mantel.correlog(physin$otu_bray,physin$otu_phydist)
write.table(man2$mantel.res,"~/Desktop/phylogenetic signal.txt",sep = "\t",quote=FALSE,row.names=TRUE)

2. 重要生态位距离并与遗传发育距离比较

2.1 生态位介绍

如果某个细菌在pH=5的条件下长得最好，那么pH=5就是该细菌的最适生态位，其它的环境因子也类似。我们估算OTU生态位的步骤如下：

• 我们对OTU表进行基于距离的冗余分析；
• 向前筛选变量的方法确定影响显著的环境因子；
• 获得各OTUs在RDA三序图中显著影响的环境因子上的投影作为生态位指标(利用wascores()函数计算)；

2.2 R语言计算OTU间生态位距离

2.2.1 加载R包并读取数据

library(vegan)
nifhotu_norep <- read.delim("~/Desktop/nifhotu_no_reps.txt", row.names=1)
env_norep <- read.delim("~/Desktop/env_no_reps.txt", row.names=1)

2.2.2 进行db-RDA分析，并利用向前筛选变量的方法获得影响显著的环境因子

#distance calculation
bray_dist_norep <- vegdist(nifhotu_norep,method = "bray")
#rda analysis
otu_dbrda_norep <- dbrda(sqrt(bray_dist_norep)~.,data=env_norep,add = TRUE)
anova(otu_dbrda_norep,permutations = how(nperm = 999))
#forward selection
step_forward_norep <- ordistep(dbrda(sqrt(bray_dist_norep)~1,data=env_norep,add = TRUE),scope = formula(otu_dbrda_norep),direction="forward")

2.2.3 计算OTU间生态位距离

env_niche <- wascores(env_norep[,c(1:2,6)],nifhotu_norep)
pH.dist = as.matrix(dist(env_niche[,1]), labels=TRUE)
tn.dist = as.matrix(dist(env_niche[,2]), labels=TRUE)

2.2.3 计算OTU间遗传发育距离并与OTU间生态位距离比较

同1.2

3. 物种关联性

3.1 物种关联性简介与计算

  在2.1 生态位介绍部分介绍的生态位主要反映的是基于环境因子的生态位，但物种交互作用也是生态位的一部分，因此物种关联矩阵可以更全面地衡量微生物间的生态位差异。物种关联矩阵可以通过spearman相关，或者sparCC软件获得。

3.2 关联矩阵与遗传发育距离比较

同1.2部分

测试数据及R代码可参考资源:
面包多:检测OTU序列遗传发育信号的测试数据与R代码
CSDN:检测OTU序列遗传发育信号的测试数据与R代码
运行效果可见：Test phylogenetic signal
更多R语言分析微生物生态学的资源可参考如下链接：
https://mbd.pub/o/bread/mbd-YpmTlZpr