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

1.加载R包和测试数据

library(vegan)
library(ape)
library(picante)
library(vegan)
nifhotu_norep <- read.delim("./data/nifhotu_no_reps.txt", row.names=1)
env_norep <- read.delim("./data/env_no_reps.txt", row.names=1)
tree <- read.tree("./data/tree")

2.让树文件与OTU表相匹配

tree$tip.label <- stringr::str_replace_all(tree$tip.label,"'","")
prune_tree<-prune.sample(nifhotu_norep,tree)
tip<-prune_tree$tip.label
coln<-colnames(nifhotu_norep)
m<-NULL
for(i in 1:length(coln)){
  if(!coln[i]%in%tip){
    print(paste("delete this OTU from tree:",coln[i]))
    m<-cbind(m,coln[i])
  }
}

[1] "delete this OTU from tree: OTU_10179"
[1] "delete this OTU from tree: OTU_10329"
[1] "delete this OTU from tree: OTU_113977"
[1] "delete this OTU from tree: OTU_122"
[1] "delete this OTU from tree: OTU_133570"
[1] "delete this OTU from tree: OTU_159"
[1] "delete this OTU from tree: OTU_22258"
[1] "delete this OTU from tree: OTU_24594"
[1] "delete this OTU from tree: OTU_26066"
[1] "delete this OTU from tree: OTU_267595"
[1] "delete this OTU from tree: OTU_35216"
[1] "delete this OTU from tree: OTU_38790"
[1] "delete this OTU from tree: OTU_43486"
[1] "delete this OTU from tree: OTU_45378"
[1] "delete this OTU from tree: OTU_53268"
[1] "delete this OTU from tree: OTU_57855"
[1] "delete this OTU from tree: OTU_87080"

m<-as.vector(m)
nifhotu_norep<-nifhotu_norep[,!colnames(nifhotu_norep)%in%m]

3.获得OTU间遗传发育距离

otu_phydist <- cophenetic(prune_tree)

4.基于生态位距离的遗传发育信号检验

4.1 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))

Permutation test for dbrda under reduced model
Permutation: free
Number of permutations: 999

Model: dbrda(formula = sqrt(bray_dist_norep) ~ pH + tn + tp + tk + ep + ek + emo + cec + den + nl + fl + sl + wc, data = env_norep, add = TRUE)
         Df SumOfSqs     F Pr(>F)    
Model    12   5.5318 1.184  0.001 ***
Residual 15   5.8400                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#forward selection
step_forward_norep <- ordistep(dbrda(sqrt(bray_dist_norep)~1,data=env_norep,add = TRUE),scope = formula(otu_dbrda_norep),direction="forward")

Start: sqrt(bray_dist_norep) ~ 1 

      Df    AIC      F Pr(>F)   
+ pH   1 68.271 2.7168  0.005 **
+ tp   1 69.375 1.6063  0.005 **
+ nl   1 69.639 1.3473  0.015 * 
+ tn   1 69.727 1.2613  0.040 * 
+ fl   1 69.815 1.1756  0.060 . 
+ emo  1 69.897 1.0967  0.200   
+ den  1 69.931 1.0633  0.245   
+ wc   1 69.958 1.0380  0.285   
+ tk   1 69.959 1.0363  0.335   
+ cec  1 69.971 1.0246  0.370   
+ sl   1 70.008 0.9889  0.465   
+ ep   1 70.144 0.8588  0.955   
+ ek   1 70.194 0.8106  0.990   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step: sqrt(bray_dist_norep) ~ pH 

      Df    AIC      F Pr(>F)   
+ tn   1 68.827 1.3225  0.005 **
+ tp   1 68.994 1.1668  0.060 . 
+ den  1 69.017 1.1445  0.075 . 
+ fl   1 69.055 1.1095  0.115   
+ emo  1 69.046 1.1174  0.125   
+ wc   1 69.071 1.0949  0.145   
+ cec  1 69.083 1.0832  0.155   
+ tk   1 69.068 1.0970  0.185   
+ nl   1 69.083 1.0836  0.205   
+ sl   1 69.114 1.0548  0.260   
+ ep   1 69.284 0.8967  0.865   
+ ek   1 69.321 0.8624  0.930   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step: sqrt(bray_dist_norep) ~ pH + tn 

      Df    AIC      F Pr(>F)  
+ tp   1 69.528 1.1400  0.060 .
+ emo  1 69.537 1.1319  0.090 .
+ den  1 69.541 1.1287  0.115  
+ cec  1 69.583 1.0912  0.175  
+ sl   1 69.602 1.0735  0.220  
+ fl   1 69.591 1.0840  0.225  
+ wc   1 69.597 1.0778  0.235  
+ tk   1 69.612 1.0651  0.285  
+ nl   1 69.701 0.9853  0.525  
+ ep   1 69.780 0.9145  0.805  
+ ek   1 69.827 0.8727  0.905  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.2 获得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)
pH.dist <- pH.dist[colnames(otu_phydist),colnames(otu_phydist)]
tn.dist <- tn.dist[colnames(otu_phydist),colnames(otu_phydist)]

4.3 可视化观察大致趋势

occor.r <- pH.dist
#occor.r <- tn.dist
data1 <- data.frame(cor = c(unlist(occor.r)),phydist=c(unlist(otu_phydist)))
data1 <- data1[!data1$cor==1,]
fac1<-cut(data1$phydist,2000)
data2<-cbind.data.frame(fac1,data1)
library(plyr)
data3<-ddply(data2,.(fac1),colwise(median))
plot(data3$phydist,data3$cor)

4.4 Mantel correlogram检验并输出结果

