## Figure 4 Code...
rm(list=ls())
setwd("/data/Users/kmccauley/MUPPITS/Microbiome_Manuscript/Manuscript/Figure4/")
library(limma)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(lmerTest)
library(emmeans)
library(gridExtra)
theme_set(theme_bw())

get_cis <- function(model) {
  OR= exp(summary(model)$coef[2,1])
  LCL = exp(summary(model)$coef[2,1] - 2*summary(model)$coef[2,2])
  UCL = exp(summary(model)$coef[2,1] + 2*summary(model)$coef[2,2])
  return(c(OR=OR, LCL=LCL, UCL=UCL))
}

load("/data/Users/kmccauley/MUPPITS/DataFiles/totalNasalVoom.RData")
load("/data/Users/kmccauley/MUPPITS/DataFiles/nasalVoom0.RData")
modules <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/nasalGeneSets_374samples.csv", row.names=1)
fungal <- read.table("/data/Users/kmccauley/MUPPITS/OTUtables/Fungal OTU Tables/fungal/otu_table_reprare1000_nr.txt", header=TRUE, check.names=F, sep="\t", comment="", row.names=1)
key <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/KeyFile_OneLinePerEvent.csv", row.names=1)
data <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/MUPPITS_Merged_Mapping_UCSF_withDomGen.csv", row.names = 1)
module_lookup <- read.csv("/data/Users/kmccauley/MUPPITS/DataFiles/Nasal module summary annotation.csv")
#modify module names to be consistent with data:
module_lookup$Module.name <- sub(".n","", module_lookup$Module.name)
module_lookup$Module.name[!grepl("mac", module_lookup$Module.name)] <- sub("^m","unma", module_lookup$Module.name[!grepl("mac", module_lookup$Module.name)])
ls() # Confirm that the two objects made it into the environment  -- and confirm their names
str(nasalVoom0)
str(totalNasalVoom)
network_lookup <- data.frame(networks=paste0("Mod", 1:19), description=c("Strep/Haemophilus","Moraxella","Staphylococcus","Alloiococcus/Coryne","Enterococcus/Alicyclobacillus","Corynebacterium","Haemophilus","Fusobacterium","Neisseriaceae","Streptophyta","Prevotella","Acinetobacter","Leptotrichia","Corynebacterium","Moraxellaceae","Anaerococcus","Rothia mucilaginosa","Anaerococcus","Aggregatibacter"))
## Read in the microbiome data --
MBdata <- read.table("/data/Users/kmccauley/MUPPITS/OTUtables/MUPPITS_OTUtable_initial.txt", header=TRUE, check.names=F, sep="\t", comment="", row.names = 1, skip=1)
MBdata_mod <- read.csv("/data/Users/kmccauley/MUPPITS/Microbiome_Manuscript/Manuscript/NetworkGeneration/ModuleAggregated_OTU_SparCCOnly_seed123.csv", row.names = 1)
MBdata$row.names <- NULL
domOTU <- t(data.frame(lapply(MBdata[,1:(ncol(MBdata)-1)], function(x) rownames(MBdata)[which.max(x)]), check.names = F))
data <- merge(data, domOTU, by=0)
rownames(data) <- data$Row.names
data$Row.names <- NULL
totalVoom <- totalNasalVoom$E ## This is the voom-transformed data
dim(totalVoom)

make_data <- function(visits, season=NULL, keep.dups = FALSE) {
  key_v1a <- key[key$Analysis_Visit %in% visits, ]
  if(!is.null(season)) {
    keep.data <- data[data$Season.of.Cold %in% season,]
  } else {
    keep.data <- data
  }
  if(keep.dups==FALSE) {
    dup.ids <- keep.data$Subject.Identifier.for.the.Study[duplicated(keep.data$Subject.Identifier.for.the.Study)]
    newData <- keep.data[!(keep.data$Subject.Identifier.for.the.Study %in% dup.ids & keep.data$Analysis.Visit %in% "Visit 2a"),] ## subset to individuals if there are duplicates
  } else {
    newData <- keep.data
  }
  newVoom <- modules %>% 
    rownames_to_column("modules") %>% 
    filter(!modules %in% c("neut0","lymp0","eos0","mac0","epi0","unma0","squa0")) %>% 
    column_to_rownames(var="modules")
  ## This is where I change the transcriptomic input (modules or totalVoom)
  colnames(newVoom) <- key_v1a$specimenid[match(colnames(newVoom), c(as.character(key_v1a$gex_nasal_id_bri)))]
  gex_withMB <- newVoom[, !is.na(colnames(newVoom))]
  sum(duplicated(colnames(gex_withMB)))
  
  MB_withGEX <- as.data.frame(t(MBdata_mod)) %>% 
    dplyr::select(matches(colnames(gex_withMB))) %>% 
    dplyr::select(matches(rownames(newData)))
  
  gex_withMB2 <- as.data.frame(gex_withMB) %>% 
    dplyr::select(matches(names(MB_withGEX))) %>% 
    dplyr::select(matches(rownames(newData)))
  
  combined_data_list <- list(GEX=t(gex_withMB2), MB=t(as.matrix(MB_withGEX)))
  anly.data <- newData[colnames(gex_withMB2), ]
  anly.data$viral_pos <- !anly.data$Viral.Type.at.Visit %in% "Non-viral"
  anly.data$Viral.Exac <- paste(anly.data$viral_pos, anly.data$Case.or.Control.Status.Full.Cohort)
  anly.data$Viral.Exac <- factor(anly.data$Viral.Exac, levels=c("FALSE Control","FALSE Case","TRUE Control","TRUE Case"), labels=c("V- E-","V- E+","V+ E-","V+ E+"))
  anly.data$Case.or.Control.Status.Full.Cohort <- relevel(anly.data$Case.or.Control.Status.Full.Cohort, ref="Control")
  anly.data$fungal_pos <- rownames(anly.data) %in% names(fungal)
  combined_data_list[["anly.data"]] <- anly.data
  return(combined_data_list)
}
table(data$Season.of.Cold) # 2=FALL, 3=Spring, 4=Summer; 5=Winter
all_anly <- make_data(visits=c("Visit 1a"), keep.dups = FALSE) ## 87 V1a visits have microbiome + transcriptome. 10 have microbiome only, some amount have transcriptome only but it's difficult to quantify.
## See "GEX_MB_Comb/VoomData_Microbiome_Modules.R"  (or something similar) if I need to go back to the full list of statistics


