
library(qiimer)
library(vegan)
library(permute)
library(lattice)
library(tidyverse)
library(qiime2R)
library(phyloseq)
library(microbiomeMarker)
library(RColorBrewer)
library(data.table)
library(ggpubr)
library(cowplot)
library(MicEco)
library(FEAST)
library(Rcpp)
library(RcppArmadillo)
library(ggthemes)
library(gridExtra)
library(reshape2)
library(conflicted)

####BETA DIVERSITY ANALYSES####
####load distance matrices and metadata####

uw_dm<-read_qiime_distmat('unw-distance-matrix.tsv')##unweighted unifrac distance matrix
w_dm<-read_qiime_distmat('we-distance-matrix.tsv') ##weighted unifrac distance matrix
map <- read_qiime_mapping_file("metadata_bromeliads.tsv")##metadata

####code for running permanovas####
perm <- how(nperm = 999)
setBlocks(perm) <- with(map, Individual)

adonis2(uw_dm ~ Species, data=map, 
        method= "unweighted", by = "margin", permutations=999)#for unweighted unifrac distance matrix

adonis2(w_dm ~ Species, data=map, 
        method= "weighted", by = "margin", permutations=999)#for weighted unifrac distance matrix

####within-community dispersion analyses with betadisper####
bet_uw=betadisper(uw_dm, map$Species)
permutest(bet_uw, pairwise=FALSE, permutations=1000)#unweighted unifrac


bet_we=betadisper(w_dm,map$Species)
permutest(bet_we, pairwise=FALSE, permutations=1000)#weighted uniFrac


####PHYLOSEQ: TAXONOMIC ANALYSES AND GRAPHS####

phyb<-qza_to_phyloseq("merged-table-final.qza", 
                      "bromeliads-rooted-tree.qza", 
                      "bromeliads-taxonomy.qza",
                      "metadata_bromeliads.tsv")#import working files to create 
                                                #a phyloseq object from qiime2 files (.qza)

phyb ##phyloseq object

####differential abundance test with lfse####
#lfse test at phylum level

m1<-run_lefse(phyb, group = "Species", taxa_rank = "Phylum",
              transform = "log10", norm= "CPM",
              kw_cutoff = 0.05, lda_cutoff = 2)

m1

df.lfse<-marker_table(m1) %>% data.frame()
head(df.lfse)
print(df.lfse)

#lfse test at family level
m2<-run_lefse(phyb, group = "Species", taxa_rank = "Family",
              transform = "log10", norm= "CPM",
              kw_cutoff = 0.05, lda_cutoff = 2)
m2
df.lfse.fam<-marker_table(m2) %>% data.frame()
head(df.lfse.fam)
print(df.lfse.fam)

####taxonomy plots####
##code for creating phylum stackbar
pstrans <- transform_sample_counts(phyb, function(otu) otu/sum(otu))#transform reads to relative abundance

glom_taxa <- tax_glom(pstrans, taxrank = "Phylum")
dat_taxa <- data.table(psmelt(glom_taxa))
dat_taxa$Phylum <- as.character(dat_taxa$Phylum)
medians <- dat_taxa[, median := median(Abundance, na.rm = TRUE), by = "Phylum"]
remainder <- dat_taxa[(median <= 0.5e-25), Phylum := "Other"]
colourCount = length(unique(dat_taxa$Phylum))

phyl.plot <- ggplot(dat_taxa, aes(x = Sample, y=Abundance, fill=Phylum, order=as.factor(Phylum)))
phyl.plot <- phyl.plot + geom_bar(aes(), stat="identity", position="fill") + 
    facet_wrap(~ factor(sample_Species, levels = c("Bromeliad", "Austrophorocera", "Copestylum")),  scales = "free_x") + 
    theme_bw() + theme(legend.position="bottom", axis.text.y = element_text(size=10, colour = "black"),
                       axis.title.x = element_text(size=17), axis.title.y = element_text(size=17), 
                       axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, color = "black"), 
                       panel.grid.minor = element_blank(), strip.text.x = element_text(size = 17),
                       legend.text=element_text(size=14),
                       legend.title=element_text(size=14)) + 
    scale_fill_brewer(palette = "Paired", labels=c('Acidobacteriota', 'Actinomycetota',
                                                   'Bacteroidota', 'Epsilonbacteraeota',
                                                   'Bacillota', 'Other', 'Patescibacteria',
                                                   'Pseudomonadota','Spirochaetota',
                                                   'Verrucomicrobiota')) + labs(y = "Proportion") 
