rm(list=ls())
t0=Sys.time() # to calculate time cost

####################
# 1 # set folder paths and file names, please change according to the folder paths and file names in your computer.
####################
# the folder saving the input files
wd="/Data/"

# the OTU table file (Tab delimited txt file)
com.file="100WSc.OTUtable.csv"

# the phylogenetic tree file
tree.file="100WSc.Tree.nwk"

# the classification (taxonomy) information
clas.file="100WSc.Classif.QIIME2.Silva138.csv"

# the treatment informaiton table
treat.file="/Code/StressLevel/grouping.stressA.msi.csv"

# the environmental varialbes
env.file="/Code/CommunityAssembly/100WSc.EnvFactor.csv"

# estimated relative abundance in metacommunity
metaab.file="100WSc.estimateLog.meta.abs.csv"

# the folder to save the output. please change to a new folder even if you are just testing the example data.
save.wd="/Code/CommunityAssembly"

####################
# 2 # key parameter setting
####################
prefix="Q100W"  # prefix of the output file names. usually use a project ID.
rand.time=1000  # randomization time, 1000 is usually enough. For example test, you may set as 100 or less to save time.
nworker=25 # nworker is thread number for parallel computing, which depends on the CPU core number of your computer.
memory.G=100 # to set the memory size as you need (but should be less than the available space in your hard disk), so that calculation of large tree will not be limited by physical memory. unit is Gb.

####################
# 3 # load R packages and data
####################
library(iCAMP)
library(ape)
source("/Code/handytool.r")
setwd(wd)
comm=t(lazyopen(com.file))
tree=lazyopen(tree.file)
clas=lazyopen(clas.file)
treat=lazyopen(treat.file)
env=lazyopen(env.file)
metaab=t(lazyopen(metaab.file))
####################
# 4 # match sample IDs in OTU table and treatment information table
####################
sampid.check=match.name(rn.list=list(comm=comm,treat=treat,env=env))
treat=sampid.check$treat
comm=sampid.check$comm
comm=comm[,colSums(comm)>0,drop=FALSE]
env=sampid.check$env

####################
# 5 # match OTU IDs in OTU table and tree file
####################
spid.check=match.name(cn.list=list(comm=comm),rn.list=list(clas=clas),tree.list=list(tree=tree))
comm=spid.check$comm
clas=spid.check$clas
tree=spid.check$tree

# 6 # calculate pairwise phylogenetic distance matrix.
# since microbial community data usually has a large number of species (OTUs or ASVs), we use "big.matrix" in R package "bigmemory" to handle the large phylogenetic distance matrix. 
setwd(save.wd)
if(!file.exists("pd.desc")) 
{
  pd.big=iCAMP::pdist.big(tree = tree, wd=save.wd, nworker = nworker, memory.G = memory.G)
  # output files:
  # path.rda: a R object to list all the nodes and  edge lengthes from root to every tip. saved in R data format. an intermediate output when claculating phylogenetic distance matrix.
  # pd.bin: BIN file (backingfile) generated by function big.matrix in R package bigmemory. This is the big matrix storing pairwise phylogenetic distance values. By using this bigmemory format file, we will not need memory but hard disk when calling big matrix for calculation.
  # pd.desc: the DESC file (descriptorfile) to hold the backingfile (pd.bin) description.
  # pd.taxon.name.csv: comma delimited csv file storing the IDs of tree tips (OTUs), serving as the row/column names of the big phylogenetic distance matrix.
}else{
  # if you already calculated the phylogenetic distance matrix in a previous run
  pd.big=list()
  pd.big$tip.label=read.csv(paste0(save.wd,"/pd.taxon.name.csv"),row.names = 1,stringsAsFactors = FALSE)[,1]
  pd.big$pd.wd=save.wd
  pd.big$pd.file="pd.desc"
  pd.big$pd.name.file="pd.taxon.name.csv"
}

####################
# 6 # Phylogenetic Signal Test
####################
# 6.1 # env factor transformation and standardization
################
envin=env
env=envin[,3:ncol(envin)]
lnx<-function(v)
{
  minv=min(v)
  if(minv<=0){minv=min(density(v)$x);v=v-minv}
  log(v)
}
envln=sapply(1:ncol(env),function(i){lnx(env[,i])})
colnames(envln)=colnames(env)
rownames(envln)=rownames(env)
envln[,which(colnames(envln)=="pH")]=env[,"pH"]


################
# 6.2 # niche difference
################
envi="pH" # test representative factors
iwd(save.wd)
library(iCAMP)
envim=envln[,envi,drop=FALSE]
nichedi=iCAMP::dniche(env = envim,comm = comm,method = "niche.value",
                      nworker = nworker,out.dist=FALSE,bigmemo=TRUE,
                      nd.wd=save.wd)

################
# 6.3 # root the tree
################
treert=iCAMP::midpoint.root.big(tree=tree,pd.desc = pd.big$pd.file,pd.spname = pd.big$tip.label,
                                pd.wd = pd.big$pd.wd,nworker = nworker)

tree.rt=treert$tree
(mpd=treert$max.pd)

################
# 6.4 # phylogenetic signal test
################
sum(nichedi$names!=pd.big$tip.label) # should be zero
t1=Sys.time()
mcrl1=iCAMP::big.mantel.correlog(x.desc=pd.big$pd.file,x.name=pd.big$tip.label,x.wd=pd.big$pd.wd,
                                 y.desc=nichedi$nd[[1]],y.name=nichedi$names,y.wd=nichedi$nd.wd,
                                 nworker=nworker,interval=0.02,break.pts=NULL,
                                 permutations=999,strata=NULL,
                                 diag.in=FALSE,temp.wd=save.wd,stepw=10^6,padjust.method="fdr")