library(ggsci)
## Adding Subgroup because Sue wanted only the significant group to be shown, and the significant subgroup is different for both panels
interaction_plots <- function(module, network, mod_name, net_name, dataset=all_anly, adjmodel, orlabels, subgroup) {
  all_data <- cbind(dataset[["anly.data"]], dataset[["GEX"]], dataset[["MB"]])
  all_data$split_net <- all_data[,network] > median(all_data[,network])
  all_data$module <- all_data[,module]
  myform <- as.formula(paste0("Case.or.Control.Status.Full.Cohort ~ module*split_net", adjmodel))
  model <- glm(myform, data=all_data, family=binomial)
  
  ## Stratified
  myformA <- as.formula(paste0("Case.or.Control.Status.Full.Cohort ~ module ", adjmodel))
  model.ov <- glm(myformA, data=all_data, family=binomial)
  model.hi <- glm(myformA, data=all_data[all_data$split_net,], family=binomial)
  model.lo <- glm(myformA, data=all_data[!all_data$split_net,], family=binomial)
  print("Overall:")
  print(get_cis(model.ov))
  print("Above Median:")
  print(get_cis(model.hi))
  print("Below Median:")
  print(get_cis(model.lo))
  
  myform2 <- as.formula(paste0("module ~ split_net"))
  mylist <- list(module=seq(min(all_data[,module]), max(all_data[,module]), length.out = 100),split_net=c(FALSE, TRUE))
  

  mypal <- rev(pal_npg("nrc")(2))
  names(mypal) <- c("Low Abundance","High Abundance")
  emmip_res <- emmip(model, myform2, CIs=TRUE, at=mylist, cov.reduce=range, type="response", plotit=FALSE)
  emmip_res$split_net <- factor(emmip_res$split_net, labels=c("Low Abundance","High Abundance"))
  emmip_res <- subset(emmip_res, split_net %in% subgroup)
  pltFinal <- ggplot(emmip_res, aes(x=module, y=yvar, color=split_net, fill=split_net)) +
    geom_line() +
    geom_ribbon(aes(y=yvar, x=module, ymin=LCL, ymax=UCL, alpha=0.3)) +
    ylab("Probability of Exacerbation") +
    xlab(mod_name) +
    scale_color_manual("Bacterial Network Abundance", values=mypal, drop=FALSE) +
    scale_fill_manual("Bacterial Network Abundance", values=mypal, drop=FALSE) +
    annotate("text", x=3, y=c(0.95, 0.9), label=orlabels, color=mypal, size=3) +
    guides(alpha=FALSE) +
    expand_limits(y=c(0,1)) +
    ggtitle(net_name)
  pltFinal
  return(pltFinal)
}

partA <- interaction_plots("unma27","Mod1","SMAD3 Expression",expression(~italic("Streptococcus")~ "Network"), adjmodel="+Age.in.years + Season.of.Cold + Sex", dataset=all_anly, orlabels=c("OR=1.78, P=0.337","OR=14.71, P=0.018"), subgroup="High Abundance") + 
  labs(tag="a)") +
  theme(text=element_text(size=9))
partB <- interaction_plots("unma27","Mod3","SMAD3 Expression",expression(~italic("Staphylococcus")~ "Network"), adjmodel="+Age.in.years + Season.of.Cold + Sex", dataset=all_anly, orlabels=c("OR=39.17, P=0.008","OR=1.46, P=0.47"), subgroup="Low Abundance")+ 
  labs(tag="b)") +
  theme(text=element_text(size=9))

completefig <- ggpubr::ggarrange(partA, partB, nrow=1, common.legend = TRUE, legend="right")
ggsave("SMAD3_Interactions3.jpeg", completefig, device="jpeg", height=90, width=180, units="mm")