BatchAdjust_clean.Rmd

---
title: "BatchAdjust"
output: html_document
editor_options: 
  chunk_output_type: console
---
```{r load libraries}
#Load libraries
#install.packages("easypackages")
#if (!requireNamespace("BiocManager", quietly = TRUE))
 #   install.packages("BiocManager")

library(easypackages) #can also install multiples with packages()
x<-c("ggplot2", "reshape2", "broom", "dplyr", "tidyverse", "GUniFrac", "phangorn", "doParallel", "clustsig","scales", "grid", "vegan", "survival",  "data.table","ape", "Biostrings", "RColorBrewer", "devtools","ampvis2", "metacoder", "VennDiagram", "vegan", "limma","taxa",  "readr", "stringr", "phyloseq", "DESeq2", "microbiomeSeq", "morpheus", "htmltools")
libraries(x)
rm(x)
#BiocManager::install("MMUPHin")
library(MMUPHin); library(Rglpk); library(igraph)
#load("BatchAdjust_clean.RData")
```


#loading paths and defining data
```{r paths and input}
input0=NA
input0$path="~/Documents/Collaborations/HelminthsMetaAnalysis/input-txts/output-merged/"
input0$counts=paste0(input0$path, "mergedCounts.txt")
input0$meta=paste0(input0$path, "mergedMETA_04122021.txt")
input0$tax=paste0(input0$path, "mergedTax.txt")
input0$amp=paste0(input0$path, "mergedTableWlineage.txt")

output0=NA
output0$path="~/Documents/Collaborations/HelminthsMetaAnalysis/input-txts/output-merged/"
output0$flt=paste0(input0$path,"Filtering/")
dir.create(output0$flt)
print(paste("creating output at:", output0$flt))
```

### Loading data
```{r load data}
data0=NA
#load counts table (USED FULL table, not the filtered _flt4 one)
data0$counts_all=read.table(input0$counts, header=T, check.names=F, row.names=1) ; head(data0$counts)
data0$counts_all<-data0$counts_all[ , order(names(data0$counts_all))] #order by column names
dim(data0$counts_all)


#metadata
data0$meta_all=read.table(input0$meta, header=T, check.names = F, sep='\t', row.names=1)
head(data0$meta_all); dim(data0$meta_all)
data0$meta_all <- data0$meta_all[ order(row.names(data0$meta_all)), ];  head(data0$meta_all)
dim(data0$meta_all) #Order by row names:

##remove 3 samples that have no Age or Gender info:
data0$meta_all=subset(data0$meta_all, data0$meta_all$AgeCategory!="NA")
data0$meta_all=subset(data0$meta_all, data0$meta$Gender!="NA")
dim(data0$meta_all)

#compare metadata sample names to counts sample names and remove problematic:
setdiff(colnames(data0$counts_all), row.names(data0$meta_all)) #Samples that will not be analyzed (to be removed from counts)
setdiff(row.names(data0$meta_all), colnames(data0$counts_all)) #samples that cannot be analyzed (to be removed from meta)
   #some samples in metadata did not have proper seqdata (were not processed). 
        #Lee_2019 samples had no seq data and need removal from metadata
 data0$meta_all=subset(data0$meta_all, StudyID!="Lee_2019")
 dim(data0$meta_all)
  #and 1 sample in metadata had no sequence files: SRR8356208
 data0$meta_all= data0$meta_all[!(row.names(data0$meta_all) %in% "SRR8356208"), ]
 dim(data0$meta_all)

#extract from "counts" only the relevant samples -> those from "meta"
data0$counts_all=data0$counts_all[, row.names(data0$meta_all) ] #seems im missing 115 Rubel2020 samples= fixed
dim(data0$counts_all)


#Now load the taxonomy (taxonomy_ NOT filtered to flt4)
data0$tax_all<-read.table(input0$tax, header=T, row.names=1, check.names=FALSE, sep="\t")
dim(data0$tax_all)

data0$taxdf=data0$tax_all
data0$taxdf$Species[is.na(data0$taxdf$Species)]<- "sp."; head(data0$taxdf)
data0$taxdf$Species<-paste(data0$taxdf$Genus, data0$taxdf$Species); head(data0$taxdf)
data0$taxdf$Species<-with(data0$taxdf, ifelse(Species=="NA sp.",paste(data0$taxdf$Family, "sp."), Species )); head(data0$taxdf)
data0$taxdf$Species<-with(data0$taxdf, ifelse(Species=="NA sp.", paste(data0$taxdf$Order, "sp."), Species )); head(data0$taxdf)
data0$taxdf$Species<-with(data0$taxdf, ifelse(Species=="NA sp.", paste(data0$taxdf$Class, "sp."), Species )); head(data0$taxdf)
data0$taxdf$Species<-with(data0$taxdf, ifelse(Species=="NA sp.", paste(data0$taxdf$Phylum, "sp."), Species )); head(data0$taxdf)
data0$taxdf$Species<-with(data0$taxdf, ifelse(Species=="NA sp.", paste(data0$taxdf$Kingdom, "sp."), Species )); head(data0$taxdf)
#View(data0$taxdf)
data0$tax_all=data0$taxdf
dim(data0$tax_all)
```