phyl.plot

##code for creating a family plot
glom_taxa2 <- tax_glom(pstrans, taxrank = "Family")
dat_taxa2 <- data.table(psmelt(glom_taxa2))
dat_taxa2$Family <- as.character(dat_taxa2$Family)
medians2 <- dat_taxa2[, median := median(Abundance, na.rm = TRUE), by = "Family"]
remainder2 <- dat_taxa2[(median <= 0.5e-60), Family := "Other"]
colourCount = length(unique(dat_taxa2$Family))

br.factor2 <- factor(dat_taxa2$Sample, c("CR1", "CR2", "CR3","CR4", "CR5", "CR6", "CR7", "CR8", "CR9","CR10","CR11", "CR12", 
                                         "CR13", "CR14", "CR15", "CR16", "CR17", "CR18", "CR19", "CR20", "CR21", "CR22",
                                         "CR23", "CR24", "CR25", "CR27", "CR28", "CR29", "CR31", "CR33", "CR34", "CR35", 
                                         "CR36", "CR37", "CR38", "CR41", "CR42"))##This is to reorder samples

color.axis<-c("darkorchid", "darkorchid", "darkorchid","darkorchid","darkorchid",
              "darkorchid","darkorchid","darkorchid","darkorchid","darkorchid",
              "darkorchid","darkgoldenrod4","darkgoldenrod4","darkgoldenrod4",
              "darkgoldenrod4","darkgoldenrod4", "darkgoldenrod4","darkgoldenrod4","darkgoldenrod4", "darkgoldenrod4", 
              "darkgoldenrod4","darkgreen","darkgreen","darkgreen","darkgreen","darkgreen","darkgreen",
              "darkgreen","darkgreen","darkgreen","darkgreen","darkgreen","darkgreen",
              "darkgreen","darkgreen","darkgreen")

fam.plot <- ggplot(dat_taxa2, aes(br.factor2, y=Family, fill=Abundance))
fam.plot <- fam.plot + geom_tile() + theme(axis.text.y = element_text(size=9, colour = "black", face= "italic"),
                                           axis.title.x = element_blank(), 
                                           axis.title.y = element_text(size=17), 
                                           axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=7, color = color.axis),
                                           legend.text=element_text(size=12),
                                           legend.title=element_text(size=12)) +
            scale_fill_gradient(low = "white", high = "black") +
            labs(y= "Microbial families", fill = "Proportion")


fam.plot


###CODES FOR CREATING DIVERSITY PLOTS####
####beta diversity (pcoa) plots####
pstrans

ps.un <- ordinate(pstrans, "PCoA", "unifrac", weighted = FALSE)
ps.we <- ordinate(pstrans, "PCoA", "unifrac", weighted = TRUE)

un.fig<-plot_ordination(pstrans, ps.un, color="Species") +
    geom_point(size = 4, shape= 1, alpha=0.9) + theme_bw() +
    scale_color_manual(values = c("darkorchid2","chartreuse4","darkgoldenrod3")) + 
    theme(axis.text = element_text(size=14, colour = "black"), 
          axis.title.x = element_text(size=14, colour = "black"),
          axis.title.y = element_text(size=14, colour = "black"),
          legend.position = "none",
          legend.text = element_text(size=14),
          plot.title = element_text(size = 11, face = "bold")) + 
    ggtitle("Unweighted UniFrac") + labs(x="PC1 25.9%", y="PC2 18.4%")

un.fig

we.fig<-plot_ordination(pstrans, ps.we, color="Species") +
    geom_point(size = 4, shape= 1, alpha=0.9) + theme_bw() +
    scale_color_manual(values = c("darkorchid2","chartreuse4","darkgoldenrod3")) + 
    theme(axis.text = element_text(size=14, colour = "black"), 
          axis.title.x = element_text(size=14, colour = "black"),
          axis.title.y = element_text(size=14, colour = "black"),
          legend.position = "none",
          legend.text = element_text(size=14),
          plot.title = element_text(size = 11, face = "bold")) + 
    ggtitle("Weighted UniFrac") + labs(x="PC1 21.8%", y="PC2 19.3%")

we.fig


####alpha diversity plots####
df.brom<-read.csv("bromeliads_df.csv", header=T)#import working file
names(df.brom)

