phylogenetic properties of cluster in co-occurrence network

Load R package

library(igraph)
library(vegan)
library(picante)
library(ape)

Load test data

igraph <- read_graph("~/Desktop/htl_data_analysis/R code/netwok and env/mygraph.graphml",format="graphml")
tree <- read.tree("~/Desktop/htl_data_analysis/R code/netwok and env/tree")
nifh_otu2 <- read.delim("~/Desktop/htl_data_analysis/R code/netwok and env/nifh_otu.txt", row.names=1)

Modular division and get the top 4 cluster ID

wtc <- cluster_louvain(igraph,NA)
modularity(wtc)

## [1] 0.610081

V(igraph)$module<-membership(wtc)
modtop4 <- order(tabulate(membership(wtc)),decreasing=TRUE)[1:4] 

Build new OTU table for NTI analysis

rlen <- length(modtop4)+1
#clen <- sum(tabulate(membership(wtc))[modtop4])
clen <- ncol(nifh_otu2)
#otu_new <- matrix(rep(0,rlen*clen),rlen,clen,dimnames=list(modtop4,V(igraph)$name[V(igraph)$module%in%modtop4]))
otu_new <- matrix(rep(0,rlen*clen),rlen,clen,dimnames=list(c(modtop4,"other"),colnames(nifh_otu2)))
for (i in 1:rlen){
  if (i < rlen){
    otu_new[i,] = as.integer(colnames(otu_new) %in% V(igraph)$name[V(igraph)$module==modtop4[i]])
  }
  else{
    otu_new[i,] = as.integer(!(colnames(otu_new) %in% V(igraph)$name[V(igraph)$module %in% modtop4]))
  }
}
otu_new[1:5,1:5]

##       OTU_1 OTU_10056 OTU_10057 OTU_10058 OTU_10059
## 6         0         0         0         0         0
## 1         1         0         0         0         0
## 2         0         1         0         0         1
## 7         0         0         0         0         0
## other     0         0         1         1         0

Tidy the tree and OTU table data for NTI analysis

tree$tip.label <- stringr::str_replace_all(tree$tip.label,"'","")

prune_tree<-prune.sample(nifh_otu2,tree)
tip<-prune_tree$tip.label
coln<-colnames(otu_new)
m<-NULL
print("These OTUs don't exist in phylogenetic tree")

## [1] "These OTUs don't exist in phylogenetic tree"

for(i in 1:length(coln)){
  if(!coln[i]%in%tip){
    print(coln[i])
    m<-cbind(m,coln[i])
  }
}

## [1] "OTU_10065"
## [1] "OTU_10179"
## [1] "OTU_10329"
## [1] "OTU_113977"
## [1] "OTU_122"
## [1] "OTU_159"
## [1] "OTU_22258"
## [1] "OTU_24594"
## [1] "OTU_26066"
## [1] "OTU_267595"
## [1] "OTU_35216"
## [1] "OTU_38790"
## [1] "OTU_43486"
## [1] "OTU_45378"
## [1] "OTU_53268"
## [1] "OTU_57855"
## [1] "OTU_87080"

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

Get the NTI result of each cluster

otu_phydist <- cophenetic(prune_tree)
mntd_wei<-ses.mntd(otu_new2, otu_phydist, null.model="taxa.labels", abundance.weighted=FALSE, runs=999)
mntd_wei

##       ntaxa   mntd.obs mntd.rand.mean mntd.rand.sd mntd.obs.rank  mntd.obs.z
## 6        84 0.22320524     0.23262750  0.036416803           423 -0.25873398
## 1        79 0.23558329     0.23586409  0.037435900           519 -0.00750086
## 2        76 0.10221605     0.24237875  0.038114871             1 -3.67737571
## 7        78 0.28986359     0.23677619  0.037524842           911  1.41472682
## other  1661 0.05257889     0.05291062  0.001674742           428 -0.19807948
##       mntd.obs.p runs
## 6          0.423  999
## 1          0.519  999
## 2          0.001  999
## 7          0.911  999
## other      0.428  999