```r
occor.r[upper.tri(occor.r, diag=TRUE)] = NA
otu_phydist[upper.tri(otu_phydist, diag=TRUE)] = NA
man2 <- mantel.correlog(occor.r,otu_phydist)
write.table(man2$mantel.res,\~/Desktop/abiotic niche to test phylogenetic signal.txt\,sep = \\t\,quote=FALSE,row.names=TRUE)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### 5.基于栖息地距离的遗传发育信号检验
#### 5.1 计算栖息地距离并与遗传距离矩阵顺序保持一致

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuaGFiaXRhdF9icmF5IDwtIGFzLm1hdHJpeCh2ZWdkaXN0KHQobmlmaG90dV9ub3JlcCksIG1ldGhvZD1cImJyYXlcIikpXG5oYWJpdGF0X2JyYXkgPC0gaGFiaXRhdF9icmF5W2NvbG5hbWVzKG90dV9waHlkaXN0KSxjb2xuYW1lcyhvdHVfcGh5ZGlzdCldXG5gYGAifQ== -->

```r
habitat_bray <- as.matrix(vegdist(t(nifhotu_norep), method="bray"))
habitat_bray <- habitat_bray[colnames(otu_phydist),colnames(otu_phydist)]

5.2 可视化观察大致趋势

occor.r <- habitat_bray
#occor.r <- tn.dist
data1 <- data.frame(cor = c(unlist(occor.r)),phydist=c(unlist(otu_phydist)))
data1 <- data1[!data1$cor==1,]
fac1<-cut(data1$phydist,2000)
data2<-cbind.data.frame(fac1,data1)
library(plyr)
data3<-ddply(data2,.(fac1),colwise(median))
plot(data3$phydist,data3$cor)

5.3 Mantel correlogram检验并输出结果

```r
#该步骤耗时较长，建议输出变量保留结果。
occor.r[upper.tri(occor.r, diag=TRUE)] = NA
otu_phydist[upper.tri(otu_phydist, diag=TRUE)] = NA
man2 <- mantel.correlog(occor.r,otu_phydist)
write.table(man2$mantel.res,\~/Desktop/habitat similarity to test phylogenetic signal.txt\,sep = \\t\,quote=FALSE,row.names=TRUE)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### 6.基于物种关联矩阵的遗传发育信号检验
#### 6.1 计算物种关联矩阵并与遗传距离矩阵顺序保持一致

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuZmFzdF9zcGVhcm1hbl9mZHIgPC0gZnVuY3Rpb24ob3R1X25pY2hlKXtcbiAgciA8LSBjb3Iob3R1X25pY2hlLG1ldGhvZD1cInNwZWFybWFuXCIpXG4gIG4gPSBucm93KG90dV9uaWNoZSlcbiAgdCA8LSAociAqIHNxcnQobiAtIDIpKS9zcXJ0KDEgLSByXjIpXG4gIHAgPC0gLTIgKiBleHBtMShwdChhYnModCksIChuIC0gMiksIGxvZy5wID0gVFJVRSkpXG4gIHBbdXBwZXIudHJpKHAsIGRpYWcgPSBGQUxTRSldIDwtIHAuYWRqdXN0KHBbdXBwZXIudHJpKHAsZGlhZyA9IEZBTFNFKV0sIFwiZmRyXCIpXG4gIHBbbG93ZXIudHJpKHAsIGRpYWcgPSBGQUxTRSldID0gdChwKVtsb3dlci50cmkocCwgZGlhZyA9IEZBTFNFKV1cbiAgcmV0dXJuKGxpc3Qocj1yLHA9cCkpXG59XG5vY2NvciA8LSBmYXN0X3NwZWFybWFuX2ZkcihuaWZob3R1X25vcmVwKVxub2Njb3IuciA8LSBvY2NvciRyXG5vY2Nvci5yIDwtIG9jY29yLnJbY29sbmFtZXMob3R1X3BoeWRpc3QpLGNvbG5hbWVzKG90dV9waHlkaXN0KV1cbmBgYCJ9 -->

```r
fast_spearman_fdr <- function(otu_niche){
  r <- cor(otu_niche,method="spearman")
  n = nrow(otu_niche)
  t <- (r * sqrt(n - 2))/sqrt(1 - r^2)
  p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE))
  p[upper.tri(p, diag = FALSE)] <- p.adjust(p[upper.tri(p,diag = FALSE)], "fdr")
  p[lower.tri(p, diag = FALSE)] = t(p)[lower.tri(p, diag = FALSE)]
  return(list(r=r,p=p))
}
occor <- fast_spearman_fdr(nifhotu_norep)
occor.r <- occor$r
occor.r <- occor.r[colnames(otu_phydist),colnames(otu_phydist)]

6.2 可视化观察大致趋势

data1 <- data.frame(cor = c(unlist(occor.r)),phydist=c(unlist(otu_phydist)))
data1 <- data1[!data1$cor==1,]
fac1<-cut(data1$phydist,2000)
data2<-cbind.data.frame(fac1,data1)
library(plyr)
data3<-ddply(data2,.(fac1),colwise(median))
plot(data3$phydist,data3$cor)

6.3 Mantel correlogram检验并输出结果

```r
#该步骤耗时较长，建议输出变量保留结果。
occor.r[upper.tri(occor.r, diag=TRUE)] = NA
otu_phydist[upper.tri(otu_phydist, diag=TRUE)] = NA
man2 <- mantel.correlog(occor.r,otu_phydist)
write.table(man2$mantel.res,\~/Desktop/habitat similarity to test phylogenetic signal.txt\,sep = \\t\,quote=FALSE,row.names=TRUE)

```