level_order <- c('Bromeliad',
                 'Austrophorocera',
                 'Copestylum') 
#observed ASV plot
pb<-ggplot(df.brom, aes(x = Species, y= ASVs)) +  geom_boxplot(color = c("chartreuse4", 
                                                                         "darkorchid2", 
                                                                         "darkgoldenrod3"), lwd= 0.5)

asvplot <-pb + scale_x_discrete(limits = level_order) + theme_bw() +
                labs(y = "Observed ASVs") + theme(axis.text.x = element_text(size=11, colour = "black"),
                                      axis.text.y = element_text(size=11, colour = "black"),
                                      axis.title.x = element_blank(),
                                      axis.title.y = element_text(size=14),
                                      panel.grid.major = element_blank(), 
                                      panel.grid.minor = element_blank(),
                                      panel.border = element_blank()) +
                geom_jitter(width=0.25, alpha=0.5) +
                theme(axis.line = element_line(color = 'black'))

asvplot

#Faith's phylogenetic diversity plot
ppd<-ggplot(df.brom, aes(x = Species, y= PD)) +  geom_boxplot(color = c("chartreuse4", 
                                                                        "darkorchid2", 
                                                                        "darkgoldenrod3"), lwd= 0.5)

pd <-ppd + scale_x_discrete(limits = level_order) + theme_bw() +
            labs(y = "Phylogenetic diversity") + theme(axis.text.x = element_text(size=11, colour = "black"),
                                               axis.text.y = element_text(size=11, colour = "black"),
                                               axis.title.x = element_blank(),
                                               axis.title.y = element_text(size=14),
                                               panel.grid.major = element_blank(), 
                                               panel.grid.minor = element_blank(),
                                               panel.border = element_blank()) +
            geom_jitter(width=0.25, alpha=0.5) +
            theme(axis.line = element_line(color = 'black'))
pd

####panel of alpha and beta diversity plots####
ggarrange(asvplot, pd, un.fig, we.fig,
          labels = c("A", "B", "C","D"),
          ncol = 2, nrow = 2, align= "v")

####FEAST ANALYSIS AND VENN DIAGRAM####

####code for venn diagram####
sample_data(phyb)
venn.fig<-ps_venn(phyb, "Species", fraction = 0.5, weight = F,
                  relative = TRUE, plot = TRUE,
                  labels = list(cex = 1.25, col = c("black","white", "black")),
                  fill = c("white","chartreuse4","white"),
                  edges = list(col = c("darkorchid", "darkgreen","darkgoldenrod4"), lex = 2),
                  quantities = list(font = 1, cex = 1.34))

venn.fig

####code for running feast (fast expectation-maximization microbial source tracking)####
metadata2 <- Load_metadata(metadata_path = "1_metadata_for_FEAST.txt")#import working metadata
otus <-Load_CountMatrix(CountMatrix_path = "2_feature_table_bromeliads_for_FEAST.txt")#import working asv's table

FEAST_output <- FEAST(C = otus, metadata = metadata2, different_sources_flag = 0, 
                      dir_path = "~/Desktop",
                      outfile="demo")#FEAST command. "FEAST will then save the file demo_FEAST.txt - 
#A file containing an S1 by S2 matrix P, where S1 is the number sinks 
#and S2 is the number of sources (including an unknown source)". 
#Taken from https://github.com/cozygene/FEAST/tree/master

#code for importing data frame to create the graph in ggplot2

brom.df <-read.csv("3_feast_graph.csv", header=T)##import working file
names(brom.df)
br.factor <- factor(brom.df$Host, c("Austrophorocera.ID12", "Austrophorocera.ID8",
                                    "Austrophorocera.ID7","Austrophorocera.ID1",
                                    "Austrophorocera.ID10","Austrophorocera.ID5",
                                    "Austrophorocera.ID4","Austrophorocera.ID11",
                                    "Austrophorocera.ID6","Austrophorocera.ID9",
                                    "Austrophorocera.ID2","Austrophorocera.ID3","Copestylum.ID21",
                                    "Copestylum.ID18","Copestylum.ID17","Copestylum.ID20","Copestylum.ID19",
                                    "Copestylum.ID22","Copestylum.ID15","Copestylum.ID16","Copestylum.ID13",
                                    "Copestylum.ID14"))#To reorder samples 


