Network properties of stochastic microbial assemblages

Raup Crick null model

raupcrick_nullmodel = function(spXsite,
                                reps = 999) {
  ## count number of sites and total species richness across all plots (gamma)
  n_sites <- nrow(spXsite)
  gamma <- ncol(spXsite)
  ##make the spXsite matrix into a new, pres/abs. matrix:
  ceiling(spXsite / max(spXsite)) -> spXsite.inc
  ##create an occurrence vector- used to give more weight to widely distributed species in the null model:
  occur <- apply(spXsite.inc, MARGIN = 2, FUN = sum)
  ##create an abundance vector- used to give more weight to abundant species in the second step of the null model:
  abundance <- apply(spXsite, MARGIN = 2, FUN = sum)
  ##make_null:
  nulllist <- NULL
  for (i in 1:reps){
    for (null.one in 1:nrow(spXsite)) {
      ##two empty null communities of size gamma:
      com1 <- rep(0, gamma)
      ##add observed number of species to com1, weighting by species occurrence frequencies:
      com1[sample(1:gamma,
                  sum(spXsite.inc[null.one, ]),
                  replace = FALSE,
                  prob = occur)] <- 1
      com1.samp.sp = sample(which(com1 > 0),
                            (sum(spXsite[null.one, ]) - sum(com1)),
                            replace = TRUE,
                            prob = abundance[which(com1 > 0)])
      com1.samp.sp = cbind(com1.samp.sp, 1)
      # head(com1.samp.sp);
      com1.sp.counts = as.data.frame(tapply(com1.samp.sp[, 2], com1.samp.sp[, 1], FUN = sum))
      colnames(com1.sp.counts) = 'counts'
      # head(com1.sp.counts);
      com1.sp.counts$sp = as.numeric(rownames(com1.sp.counts))
      # head(com1.sp.counts);
      com1[com1.sp.counts$sp] = com1[com1.sp.counts$sp] + com1.sp.counts$counts
      # com1;
      #sum(com1) - sum(spXsite[null.one,]); ## this should be zero if everything work properly
      rm('com1.samp.sp', 'com1.sp.counts')
      spXsite[null.one,] <- com1
    }
    nulllist[[i]] <- spXsite
  }
  return(nulllist)
}

Prepare OTU table

library(vegan)

## Warning: 程辑包'vegan'是用R版本4.1.3 来建造的

## 载入需要的程辑包：permute

## Warning: 程辑包'permute'是用R版本4.1.3 来建造的

## 载入需要的程辑包：lattice

## This is vegan 2.5-7

#读取OTU表
otu_biom2 <- read.delim("./ncm/otu table.txt", row.names=1)
#查看样品序列条数
summary(colSums(otu_biom2))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    8565   16699   19754   20431   22386   36187

otu_biom2  = otu_biom2[,colSums(otu_biom2) > 8500]
#将OTU表转置为行名为样品名称，列名为OTU名称的表
otu<-t(otu_biom2)
head(rowSums(otu))

##    t1    t2    t3    t4    t5    t6 
## 22027 11132 28759 16766 17243 16677

#对OTU表进行抽平
otu = as.data.frame(rrarefy(otu, 8500))
head(rowSums(otu))

##   t1   t2   t3   t4   t5   t6 
## 8500 8500 8500 8500 8500 8500

#去除序列条数较少的OTU
otu2<-otu[,log(colSums(otu)/sum(otu),10)>-4.5]
#去除出现频率低的微生物，即去除只在少数样品中出现的OTU
otu_niche<-otu2[,colSums(otu2>0)/nrow(otu2)>0.2]

Generate null otu table

nullmodel <- raupcrick_nullmodel(otu_niche,4)

Correlation method

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

Network construction of simulated null otu table

cor <- fast_spearman_fdr(nullmodel[[2]])

## Warning in cor(otu_niche, method = "spearman"): 标准差为零

#cor <- fast_spearman_fdr(otu_niche)
corr <- cor$r
corr[is.na(corr)] = 0
corr[cor$r < 0.6] = 0
corr[cor$p > 0.05] = 0
library(igraph)

## Warning: 程辑包'igraph'是用R版本4.1.3 来建造的

## 
## 载入程辑包：'igraph'

## The following object is masked from 'package:vegan':
## 
##     diversity

## The following object is masked from 'package:permute':
## 
##     permute

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

ng <- graph_from_adjacency_matrix(corr,mode="undirected",weighted=TRUE,diag = FALSE)
ng2 <- induced_subgraph(ng,colnames(corr)[colSums(corr)>1])

Network properties of null network

histogram(degree(ng2))

Network construction of observed otu table

cor <- fast_spearman_fdr(otu_niche)
#cor <- fast_spearman_fdr(otu_niche)
corr <- cor$r
corr[is.na(corr)] = 0
corr[cor$r < 0.6] = 0
corr[cor$p > 0.05] = 0
library(igraph)
ng <- graph_from_adjacency_matrix(corr,mode="undirected",weighted=TRUE,diag = FALSE)
ng2 <- induced_subgraph(ng,colnames(corr)[colSums(corr)>1])

Network properties of otu network

histogram(degree(ng2))