```{r save work}
save.image("BatchAdjust_clean.RData")
```

##Filtering out samples and TAXA with low presence
```{r filter data}
#filtering
dim(data0$counts_all)
sort(colSums(data0$counts_all))
  #also for recording stats while filtering
  stats=NA
  stats$before=print(c(paste(c("ASVs before filtering:", "samples before filtering:") , 
                           dim(data0$counts)) , paste("minASVabund before filtering:",
                                                      min(rowSums(data0$counts)))))
  stats$before
#filtering samples with reads below minreads (sampling depth):
   stats$depth=print("filtering samples based on number of reads")
minreads=400 #choice of min reads/sample
  stats$minreads=print(paste("removing samples with less than", minreads, 
                             "reads/sample (sample depth)"))
#print(paste("min sampling depth:", min(colSums(data0$counts))))
sampdepth = max(min(colSums(data0$counts_all)) -1, minreads) #smp with at least minreads value or (the min numb of reads in a sample)-1
data0$flt=data0$counts_all[, (colSums(data0$counts_all) >= sampdepth) ]
data0$flt[1:3,(ncol(data0$flt)-3):ncol(data0$flt)]; dim(data0$flt)
a <- which(colSums(data0$flt)<=min(colSums(data0$flt))); names(a) #which sample has lowest # of reads now after flt?
  stats$fltd=print(paste("samples after read depth filtering:", dim(data0$flt)[2], 
                         "(out of", dim(data0$counts)[2], "initial samples)"))
  stats$minsmp=print(paste("after read depth filtering, sample with least reads is now", names(a), "with", 
                           min(colSums(data0$flt)), "reads"))
dim(data0$flt)

#filtering TAX appearing in less than 3 samples (prevalence) with individual abun 0
    stats$asv=print("filtering TAXa based on prevalence in samples")
nsamp=3
minabd=2

   stats$fltasv=print(paste("filtering TAXa thresholds chosen: 1) Prevalence of taxa:", nsamp, 
                          "(taxa appearing in at least this many samples); and 2) minimum abundance of taxa per sample:", minabd))
dim(data0$flt)
prevth=nsamp*100/ncol(data0$flt); prevth #is in %
  stats$prev=print(paste0(nsamp, " is ", round(prevth, 3), "% of the total ",
                          ncol(data0$flt)," samples"))
prevdf = apply(data0$flt, MARGIN=1, FUN = function(x){100*sum(x>minabd)/length(x)})
data0$flt=data0$flt[ prevdf > prevth , ]
dim(data0$flt)
  stats$prev1=print(paste("after taxa filtering, most prevalent ASV is observed in",
                          round(max(prevdf), 3), "% of the samples"))
  stats$prev2=print(paste("after taxa filtering, the least prevalent ASV is observed in", 
                          round(min(prevdf), 3), "% of the samples"))
b<- which(rowSums(data0$flt)<=min(rowSums(data0$flt))); names(b)
c<-  which(rowSums(data0$flt)==max(rowSums(data0$flt))); names(c)

stats$minasv=print(paste("ASV with lowest abunance was", names(b), 
                           ", which was represented by a total of", 
                           min(rowSums(data0$flt)), "reads accross all samples") )
stats$maxasv=print(paste("ASV with highest abunance was", names(c), 
                           ", which was represented by a total of", 
                           max(rowSums(data0$flt)), "reads accross all samples") )
  stats$prev3=print(paste("most prevalent ASV was observed in", 
                          round(max(prevdf), 9), "% of the samples"))
  stats$after1=print(paste("sample with lowest read count after filtering:", 
                         names(a), "with", min(colSums(data0$flt)), "reads"))
  stats$after2=print(paste("# ASVs after filtering:", dim(data0$flt)[[1]]) )

capture.output(stats, file=paste0(output0$flt, "stats_flt-v3.txt"))#, col.names = F, quote=F, row.names = F, sep="\n", eol="\n")

#Reads length distribution through samples & hist plot
tmp=as.data.frame(colSums(data0$flt))
tmp$SampleID=row.names(tmp)
names(tmp)=c("Reads", "SampleID")
tmp=tmp[order(tmp$Reads),]
h=hist(tmp$Reads,main="Sample Distribution after filtering") ; plot(h)
#pdf(paste0(output0$curve, "ReadLenthDistrib-AF.pdf")); plot(h); dev.off()

```