##vector of colorss for the x axis##
colores<- c("darkorchid", "darkorchid", "darkorchid","darkorchid","darkorchid",
            "darkorchid","darkorchid","darkorchid","darkorchid","darkorchid","darkorchid", 
            "darkorchid","darkgoldenrod4","darkgoldenrod4","darkgoldenrod4",
            "darkgoldenrod4","darkgoldenrod4", "darkgoldenrod4","darkgoldenrod4",
            "darkgoldenrod4", "darkgoldenrod4", "darkgoldenrod4")


p1 <- ggplot(brom.df, aes(x = br.factor, y = Proportion, fill = Source)) + 
    geom_bar(stat = "identity")

plotbrom<-p1 + theme_bw() + theme(axis.text.y = element_text(size=10, colour = "black"),
                                  axis.title.x = element_text(size=17),
                                  axis.title.y = element_text(size=17),
                                  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=10, colour = colores),
                                  axis.text.y.left = element_text(size=14, colour = "black"),
                                  legend.text=element_text(size=16),
                                  legend.title=element_text(size=17),
                                  panel.grid.major = element_blank(), 
                                  panel.grid.minor = element_blank(),
                                  panel.border = element_blank()) + 
    theme(axis.line = element_line(color = 'black')) +
    scale_fill_manual(values = c("chartreuse4", "gray88")) + labs(x="")

plotbrom

##Creating a figure with Venn diagram and Feast plot
ggarrange(plotbrom, venn.fig,
          labels = c("A", "B"),
          ncol = 2, nrow = 1, align= "v")



####SUPPLEMENTAL####
#distance to centroid plots
#it uses the objects created in within-community dispersion analyses

dist.uw<-bet_uw[["distances"]]##code to obtain distances to centroid using unweighted unifrac distance matrix
dist.we<-bet_we[["distances"]]##code to obtain distances to centroid using weighted distance matrix

dist.uw<-as.data.frame(dist.uw)#create a data frame with the unweighted distance matrix
dist.we<-as.data.frame(dist.we)#create a data frame with the weighted distance matrix

write.csv(dist.uw, "centroid_uw.csv")#data frame for unweighted unifrac distances to centroid to work in ggplot2
write.csv(dist.we, "centroid_we.csv")#data frame for weighted unifrac distances to centroid to work in ggplot2

#load file to make boxplots in ggplot2#

centbrom <-read.csv("centroid_bromeliads.csv",header = T)#import working file
names(centbrom)
level_order <- c('Bromeliad',
                 'Austrophorocera',
                 'Copestylum') 
#Unweighted plot
cenplot.uw<-ggplot(centbrom, aes(x = Species, y= dist.uw)) + 
    geom_boxplot(color = c("chartreuse4", "darkorchid2", "darkgoldenrod3"), lwd= 0.5)

cenplot.uw
cenUw <-cenplot.uw + scale_x_discrete(limits = level_order) + theme_bw() +
    labs(y = "Distance to centroid") + theme(axis.text.x = element_text(size=11, colour = "black"),
                                             axis.text.y = element_text(size=11, colour = "black"),
                                             axis.title.x = element_blank(),
                                             axis.title.y = element_text(size=14),
                                             panel.grid.major = element_blank(), 
                                             panel.grid.minor = element_blank(),
                                             panel.border = element_blank()) +
    geom_jitter(width=0.25, alpha=0.5) +
    theme(axis.line = element_line(color = 'black')) + 
    ggtitle("Unweighted")

cenUw

#Weighted
cenplot.we<-ggplot(centbrom, aes(x = Species, y= dist.we)) + 
    geom_boxplot(color = c("chartreuse4", "darkorchid2", "darkgoldenrod3"), lwd= 0.5)
cenplot.we

cenWe <-cenplot.we + scale_x_discrete(limits = level_order) + theme_bw() +
    labs(y = "Distance to centroid") + theme(axis.text.x = element_text(size=11, colour = "black"),
                                             axis.text.y = element_text(size=11, colour = "black"),
                                             axis.title.x = element_blank(),
                                             axis.title.y = element_text(size=14),
                                             panel.grid.major = element_blank(), 
                                             panel.grid.minor = element_blank(),
                                             panel.border = element_blank()) +
    geom_jitter(width=0.25, alpha=0.5) +
    theme(axis.line = element_line(color = 'black')) + 
    ggtitle("Weighted")
cenWe

##Make a panel with the two plots 
ggarrange(cenUw, cenWe,
          labels = c("A", "B"),
          ncol = 1, nrow = 2, align= "v")
