Load library
library(vegan)
library(Hmisc)
library(minpack.lm)
library(stats4)
Import test data
#读取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[1:5,1:5]
## 4479946 4479944 denovo5393 11544 623634
## t1 1 0 0 1 4
## t2 3 0 0 0 2
## t3 0 12 0 0 0
## t4 0 10 0 0 0
## t5 0 7 0 0 0
Get immigration probility \(m\) and \(R^2\)
#using Non-linear least squares (NLS) to calculate R2:
#spp: A community table with taxa as rows and samples as columns
spp<-otu
#remove other variable to save memory
#rm(otu2);rm(otu_niche);rm(otu_biom2);rm(otu);
#get the mean of abundance of each sample
N <- mean(apply(spp, 1, sum))
#get the mean of species relative abundance in metacommmunity
p.m <- apply(spp, 2, mean)
p.m <- p.m[p.m != 0]
#get the percentage of each species in metacommmunity
p <- p.m/N
#get the binary data of community abundance matrix
spp.bi <- 1*(spp>0)
#get the frequncy of species occurrence in metacommunity
freq <- apply(spp.bi, 2, mean)
freq <- freq[freq != 0]
#get a table record species percentage and occurrence frequency in metacommunity
C <- merge(p, freq, by=0)
#sort the table according to occurence frquency of each species
C <- C[order(C[,2]),]
C <- as.data.frame(C)
#delete rows containning zero
C.0 <- C[!(apply(C, 1, function(y) any(y == 0))),]
p <- C.0[,2]
freq <- C.0[,3]
names(p) <- C.0[,1]
names(freq) <- C.0[,1]
d = 1/N
##Fit model parameter m (or Nm) using Non-linear least squares (NLS)
m.fit <- nlsLM(freq ~ pbeta(d, N*m*p, N*m*(1-p), lower.tail=FALSE),start=list(m=0.1))
m.fit #get the m value
## Nonlinear regression model
## model: freq ~ pbeta(d, N * m * p, N * m * (1 - p), lower.tail = FALSE)
## data: parent.frame()
## m
## 0.3261
## residual sum-of-squares: 176.8
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 1.49e-08
m.ci <- confint(m.fit, 'm', level=0.95)
## Waiting for profiling to be done...
freq.pred <- pbeta(d, N*coef(m.fit)*p, N*coef(m.fit)*(1-p), lower.tail=FALSE)
pred.ci <- binconf(freq.pred*nrow(spp), nrow(spp), alpha=0.05, method="wilson", return.df=TRUE)
# get the R2 value
Rsqr <- 1 - (sum((freq - freq.pred)^2))/(sum((freq - mean(freq))^2))
Rsqr
## [1] 0.5413392
Data visulization
bacnlsALL <-data.frame(p,freq,freq.pred,pred.ci[,2:3])
inter.col<-rep('black',nrow(bacnlsALL))
inter.col[bacnlsALL$freq <= bacnlsALL$Lower]<-'#A52A2A'#define the color of below points
inter.col[bacnlsALL$freq >= bacnlsALL$Upper]<-'#29A6A6'#define the color of up points
library(grid)
grid.newpage()
pushViewport(viewport(h=0.6,w=0.6))
pushViewport(dataViewport(xData=range(log10(bacnlsALL$p)), yData=c(0,1.02),extension=c(0.02,0)))
grid.rect()
grid.points(log10(bacnlsALL$p), bacnlsALL$freq,pch=20,gp=gpar(col=inter.col,cex=0.7))
grid.yaxis()
grid.xaxis()
grid.lines(log10(bacnlsALL$p),bacnlsALL$freq.pred,gp=gpar(col='blue',lwd=2),default='native')
grid.lines(log10(bacnlsALL$p),bacnlsALL$Lower ,gp=gpar(col='blue',lwd=2,lty=2),default='native')
grid.lines(log10(bacnlsALL$p),bacnlsALL$Upper,gp=gpar(col='blue',lwd=2,lty=2),default='native')
grid.text(y=unit(0,'npc')-unit(2.5,'lines'),label='Mean Relative Abundance (log10)', gp=gpar(fontface=2))
grid.text(x=unit(0,'npc')-unit(3,'lines'),label='Frequency of Occurance',gp=gpar(fontface=2),rot=90)
#grid.text(x=unit(0,'npc')-unit(-1,'lines'), y=unit(0,'npc')-unit(-15,'lines'),label='Mean Relative Abundance (log)', gp=gpar(fontface=2))
#grid.text(round(coef(m.fit)*N),x=unit(0,'npc')-unit(-5,'lines'), y=unit(0,'npc')-unit(-15,'lines'),gp=gpar(fontface=2))
#grid.text(label = "Nm=",x=unit(0,'npc')-unit(-3,'lines'), y=unit(0,'npc')-unit(-15,'lines'),gp=gpar(fontface=2))
#grid.text(round(Rsqr,2),x=unit(0,'npc')-unit(-5,'lines'), y=unit(0,'npc')-unit(-16,'lines'),gp=gpar(fontface=2))
#grid.text(label = "Rsqr=",x=unit(0,'npc')-unit(-3,'lines'), y=unit(0,'npc')-unit(-16,'lines'),gp=gpar(fontface=2))
draw.text <- function(just, i, j) {
grid.text(paste("Rsqr=",round(Rsqr,3),"\n","Nm=",round(coef(m.fit)*N)), x=x[j], y=y[i], just=just)
#grid.text(deparse(substitute(just)), x=x[j], y=y[i] + unit(2, "lines"),
# gp=gpar(col="grey", fontsize=8))
}
x <- unit(1:4/5, "npc")
y <- unit(1:4/5, "npc")
draw.text(c("centre", "bottom"), 4, 1)