```{r subsetting and exporting filtered data}
#filtered counts
write.table(data0$flt, paste0(output0$flt, "ASV_counts_flt-v3.txt"),  
            sep="\t", col.names=NA, quote=F, na="")

#filtered metadata
data0$meta=data0$meta_all[ colnames(data0$flt),  ]; dim(data0$meta)
write.table(data0$meta, paste0(output0$flt, "ASV_metadata_flt-v3.txt"), sep="\t", 
              col.names=NA, quote=F, na="")

#filtered taxonomies/ASVs
data0$tax=data0$tax_all[rownames(data0$flt), ]; dim(data0$tax)
write.table(data0$tax, paste0(output0$flt, "ASV_taxonomy_flt-v3.txt"),  
            sep="\t", col.names=NA, quote=F, na="NA")
```

###Other/useless summary statistics
```{r}
library(dplyr)
stats$counts_all=NA
stats$counts_all$x=melt(data0$counts_all)
stats$counts_all$totcount=sum(stats$counts_all$x$value); stats$counts_all$totcount
stats$counts_all$std=round(sd(stats$counts_all$x$value),2); stats$counts_all$std
stats$counts_all$sumlist=as.list(round(summary(stats$counts_all$x$value), 2)); stats$counts_all$sumlist
stats$counts_all$density=density(stats$counts_all$x$value); stats$counts_all$density=round(max(summary(stats$counts_all$density$y)), 4); stats$counts_all$density
stats$counts_all$summary=as.data.frame(c(
                                    "Total number of reads in dataset:", stats$counts_all$totcount,
                                    "Min number of reads/TAX", stats$counts_all$sumlist$Min.,
                                    "Max number of reads/TAX", stats$counts_all$sumlist$Max.,
                                    "Median", stats$counts_all$sumlist$Median,
                                    "Mean", stats$counts_all$sumlist$Mean,
                                    "standard deviation:", stats$counts_all$std,
                                    "Table density", stats$counts_all$density))
capture.output(print("more stats: ###before filtering:"), file=paste0(output0$flt, "stats_flt-v3.txt"), append = T)
capture.output(stats$counts_all$summary, file=paste0(output0$flt, "stats_flt-v3.txt"), append = T)


stats$flt=NA
stats$flt$x=melt(data0$flt)
stats$flt$totcount=sum(stats$flt$x$value); stats$flt$totcount
stats$flt$std=round(sd(stats$flt$x$value),2); stats$flt$std
stats$flt$sumlist=as.list(round(summary(stats$flt$x$value), 2)); stats$flt$sumlist
stats$flt$density=density(stats$flt$x$value); stats$flt$density=round(max(summary(stats$flt$density$y)), 4); stats$flt$density
stats$flt$summary=as.data.frame(c(
                                    "Total number of reads in dataset:", stats$flt$totcount,
                                    "Min number of reads/TAX", stats$flt$sumlist$Min.,
                                    "Max number of reads/TAX", stats$flt$sumlist$Max.,
                                    "Median", stats$flt$sumlist$Median,
                                    "Mean", stats$flt$sumlist$Mean,
                                    "standard deviation:", stats$flt$std,
                                    "Table density", stats$flt$density))
capture.output(print("more stats: ###after filtering:"), file=paste0(output0$flt, "stats_flt-v3.txt"), append = T)
capture.output(stats$flt$summary, file=paste0(output0$flt, "stats_flt-v3.txt"), append = T)

```

