rm(list=ls())
library(phyloseq)
library(vegan)
library(DESeq2)
library(ggplot2)
library(RColorBrewer)

theme_set(theme_bw())

  source("/data/Users/kmccauley/MUPPITS/Archive/DESeqFunctions.R")
setwd("/data/Users/kmccauley/MUPPITS/Microbiome_Manuscript/Manuscript/Figure3/")
  rare.otu <- read.table("/data/Users/kmccauley/MUPPITS/OTUtables/MUPPITS_OTUtable_initial.txt", header=TRUE, check.names=FALSE, sep="\t", comment="", skip=1, row.names = 1)
  map <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/MUPPITS_Merged_Mapping_UCSF_Dec2018.csv", row.names = "specimenid")
  map$X <- NULL
  map$`#SampleID` <- rownames(map)
  tree <- read_tree("/data/Users/kmccauley/MUPPITS/OTUtables/otus_labeled_aligned_pfiltered_raw_otu_table.tre")
  
  visit_info <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/MUPPITS1_Timeline_Data.csv")
  
  muppits.phy <- make_phyloseq(rare.otu, map)
  muppits.phy <- merge_phyloseq(muppits.phy, phy_tree(tree)) ## Add in the phylogenetic tree
  muppits.phy
  sample_data(muppits.phy)$viral_ac <- factor(sample_data(muppits.phy)$Viral.Type.at.Visit %in% c("RV-A","RV-C"))
  sample_data(muppits.phy)$viral_pos <- factor(! sample_data(muppits.phy)$Viral.Type.at.Visit %in% c("Non-viral"))

table(sample_data(muppits.phy)$Analysis.Visit, sample_data(muppits.phy)$Season.of.Cold)

#First, if the genus is just "g", make it the family...
tax_table(muppits.phy)[,"Genus"][tax_table(muppits.phy)[,"Genus"] %in% "g"] <- tax_table(muppits.phy)[,"Family"][tax_table(muppits.phy)[,"Genus"] %in% "g"]
tax_table(muppits.phy)[,"Genus"][tax_table(muppits.phy)[,"Genus"] %in% "f"] <- tax_table(muppits.phy)[,"Order"][tax_table(muppits.phy)[,"Genus"] %in% "f"]

#Take out the first letter:
tax_table(muppits.phy) <- apply(tax_table(muppits.phy), 2, function(x) substr(x, 2, 100))

muppits.overall <- subset_samples(muppits.phy, !Case.or.Control.Status.Original %in% "")
muppits.v1 <- subset_samples(muppits.overall, Analysis.Visit %in% "Visit 1a")

muppits.binary <- muppits.phy
binary.dat <- apply(otu_table(muppits.binary),2, function(x) as.numeric(x > 0))
rownames(binary.dat) <- taxa_names(muppits.binary)
otu_table(muppits.binary) <- otu_table(binary.dat, taxa_are_rows = TRUE)
ord <- ordinate(muppits.binary, method="PCoA", distance="binary")
plot_ordination(muppits.phy, ord, type="taxa", color="Order")

sample_data(muppits.v1)$Fall <- sample_data(muppits.v1)$Season.of.Cold %in% "Autumn (Sep.-Nov.)"
sample_data(muppits.v1)$Winter <- sample_data(muppits.v1)$Season.of.Cold %in% "Winter (Dec.-Feb.)"
sample_data(muppits.v1)$Spring <- sample_data(muppits.v1)$Season.of.Cold %in% "Spring (Mar.-May)"
sample_data(muppits.v1)$Summer <- sample_data(muppits.v1)$Season.of.Cold %in% "Summer (Jun.-Aug.)"

muppits.v1.fall <- subset_samples(muppits.v1, Fall %in% TRUE)
  set.seed(100)
  print(vegan::adonis2(phyloseq::distance(muppits.v1.fall, method="wunifrac") ~ sample_data(muppits.v1.fall)$Case.or.Control.Status.Full.Cohort + sample_data(muppits.v1.fall)$viral_pos + sample_data(muppits.v1.fall)$Study.Site + sample_data(muppits.v1.fall)$Age.in.years + sample_data(muppits.v1.fall)$Sex, by="margin"))

muppits.v1.season <- subset_samples(muppits.v1, Winter %in% TRUE)
set.seed(123)
print(vegan::adonis2(phyloseq::distance(muppits.v1.season, method="unifrac") ~ sample_data(muppits.v1.season)$Case.or.Control.Status.Full.Cohort + 
                       sample_data(muppits.v1.season)$viral_pos + sample_data(muppits.v1.season)$Study.Site + sample_data(muppits.v1.season)$Age.in.years + 
                       sample_data(muppits.v1.season)$Sex, by="margin"))

  
  ## using Viral positivity as a variable instead of HRV-A or HRV-C
#phy.prune <- filter_taxa(muppits.v1.fall, function (x) {sum(x > 0)/length(x) > 0.1}, prune=TRUE)
my_formula <- as.formula("~ Case.or.Control.Status.Full.Cohort + Study.Site + viral_pos + Age.in.years + Sex")
deseqdat <- phyloseq_to_deseq2(muppits.v1.fall, my_formula)
deseqdat2 <- DESeq(deseqdat, sfType = "poscounts", fitType = "local")
caco.status <- print_res(deseqdat2, muppits.v1.fall, var = "Case.or.Control.Status.Full.Cohort", ref="Control", cont="Case", alpha=0.1)
caco.status$group <- "Exacerbation"
viral.status <- print_res(deseqdat2, muppits.v1.fall, var = "viral_pos", ref="FALSE", cont="TRUE", alpha=0.1)
viral.status$group <- "Viral Positive"

all_data.fall <- rbind(caco.status, viral.status)
all_data.fall$otus <- paste0(all_data.fall$Genus, " (", all_data.fall$OTUname, ")")

otus_lab <- sapply(
  strsplit(as.character(all_data.fall$otus), " "), 
  function(x) parse(text = paste0("italic('", x[1], "')~", x[2]))
)

mycols <- c(brewer.pal(8, "Set1"), "grey")
names(mycols) <- c("Moraxella","Staphylococcus","Haemophilus", "Corynebacterium", "Alloiococcus","Streptococcus", "Leptotrichia", "Acinetobacter", "Other")

viral_plt <- plot_res(deseqdat2, muppits.v1.fall, var = "viral_pos", ref="FALSE", cont="TRUE", alpha=0.1, 
         genus.cols = mycols, left="V-",right="V+") +
  labs(tag="b)", size="Average Read Count") +
  ggtitle("Viral Status") +
  theme(text=element_text(size=8)) + 
  guides(fill=guide_legend(override.aes=list(alpha=0.6, size=5)))
exac_plt <- plot_res(deseqdat2, muppits.v1.fall, var = "Case.or.Control.Status.Full.Cohort", ref="Control", cont="Case", alpha=0.1, 
         genus.cols = mycols, left="E-",right="E+") +
  labs(tag="a)", size="Average Read Count") +
  ggtitle("Exacerbation") +
  theme(text=element_text(size=8)) + 
  guides(fill=guide_legend(override.aes=list(alpha=0.6, size=5)))

library(scales)
ggsave("Fall_Enriched_Taxa.jpeg", ggpubr::ggarrange(exac_plt, viral_plt, common.legend = TRUE, legend="right"), device="jpeg", height=130, width=180, dpi=300, units="mm")
write.csv(caco.status, "Exacerbation_DESeq.csv")
write.csv(viral.status, "Viral_DESeq.csv")