save(mcrl1,file=paste0(envi,".MantelCorrelog.niche.value.rda"))

####################
# 7 # iCAMP analysis
####################
bin.size.limit = 48 # For real data, usually use a proper number according to phylogenetic signal test or try some settings then choose the reasonable stochasticity level. our experience is 12, or 24, or 48. but for this example dataset which is too small, have to use 5.
sig.index="SES.RC" # see other options in help document of icamp.big.
icres=iCAMP::icamp.big(comm=comm, pd.desc = pd.big$pd.file, pd.spname=pd.big$tip.label,
                       pd.wd = pd.big$pd.wd, rand = rand.time, tree=tree,
                       prefix = prefix, ds = 0.2, pd.cut = NA, sp.check = TRUE,
                       phylo.rand.scale = "within.bin", taxa.rand.scale = "across.all",
                       phylo.metric = "bMPD", sig.index=sig.index, bin.size.limit = bin.size.limit, 
                       nworker = nworker, memory.G = memory.G, rtree.save = FALSE, detail.save = TRUE, 
                       qp.save = FALSE, detail.null = FALSE, ignore.zero = TRUE, output.wd = save.wd, 
                       correct.special = TRUE, unit.sum = rowSums(comm), special.method = "Simple",
                       ses.cut = 1.96, rc.cut = 0.95, conf.cut=0.975, omit.option = "no",meta.ab = metaab[1,])
save.file(icres$bNRIiRCa,prefix=prefix,filename = "ProcessImp.Pairwise.bNRIi.RCa")

####################
# 8 # iCAMP bin level statistics
####################
treats=data.frame(CountAll="All",treats,stringsAsFactors = FALSE)
icbin=icamp.bins(icamp.detail = icres$detail,treat = treats,
                 clas = clas,boot = TRUE, rand.time = 1000)
save(icbin,file = paste0(prefix,".iCAMPBinSummaryDetail.rda"))

Ptk=icbin$Ptk
BPtk=icbin$BPtk
BRPtk=icbin$BRPtk
clas.bin=icbin$Class.Bin
bin.clas=icbin$Bin.TopClass
Binwt=icbin$Binwt
save.file(Ptk,prefix=prefix,filename = "ProcessImpBin.Ptk")
save.file(BPtk,prefix=prefix,filename = "BinContrProcess.BPtk")
save.file(BRPtk,prefix=prefix,filename = "BinRelContrProcess.BPtk")
save.file(clas.bin,prefix=prefix,filename = "OTUBinClass")
save.file(bin.clas,prefix=prefix,filename = "BinTopOTU")
save.file(Binwt,prefix=prefix,filename = "BinRelAbundance")

####################
# 9 # Other approach: QPEN (quantifying community assembly processes based on entire-community null model analysis)
####################
# 9.1 # QPEN calculation
qpout=iCAMP::qpen.cm(comm=comm,pd=pd.big$pd.file,pd.big.wd=pd.big$pd.wd,
                     pd.big.spname=pd.big$tip.label,meta.ab=metaab,ab.weight=TRUE,
                     rand.time=rand.time, nworker=nworker,project=prefix,
                     wd=save.wd, save.bNTIRC=TRUE)
# 9.2 # significance test
qptest=iCAMP::qpen.test(qpen.result = qpout,treat = treat,rand.time = rand.time,
                        between.group = TRUE,out.detail=TRUE,silent=FALSE)
save(qptest,file = paste0(prefix,".cM.QPEN.boot.detail.rda"))

####################
# 10 # Neutral taxa percentage
####################
snmout=iCAMP::snm.comm(comm = comm, treat = treat, meta.com = metaab,
                       rand=rand.time, alpha=0.05)
save(snmout, file = paste0(prefix,".cM.S7n.Rep.NPdetail.rda"))

####################
# 14 # tNST and pNST (taxonomic and phylogenetic normalized stochasticity ratio)
####################
library(NST)

# 14.1 # tNST
tnstout=NST::tNST(comm=comm, group=treat, meta.com=metaab, dist.method="jaccard",
                  abundance.weighted=TRUE, rand=rand.time,  
                  nworker=nworker, null.model="PF", output.rand = TRUE,
                  SES = TRUE, RC = TRUE)
tnst.bt=NST::nst.boot(nst.result=tnstout, group=treat,
                      rand=rand.time, nworker=nworker)
save(tnst.bt, file = paste0(prefix,".tNSTboot.rda"))

# 14.2 # pNST
pnstout=NST::pNST(comm=comm, pd.desc=pd.big$pd.file, pd.wd=pd.big$pd.wd, 
                  pd.spname=pd.big$tip.label, group=treat, abundance.weighted=TRUE,
                  rand=rand.time, phylo.shuffle=TRUE, nworker=nworker,
                  output.rand = TRUE, SES=FALSE, RC=FALSE)
save(pnstout, file = paste0(prefix,".pNST.rda"))

pnst.bt=NST::nst.boot(nst.result=pnstout, group=treat,
                      rand=rand.time, nworker=nworker)
save(pnst.bt, file = paste0(prefix,".pNST.boot.rda"))

# End #