### rarecurves exploration before and after filtering
```{r rarecurves before filtering}
library(ampvis2)
output0$curve=paste0(output0$path, "rarecurves/")
dir.create(output0$curve)

#================ rarecureves before filtering
    #PREPPING the ampvis objsect
    dim(data0$counts_all); dim(data0$tax_all)
    data0$coutax=cbind(data0$counts_all, data0$tax_all)
    dim(data0$coutax) #+ 7 columns of ranks
    data0$meta_all$SampleID=row.names(data0$meta_all)
    relocate(data0$meta_all, SampleID) -> data0$meta_all; colnames(data0$meta_all)
    dim(data0$meta_all)
data0$amp=amp_load(data0$coutax, data0$meta_all)


#rarefication curve before filtering
rarecurve <- amp_rarecurve(data0$amp, color_by = "StudyID",facet_by = "StudyID") + xlim(0,25000) +
  ylab("Number of ASVs") + xlab("Number of reads (sequencing depth)"); rarecurve
pdf(paste0(output0$curve, "rcurve_BFflt_xlim=25K_StudyIDampvis-v3.pdf"), height=10, width=14); rarecurve; dev.off()
#\------------- done with "before"


#=============== rarecureves after filtering
    #PREPPING the ampvis objsect
    dim(data0$flt); dim(data0$tax)
    data0$coutax=cbind(data0$flt, data0$tax)
    dim(data0$coutax) #+ 7 columns of ranks
    data0$meta$SampleID=row.names(data0$meta)
    relocate(data0$meta, SampleID) -> data0$meta; colnames(data0$meta)
    dim(data0$meta)
data0$amp=amp_load(data0$coutax, data0$meta)

#rarefication curve #2 (one with ampvis2)
rarecurve <- amp_rarecurve(data0$amp, color_by = "StudyID",facet_by = "StudyID") + xlim(0,25000) +
  ylab("Number of ASVs") + xlab("Number of reads (sequencing depth)"); rarecurve
pdf(paste0(output0$curve, "rcurve_AFflt_xlim=25K_StudyIDampvis-v3.pdf"), height=10, width=14); rarecurve; dev.off()
#\----------------- done with "after"
```


############### BATCH EFFECTS EVALUATION ######################
Exploring the strength of the different batches before batch adjustment
```{r}
output0$beval=paste0(output0$path,"BatchEvaluation/")
dir.create(output0$beval)

batch0=NA
batch0$all=c("StudyID","InfectionStatus", "InfectionType", "Gender", "AgeCategory", "HelminthSpecies", "HelminthCohort")

#evaluating batch effect on flt counts
D_before <- vegdist(t(data0$flt))
    #dip toe
fit_adonis <- adonis(D_before ~ data0$meta$StudyID); fit_adonis
    #full dip
for (b in batch0$all){
set.seed(1)
fit_adonis <- adonis(D_before ~ data0$meta[[b]]); fit_adonis
capture.output(fit_adonis, file=paste0(output0$beval, "gPMVA-", b , "~flt_", "_BatchAdj=NONE.txt"), quotes=F)
}
 #where differences in the factor b can expalin a total of R2*100% of the variability in microbial abundance profiles => factor b contributes R2*100% of the differences in microbial profiles
#in this case, highest contributor was STUDY with close followup of HelminthCohort (at 15% and HelminthSpecies at 13%)

#now with anosim
for (b in batch0$all){
set.seed(1)
fit_anosim <- anosim(D_before, data0$meta[[b]]); fit_anosim
capture.output(fit_anosim, file=paste0(output0$beval,"gANSM",  b , "~flt_", "_BatchAdj=NONE.txt"), quotes=F)
}
```
Here we learned that Batches of StudyID is identical to Country (makes sense) => removing "Country" as factor

############### BATCH EFFECTS adjustments ######################
Adjusting batch effect
```{r}
library(MMUPHin)
data0$meta$StudyID=factor(data0$meta$StudyID) #MMUphin needs this to be a factor
fit_adjust_batch <- adjust_batch(feature_abd = as.matrix(data0$flt),
                                 batch = "StudyID",
                                 covariates = c("InfectionType", "AgeCategory"), #explore to reduce R2 as much as possible
                                 data = data0$meta,
                                 control = list(verbose = F))
data0$adj=fit_adjust_batch$feature_abd_adj
```

