wd="/Code/AlphaBetaDiversity"
comf="/Data/100WSc.OTUtable.csv"
treef="/Data/100WSc.Tree.nwk"

##################
# 1 # load data
##################
source("/Code/handytool.r")
comm=t(lazyopen(comf))
tree=lazyopen(treef)

##################
# 2 # Bray Curtis 
##################
library(vegan)
BC=vegan::vegdist(comm,method = "bray",binary = FALSE)
save.file(as.matrix(BC),filename = "Q100WSc.Bray.wt")


##################
# 3 # phylogenetic distance 
##################
library(iCAMP)
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"
}

##################
# 4 # beta MNTD
##################

bmntd=iCAMP::bmntd.big(comm=comm, pd.desc = pd.big$pd.file,
                       pd.spname = pd.big$tip.label, pd.wd = pd.big$pd.wd,
                       abundance.weighted = TRUE)
save.file(as.matrix(bmntd),filename = "Q100WSc.bMNTD.wt")


##################
# 5 # PERMANOVA
##################
msif="/Code/StressLevel/100W.MSI.csv"
msi=lazyopen(msif)
betaif="/Code/AlphaBetaDiversity/Q100WSc.Bray.wt.csv"
betai=lazyopen(betaif)
sampc=match.name(rn.list = list(treat,msi,logmsi),both.list = list(betai)) # should be "all match very well"
library(vegan)
(permaov=adonis2(as.dist(betai)~MSI,data = msi,permutations = 1000))


betaif="/Code/AlphaBetaDiversity/Q100WSc.bMNTD.wt.csv"
betai=lazyopen(betaif)
sampc=match.name(rn.list = list(treat,msi,logmsi),both.list = list(betai)) # should be "all match very well"
library(vegan)
(permaov=adonis2(as.dist(betai)~MSI,data = msi,permutations = 1000))

# End