####LIBRARIES####
library("vegan")
library("adespatial")
library("ade4")
library("lme4)
library("ggplot2")
library("hillR")
library("FactoMineR")
library("metabaR")
library("MASS")

#####LEAF TRAITS VARIATIONS BETWEEN POSITION AND TREE SPECIES####
tmp_leaf=dplyr::select(tt,carre,arbre,binom,position,feuille,samples, sample_id,position_leave=epi_endo,x=x.utm,y=y.utm,chloro=Chloro.microg.cm2,thickness=THICKNESS.microm,ldmc=LDMC.mg.g,leaf.area=Leaf.area.cm2,lma=LMA.g.m2,ccm=CCm.gGluc.gMS,lwc=leaf.water.content,C=C.gC.gMS,N=N.gN.gMS,Ca, Ash, Cu,Fe, Mg,Mn,C.N,P,K,Na, Zn,SD=stomatal.density)

##test normality of variables
shapiro.test(tmp_leaf$chloro) # p-value = 0.2831 normal
shapiro.test(tmp_leaf$thickness)#non normal
## to do for each variable and transform accordingly with boxcox

## example with thickness - do the same for each leaf trait
bc<-boxcox(tmp_leaf$thickness~tmp_leaf$position)
(lambda <- bc$x[which.max(bc$y)])
new_model <- lm(((tmp_leaf$thickness^lambda-1)/lambda) ~ tmp_leaf$position)
qqnorm(new_model$residuals)
qqline(new_model$residuals)
norm_thickness<-(tmp_leaf$thickness^lambda-1)/lambda

##effects of position X tree species on thickness
lm_thickness<-aov(norm_thickness~position*binom, data=tmp_norm)
summary(lm_thickness)

g8<-ggplot(tmp_leaf_short,aes(x=binom, y=thickness, fill=position))+
  geom_boxplot(outlier.colour = "white")+
  scale_fill_manual(values=c("white","azure3","azure4"))+
  scale_x_discrete(breaks=c("Eperua_falcata", "Macrolobium_bifolium","Tetragastris_sp"),labels=c("E.f","M.b", "T.sp") )+
  theme_bw()+
  theme(panel.grid = element_blank(),axis.title.x = element_blank(),legend.position = "none",axis.title.y = element_text(size=12), axis.text.x = element_text(size =12))+
  ylab("Thickness")+
  annotate("text",x=2,y=290,label="ANOVA,Host ***",fontface="bold")+
  annotate("text",x=2,y=285,label="Position *",fontface="bold")+
  annotate("text",x=2,y=280,label="Host x Position NS",fontface="bold")

####FUNGI####

#FUNGI FILES
# load OTU and samples data and build metabaR files (see URL: https://github.com/metabaRfactory/metabaR)
# compile infos
d=read.delim(file="fungi_97",header=T,sep="\t")
d2=read.delim(file = "fungi_blast_Unite_V2-1.txt",header=T,sep="\t")
d3=read.delim(file="fungi_embl_ecotag_clust97_true_k.tab",header=T)
tmp<-data.frame(right_join(d,d2,by="id"))
data_f<-data.frame(left_join(tmp,d3[,c(1,417)],by="id"))

##reads
tmp2=droplevels.data.frame(data_f[,c(1,14:424)])
tmp2[is.na(tmp2)] <- 0
write.table(tmp2,file="tmp2.txt",row.names = F)
reads=t(tmp2[,-1])
colnames(reads)=tmp2$id
rownames(reads)=colnames(tmp2[,-1])

##check NA values by rows in data.frame
apply(reads, 1, function(x) any(is.na(x)))
reads<-reads[!(rownames(reads)=="f12_r16"),]
reads=as.matrix(reads)

##motus
motus=data_f[,c(1,5,13,425:433)]
rownames(motus)=motus$id

##samples
tmp5=read.delim("db_field_vertige_corrige.txt",header=T)
tmp6<-tmp5[which(tmp5$sample_id %in% rownames(reads)),]
samples=dplyr::select(tmp6,carre,arbre,binom,position,feuille,samples, sample_id,position_leave=epi_endo,x=x.utm,y=y.utm,chloro=Chloro.microg.cm2,thickness=THICKNESS.microm,ldmc=LDMC.mg.g,leaf.area=Leaf.area.cm2,lma=LMA.g.m2,ccm=CCm.gGluc.gMS,lwc=leaf.water.content,C=C.gC.gMS,N=N.gN.gMS,Ca, Ash, Cu,Fe, Mg,Mn,C.N,P,K,Na, Zn,SD=stomatal.density)
rownames(samples)<-samples$sample_id
samples<-samples[!(samples$sample_id=="f12_r16"),]

##pcrs
tmp3=read.delim2("ngsfilter_its.txt",header=T)
pcrs=dplyr::filter(tmp3, tmp3$sample_id %in% rownames(reads))
rownames(pcrs)=pcrs$sample_id

##metabarlist
fungi=metabarlist_generator(
  reads = reads,
  motus = motus,
  pcrs = pcrs,
  samples = samples)

##tagjump contamination
fungi_clean<-tagjumpslayer(fungi,method="cut",0.03)
fungi_clean_cont<-contaslayer(fungi_clean)

## Distribution of the most abundant contaminant MOTU in the PCR plate design
contaminants <- rownames(fungi_clean_cont$motus)[fungi_clean_cont$motus$not_a_contaminant == FALSE]
max.conta <- contaminants[which.max(fungi_clean_cont$motus[contaminants, "COUNT"])]

p <- ggpcrplate(fungi_clean_cont,
  legend_title = "# reads",
  FUN = function(m) {
    m$reads[, max.conta]
  }
)

p + scale_size(limits = c(1, max(fungi_clean_cont$reads[, max.conta]))) +
  ggtitle("Distribution of the most abundant contaminant")

## remove contaminants in datasets
fungi_final<-subset_metabarlist(fungi_clean_cont,table = "motus",indices = (fungi_clean_cont$motus$not_a_contaminant==TRUE))

## keep only samples for analyses
tp<-subset_metabarlist(fungi_final,table="samples",indices = (!(fungi_final$samples$samples %in% "control")))
#keep only HMBS as position (L=senescent)
fungi_final_samples<-subset_metabarlist(tp,table="samples",indices=(!(tp$samples$position %in% "L")))
#remove samples with 0 data
fungi_final_samples<-subset_metabarlist(fungi_final_samples,table="reads",indices=(rowSums(fungi_final_samples$reads) !=0))

##split endo and epiphyte
fungi_final_samples_endo<-subset_metabarlist(fungi_final_samples,table="samples",indices=((fungi_final_samples$samples$position_leave %in% "endo")))
fungi_final_samples_epi<-subset_metabarlist(fungi_final_samples,table="samples",indices=((fungi_final_samples$samples$position_leave %in% "epi")))


summary_metabarlist(fungi_final_samples)
summary_metabarlist(fungi_final_samples_endo)
summary_metabarlist(fungi_final_samples_epi)



#RELATIVE CONTRIBUTION OF ENVIRONMENT, LEAF TRAITS AND SPACE ON COMMUNITIES
##compute MEM (spatial variability)
##ENDOPHYTES
tmp_leaf_endo<-fungi_final_samples_endo$samples

mem_endo<-dbmem(tmp_leaf_endo[,9:10],MEM.autocor = "non-null")
##Compute and test associated Moran's I values
## Eigenvalues are proportional to Moran's 
test <- moran.randtest(mem_endo, nrepet = 999) # 7 first are significant (positive eigenvalues)
plot(test$obs, xlab="MEM rank", ylab="Moran's I")
abline(h=-1/(nrow(tmp_leaf_endo[,9:10]) - 1), col="red")


## compute alpha and betadiversities
##endophytes - rarefaction curves - do the same for fungal epiphytes communities
alpha_fungi<-hill_rarefaction(fungi_final_samples,nboot=20,nsteps=10)

p<-gghill_rarefaction(alpha_fungi, group=fungi_final_samples$samples$position_leave)

##endophytes - alpha diversity (Hill numbers, with q=1 for Shannon) - do the same with fungal epiphytes communities
alpha_fungi_endo<-hill_taxa(fungi_final_samples_endo$reads,q=1)

##endophytes - betadiversity (Hill numbers, with q=1 , equivalent to Morisita-Horn dissimilarity matrix) -do the same with fungal epiphytes communities
fun_dist_endo<-hill_taxa_parti_pairwise(fungi_final_samples_endo$reads,q=1,output="matrix", pairs = "full")

##endophytes - effects on alphadiversity with a linear model and stepwise variable selection - do the same with fungal epiphytes communities
fun_lm<-lm(norm_Shannon~., env2)
fun_lm.step<-stepAIC(fun_lm,direction="both")

##endophytes - effects on betadiversity with a stepwise model for constrained ordination methods
fun_rda.vasc.0 <- dbrda (fun_dist_endo$TD_beta~ 1, data = env2) # model containing only species matrix and intercept
fun_rda.vasc.all <- dbrda (fun_dist_endo$TD_beta ~ ., data = env2[,-10]) # model including all variables from matrix chem1 
fun_sel.osR2 <- ordiR2step (fun_rda.vasc.0, scope = formula (fun_rda.vasc.all), direction = c("both", "backward", "forward"))
fun_sel.osR2$anova

#RELATIVE CONTRIBUTION OF DETERMINISTIC AND STOCHASTIC PROCESSES IN THE ASSEMBLY : NDV - R code proposed by Luan et al 2020
# to do for all Fungi, and endophytes and epiphytes communities separately
source("MetacommunityDynamicsFctsOikos.R")
source("PANullDevFctsOikos.R")

##all fungi
### Prepare and calculate abundance beta-null deviation metric
## Adjusted from Stegen et al 2012 GEB
comm.t=t(fungi_final_samples$reads)
bbs.sp.site <- comm.t
patches=nrow(bbs.sp.site)
rand <- 999
null.alphas <- matrix(NA, ncol(comm.t), rand)
null.alpha <- matrix(NA, ncol(comm.t), rand)
expected_beta <- matrix(NA, 1, rand)
null.gamma <- matrix(NA, 1, rand)
null.alpha.comp <- numeric()
bucket_bray_res <- matrix(NA, patches, rand)

bbs.sp.site = ceiling(bbs.sp.site/max(bbs.sp.site)) 
mean.alpha = sum(bbs.sp.site)/nrow(bbs.sp.site) #mean.alpha
gamma <- ncol(bbs.sp.site) #gamma
obs_beta <- 1-mean.alpha/gamma
obs_beta_all <- 1-rowSums(bbs.sp.site)/gamma

for (randomize in 1:rand) {  
  null.dist = comm.t
  for (species in 1:ncol(null.dist)) {
    tot.abund = sum(null.dist[,species])
    null.dist[,species] = 0
    for (individual in 1:tot.abund) {
      sampled.site = sample(c(1:nrow(bbs.sp.site)), 1)
      null.dist[sampled.site, species] = null.dist[sampled.site, species] + 1
    }
  }
    ##Calculate null deviation for null patches and store
  null.alphas[,randomize] <- apply(null.dist, 2, function(x){sum(ifelse(x > 0, 1, 0))})
  null.gamma[1, randomize] <- sum(ifelse(rowSums(null.dist)>0, 1, 0))
  expected_beta[1, randomize] <- 1 - mean(null.alphas[,randomize]/null.gamma[,randomize])
  null.alpha <- mean(null.alphas[,randomize])
  null.alpha.comp <- c(null.alpha.comp, null.alpha)

  bucket_bray <- as.matrix(vegdist(null.dist, "bray"))
    diag(bucket_bray) <- NA
    bucket_bray_res[,randomize] <- apply(bucket_bray, 2, FUN="mean", na.rm=TRUE)
} ## end randomize loop

## Calculate beta-diversity for obs metacommunity
beta_comm_abund <- vegdist(comm.t, "bray")
res_beta_comm_abund <- as.matrix(as.dist(beta_comm_abund))
diag(res_beta_comm_abund) <- NA 
#output beta diversity (Bray)
beta_div_abund_stoch <- apply(res_beta_comm_abund, 2, FUN="mean", na.rm=TRUE)

# output abundance beta-null deviation
bray_abund_null_dev_fun <- beta_div_abund_stoch - mean(bucket_bray_res)
bray_abund_null_dev_fun

#ESTIMATING NICHE BREADTH WITH OMI (R code Luan and al 2020)
##to do for all fungi, and endophytes and epiphytes separately
t<-fungi_final_samples$samples
tbis <- t[,c(11:31)] %>% #replace NA by mean for each traits
     mutate(chloro=replace_na(chloro, mean(chloro, na.rm=TRUE)),
            thickness=replace_na(thickness,179), #moyenne calculee a part et garder sans decimale
            ldmc=replace_na(ldmc, mean(ldmc, na.rm=TRUE)),
            leaf.area=replace_na(leaf.area, mean(leaf.area, na.rm=TRUE)),
            lma=replace_na(lma, mean(lma, na.rm=TRUE)),
            ccm=replace_na(ccm, mean(ccm, na.rm=TRUE)),
            lwc=replace_na(lwc, mean(lwc, na.rm=TRUE)),
            C=replace_na(C, mean(C, na.rm=TRUE)),
            N=replace_na(N, mean(N, na.rm=TRUE)),
            Ca=replace_na(Ca, mean(Cu, na.rm=TRUE)),
            Ash=replace_na(Ash, mean(Ash, na.rm=TRUE)),
            Cu=replace_na(Cu, mean(Ca, na.rm=TRUE)),
            Fe=replace_na(Fe, mean(Fe, na.rm=TRUE)),
            Mg=replace_na(Mg, mean(Mg, na.rm=TRUE)),
            Mn=replace_na(Mn, mean(Mn, na.rm=TRUE)),
            C.N=replace_na(C.N, mean(C.N, na.rm=TRUE)),
            P=replace_na(P, mean(P, na.rm=TRUE)),
            K=replace_na(K, mean(K, na.rm=TRUE)),
            Na=replace_na(Na, mean(Na, na.rm=TRUE)),
            Zn=replace_na(Zn, mean(Zn, na.rm=TRUE)),
            SD=replace_na(SD, mean(SD, na.rm=TRUE)))
tbis<-data.frame(tbis,t[,c(3,4)])
tbis$position<-as.factor(tbis$position)
tbis$binom<-as.factor(tbis$binom)

OTU_all=data.frame(fungi_final_samples$reads)
dudi1_all <- dudi.mix(tbis, scannf = FALSE, nf = 3)
nic1_fun_all <- niche(dudi1_all, OTU_all, scann = FALSE)
niov_fun_all=niche.param(nic1_fun_all)
niche_fun_all=rtest(nic1_fun_all,19)

####BACTERIA####
#use the same script to analyze endophytic and epiphytic bacterial communities