############### BATCH EFFECTS EVALUATION ######################
#after batch adjustment
```{r}
output0$batch=paste0(output0$path,"BatchADJusted/")
dir.create(output0$batch)

write.table(data0$adj, paste0(output0$batch, "ASV_counts_ADJ-v3.txt"),  
            sep="\t", col.names=NA, quote=F, na="")

##after adjustment
E_after <- vegdist(t(data0$adj))#, method="euclidean")
 #dip toe
fit_adonis <- adonis(E_after ~ data0$meta$StudyID); fit_adonis; hist(fit_adonis$f.perms)
fit_adonis <- adonis(E_after ~ data0$meta$InfectionType); fit_adonis; hist(fit_adonis$f.perms)
fit_adonis <- adonis(E_after ~ data0$meta$AgeCategory); fit_adonis; hist(fit_adonis$f.perms)
 
 #full dip
#with PERMANOVA R2 (adonis)
for (b in batch0$all){
set.seed(1)
fit_adonis <- adonis(E_after ~ data0$meta[[b]]); fit_adonis
capture.output(fit_adonis, file=paste0(output0$batch, "gPMVA-", b , "~adj_", "_BatchAdj=mmuphin.txt"), quotes=F)
}

#or with anosim
for (b in batch0$all){
set.seed(1)
fit_anosim <- anosim(E_after, data0$meta[[b]]); fit_anosim
capture.output(fit_anosim, file=paste0(output0$batch,"gANSM",  b , "~adj_", "_BatchAdj=mmuphin.txt"), quotes=F)
}
```


#prep colors
```{r colors}
#Create a distinctive color pallete (a color_vector)
library(RColorBrewer)
  qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',] 
  col_vector = unique(unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals))))
  n=25 #pick between 2 and 70
  pie(rep(1,n), col=col_vector, labels = col_vector, cex=0.8)
  mycol=c("orangered","green", "yellow4", "palegreen4", "orange","blue", "paleturquoise3", "peru", "yellow", 
          "purple", "royalblue1", "salmon",  "seagreen1" , "seashell3" , "salmon4","pink", "skyblue2", "gray40")
  pie(rep(1,length(mycol)), col=mycol, labels = mycol, cex=0.8);
  col=c(mycol, col_vector)
  pie(rep(1,n), col=col, labels = col, cex=0.8);
  #col=replace(col_vector, col_vector %in% c("#FFFF99", "#D95F02", "#A6CEE3"), c("yellow","orange", "peachpuff" ))
  #pie(rep(1,n), col=col, labels = col, cex=0.8);
  #tol=sample(col_vector, n);   tol
#Create pch vector that I like (up to 19 values):
  pich=c(2:20, 1)
```


```{r PCoA for batch effect before and after adjustment}
data0$meta$group1=NULL; data0$meta$group2=NULL #cleaning custom colums, if such exist

#================ PcOA before batch adjustments (after filtering)
     #PREPPING the ampvis objsect
    dim(data0$flt); dim(data0$tax)
    data0$flttax=cbind(data0$flt, data0$tax)
    dim(data0$flttax) #+ 7 columns of ranks
    data0$meta$SampleID=row.names(data0$meta)
    relocate(data0$meta, SampleID) -> data0$meta; colnames(data0$meta)
    dim(data0$meta)
data0$amp=amp_load(data0$flttax, data0$meta)

#PCOA plots
set.seed(123)
pcoa_bray <- amp_ordinate(data0$amp, filter_species = 0.01, type = "PCOA",
    sample_color_by = "StudyID", sample_shape_by = "InfectionStatus", 
    detailed_output = TRUE, transform = "none", distmeasure = "bray")
pcoa_bray$plot= pcoa_bray$plot + scale_color_manual(values=col) + scale_shape_manual(values=pich) ; pcoa_bray$plot
pdf(paste0(output0$batch, "PCoA_StudyID~InfStatus_AFflt_beforeADJ",".pdf"), width=10, height=8); pcoa_bray$plot; dev.off()
# pcoa_bray$figcaption

#NMDS
nmds_bray <- amp_ordinate(data0$amp, filter_species = 0.01, type = "NMDS",
     sample_color_by = "StudyID",sample_shape_by = "InfectionStatus", 
    detailed_output = TRUE, transform = "none", distmeasure = "bray")
nmds_bray$plot=nmds_bray$plot + scale_color_manual(values=col) + scale_shape_manual(values=pich); nmds_bray$plot
pdf(paste0(output0$batch, "NMDS_StudyID~InfStatus_AFflt_beforeADJ", ".pdf")); nmds_bray$plot; dev.off()
#\------------------- done with "before" ordination


#subset if needed
    dim(data0$adj); dim(data0$tax)
    data0$adjtax=cbind(data0$adj, data0$tax)
    dim(data0$adjtax) #+ 7 columns of ranks
    data0$meta$SampleID=row.names(data0$meta)
    relocate(data0$meta, SampleID) -> data0$meta; colnames(data0$meta)
    dim(data0$meta)
data0$amp=amp_load(data0$adjtax, data0$meta)


#PCOA plots
set.seed(123)
pcoa_bray <- amp_ordinate(data0$amp, filter_species = 0.01, type = "PCOA",
    sample_color_by = "StudyID", sample_shape_by = "InfectionStatus", 
    detailed_output = TRUE, transform = "none", distmeasure = "bray", scale_x_reverse())
pcoa_bray$plot= pcoa_bray$plot + scale_color_manual(values=col) + scale_shape_manual(values=pich) ; pcoa_bray$plot
pcoa_bray$plot= pcoa_bray$plot + scale_x_reverse()+ scale_y_reverse() ; pcoa_bray$plot
pdf(paste0(output0$batch, "PCoA_StudyID~InfStatus_AFflt_afterADJ",".pdf"), width=10, height=8); pcoa_bray$plot; dev.off()
# pcoa_bray$figcaption

#NMDS
nmds_bray <- amp_ordinate(data0$amp, filter_species = 0.01, type = "NMDS",
     sample_color_by = "StudyID",sample_shape_by = "InfectionStatus", 
    detailed_output = TRUE, transform = "none", distmeasure = "bray")
nmds_bray$plot=nmds_bray$plot + scale_color_manual(values=col) + scale_shape_manual(values=pich); nmds_bray$plot
pdf(paste0(output0$batch, "NMDS_StudyID~InfStatus_AFflt_afterADJ", ".pdf")); nmds_bray$plot; dev.off()
```


#DifAbund with MMUphin
```{r diff abundance with mmuphin}
library(MMUPHin); library(Maaslin2); library(vegan)
#lm_meta() only works on variates with 2 categories, so that disqualifies it from a lot of analyses


data0$meta$StudyID=factor(data0$meta$StudyID)
data0$meta$InfectionType=factor(data0$meta$InfectionType, levels = c("None", "Single-species", "Multi-species"))
data0$meta$Gender=factor(data0$meta$Gender) #the factor of interest has to have 2 levels and be factored
#=========== uninfected 
data0$subset$meta=subset(data0$meta, InfectionStatus=="Uninfected", drop=T);
      data0$subset$meta[1:4, 3:8]; dim(data0$subset$meta)
data0$subset$flt=data0$flt[, rownames(data0$subset$meta) ]; 
#lm_meta() only works on variates with 2 categories
fit_lm_meta <- lm_meta(feature_abd = as.matrix(data0$subset$flt),
                       batch = "StudyID", exposure = "Gender",
                       covariates = "AgeCategory", #Age.Category has only 1 category for Yang_2017
                       data = data0$subset$meta,
                       control = list(verbose = T))
meta_fits <- fit_lm_meta$meta_fits
dim(meta_fits); meta_fits[1:6, ]
#replacing TAXid with lineage info
data0$tax[1:6, 1:7]
meta_fits$tax=data0$tax$Species
dim(meta_fits); meta_fits[1:6,]

#this will create a "lm_meta" folder in the R-path directory. Use if needed

da=meta_fits %>% 
  filter(qval.fdr < 0.05) %>% 
  arrange(coef) %>% 
  mutate(feature = factor(feature, levels = feature)) %>% 
  ggplot(aes(y = coef, x = tax)) +
  geom_bar(stat = "identity") + labs(title="DiffAbund of Taxa in various infected Samples ~ Uninfected") + 
  coord_flip() + theme_bw(); plot(da) 
#empty plot cuz no values
#pdf(paste0(output0$batch, "DiffAbubd_adj-variousInfected~Uninfected.pdf"), height=8, width=12); plot(da); dev.off()


#============================= infected
data0$subset$meta=subset(data0$meta, InfectionStatus=="Infected", drop=T);
      data0$subset$meta[1:4, 3:8]; dim(data0$subset$meta)
data0$subset$meta$InfectionType=factor(data0$subset$meta$InfectionType)
data0$subset$flt=data0$flt[, rownames(data0$subset$meta) ]; 
#lm_meta() only works on variates with 2 categories
fit_lm_meta <- lm_meta(feature_abd = as.matrix(data0$subset$flt),
                       batch = "StudyID", exposure = "InfectionType",
                       #covariates = "AgeCategory",
                       data = data0$subset$meta,
                       control = list(verbose = T))
meta_fits <- fit_lm_meta$meta_fits
dim(meta_fits); meta_fits[1:6, ]
#replacing TAXid with lineage info
data0$tax[1:6, 1:7]
meta_fits$tax=data0$tax$Species
dim(meta_fits); meta_fits[1:6,]

#this will create a "lm_meta" folder in the R-path directory. Use if needed

da=meta_fits %>% 
  filter(qval.fdr < 0.05) %>% 
  arrange(coef) %>% 
  mutate(feature = factor(feature, levels = feature)) %>% 
  ggplot(aes(y = coef, x = tax)) +
  geom_bar(stat = "identity") + labs(title="DiffAbund of Taxa in various infected Samples ~ Uninfected") + 
  coord_flip() + theme_bw(); plot(da)
#getting NA results
#pdf(paste0(output0$batch, "DiffAbubd_adj-variousInfected~Uninfected.pdf"), height=8, width=12); plot(da); dev.off()


```

#Identify descrete structures:
```{r}
output0$dstr=paste0(output0$batch,"DiscreteStr/")
dir.create(output0$dstr)
# all studies, all variables

##=========== before batch adjustment
E_subset <- vegdist(t(data0$flt))
fit_discrete <- discrete_discover(D = E_subset,
                                  batch = "StudyID", data = data0$meta,
                                  control = list(k_max = 8, verbose = FALSE))
#this will create a "discrete_diagnostic.pdf" file in R-path directory. Compare internal to external. 
#"A cluster number with good support from both internal and external evaluations provides
#meta-analytical evidence for discrete structures in the microbial abundance profiles."

#note the "set K to #" line, or do
fit_discrete$internal_mean #count the clusters
fit_discrete$external_mean

#visualization
k=max(row.names(fit_discrete$internal_mean))
 st=colnames(fit_discrete$internal_mean); st
#st="Easton_2019"
for (i in st){
internal <- data.frame(
  # By default, fit_discrete evaluates cluster numbers 2-10
  K = 2:k,
  statistic = 
    fit_discrete$internal_mean[, i],
  se = 
    fit_discrete$internal_se[, i],
  type = "internal")
external <- data.frame(
  # By default, fit_discrete evaluates cluster numbers 2-10
  K = 2:k,
  statistic = 
    fit_discrete$external_mean[, i],
  se = 
    fit_discrete$external_se[, i],
  type = "external")
p=rbind(internal, external) %>% 
  ggplot(aes(x = K, y = statistic, color = type)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  geom_line(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = statistic - se, ymax = statistic + se),
                position = position_dodge(width = 0.5), width = 0.5) +
  ggtitle(paste("Evaluation of discrete structure in microbiomes for", i)); p
pdf(paste0(output0$dstr, "DctStrs_flt~", i, "_BFadj.pdf"), height=8, width=12); 
plot(p); dev.off()
}
#\------------------------------
 
 
#=================== after batch adjustment
 E_subset <- vegdist(t(data0$adj))
fit_discrete <- discrete_discover(D = E_subset,
                                  batch = "StudyID", data = data0$meta,
                                  control = list(k_max = 8, verbose = FALSE))
#this will create a "discrete_diagnostic.pdf" file in R-path directory. Compare internal to external. 
#"A cluster number with good support from both internal and external evaluations provides
#meta-analytical evidence for discrete structures in the microbial abundance profiles."

#note the "set K to #" line, or do
fit_discrete$internal_mean #count the clusters
fit_discrete$external_mean

#visualization
k=max(row.names(fit_discrete$internal_mean))
 st=colnames(fit_discrete$internal_mean); st
#st="Easton_2019"
for (i in st){
internal <- data.frame(
  # By default, fit_discrete evaluates cluster numbers 2-10
  K = 2:k,
  statistic = 
    fit_discrete$internal_mean[, i],
  se = 
    fit_discrete$internal_se[, i],
  type = "internal")
external <- data.frame(
  # By default, fit_discrete evaluates cluster numbers 2-10
  K = 2:k,
  statistic = 
    fit_discrete$external_mean[, i],
  se = 
    fit_discrete$external_se[, i],
  type = "external")
p=rbind(internal, external) %>% 
  ggplot(aes(x = K, y = statistic, color = type)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  geom_line(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = statistic - se, ymax = statistic + se),
                position = position_dodge(width = 0.5), width = 0.5) +
  ggtitle(paste("Evaluation of discrete structure in microbiomes for", i)); p
pdf(paste0(output0$dstr, "DctStrs_adj~", i, "_AFadj.pdf"), height=8, width=12); 
plot(p); dev.off()
}
```
Evidence for existence of discrete structures would be a “peaking” of the mean statistics at a particular cluster number. The decreasing trend for both the internal and external statistics along with number of clusters (K) suggests that discrete structures cannot be well-established.

discrete_discover takes as input sample-by-sample dissimilarity measurements (generated from microbial abundance profiles), and performs unsupervised clustering within each batch across a range of cluster numbers. ***It then evaluates the support for each cluster number with both internal (i.e., samples within the batch) and external (i.e., samples in other batches) data. Internal evaluation is realized with prediction.strength and external evaluation is based on a generalized version of the same method.*** discrete_discover generates as output the evaluation statistics for each cluster number. A cluster number with good support from both internal and external evaluations provides meta-analytical evidence for discrete structures in the microbial abundance profiles.


##continuious structures discovery
# Much like adjust_batch and lm_meta, continuous_discover also takes
# as input feature-by-sample abundances. control offers many tuning parameters
# and here we set one of them, var_perc_cutoff, to 0.5, which asks the method
# to include top principal components within each batch that in total explain
# at least 50% of the total variability in the batch. See 
# help(continuosu_discover) for more details on the tuning parameters and 
# their interpretations.
```{r}
output0$cstr=paste0(output0$batch,"ContStr/")
dir.create(output0$cstr)
# First subset both feature abundance table and metadata to only control samples
subset="Ininfected"
data0$subset$meta <- subset(data0$meta, InfectionStatus == subset)
data0$subset$adj <- data0$adj[, rownames(data0$subset$meta)]
dim(data0$subset$meta); dim(data0$meta)

#install.packages("igraph")
library(igraph); packageVersion("igraph")
#install.packages("Rglpk", pkgs = "~/Documents/MyApps/GLPK-5.0/Rglpk_0.6-4.tar.gz" )
library(Rglpk); packageVersion("Rglpk")
library(phyloseq); library(MMUPHin); library(Maaslin2); library(vegan)
#real data
        data0$subset$meta$StudyID=factor(data0$subset$meta$StudyID)
fit_continuous <- continuous_discover(feature_abd = data0$subset$adj,
                                      batch = "StudyID",
                                      data = data0$subset$meta,
                                      control = list(var_perc_cutoff = 0.5,
                                                     verbose = T))
fcld=fit_continuous$consensus_loadings; fcld
row.names(fcld)=data0$tax$Species

which_loading=1
loading <- data.frame(feature = rownames(fcld),
                      loading = fcld[, which_loading])
(head(loading))
#the features / TAXA gradient
ld=loading %>% dplyr::arrange(-abs(loading)) %>%
    dplyr::slice(1:20) %>%  dplyr::arrange(loading) %>%
    dplyr::mutate(feature = factor(feature, levels = feature)) %>%
    ggplot(aes(x = feature, y = loading)) +
    geom_bar(stat = "identity") + theme_bw() +
    coord_flip() + ggtitle("Top 20 taxonomic features that produce the observed gradient"); plot(ld)
pdf(paste0(output0$batch, paste0("ContStruTaxaLoading",which_loading, "-", subset, ".pdf")), height=8, width=12); plot(ld); dev.off()

#The PCoA plot
E_subset <- vegdist(t(data0$subset$adj))
mds <- cmdscale(d = E_subset)
colnames(mds) <- c("Axis1", "Axis2")
fcsc=fit_continuous$consensus_scores[, which_loading]

pcoa=as.data.frame(mds) %>% 
  dplyr::mutate(score1 = fcsc) %>% 
  ggplot(aes(x = Axis2, y = Axis1, color = score1, shape=as.factor(data0$subset$meta$StudyID))) +
  geom_point() + theme_bw() + scale_shape_manual(values=pich, name="StudyID")  +
  coord_flip(); plot(pcoa)# 
pdf(paste0(output0$batch, "ContStruPCoA-", subset ,".pdf"), height=8, width=10); plot(pcoa); dev.off()

```
#test
        data("CRC_abd", "CRC_meta")
        fit_adjust_batch <- adjust_batch(feature_abd = CRC_abd, batch = "studyID", covariates = "study_condition", 
                                         data = CRC_meta, control = list(verbose = FALSE))
        CRC_abd_adj <- fit_adjust_batch0$feature_abd_adj
        control_meta <- subset(CRC_meta, study_condition == "control")
        control_abd_adj <- CRC_abd_adj[, rownames(control_meta)]
        fit_continuous <- continuous_discover(feature_abd = control_abd_adj, batch = "studyID",    
                            data = control_meta, control = list(var_perc_cutoff = 0.5, verbose = FALSE))
```{r save work}
save.image("BatchAdjust_clean.RData")
```