#################################################################################################################################
#Section 2A - Hits visualizations (Beeswarm, UpSet)
#################################################################################################################################
#following section utilizes lists of hits (gain or loss) generated from previous code 
#The following codes are used across multiple backgrounds and inform mutants that show phenotype in multiple background
library(plyr)
library(dplyr)
library(tibble)
library(viridisLite)
library(UpSetR)
library(gplots)
library(ggplot2)
library(beeswarm)

filepath <- "~/"
setwd(filepath) 
#function to column bind regardless of column length
cbind.fill <- function(...){
  nm <- list(...) 
  nm <- lapply(nm, as.matrix)
  n <- max(sapply(nm, nrow)) 
  do.call(cbind, lapply(nm, function (x) 
    rbind(x, matrix(, n-nrow(x), ncol(x))))) 
}
#
hist2d_breaks = function (x, y = NULL, nbins = 200,same.scale = FALSE, na.rm = TRUE, 
                          show = TRUE, col = c("black", heat.colors(12)), FUN = base::length, 
                          xlab, ylab,x.breaks,y.breaks, ...) 
{
  if (is.null(y)) {
    if (ncol(x) != 2) 
      stop("If y is ommitted, x must be a 2 column matirx")
    y <- x[, 2]
    x <- x[, 1]
  }
  if (length(nbins) == 1) 
    nbins <- rep(nbins, 2)
  nas <- is.na(x) | is.na(y)
  if (na.rm) {
    x <- x[!nas]
    y <- y[!nas]
  }
  else stop("missinig values not permitted if na.rm=FALSE")
  if(same.scale){
    x.cuts = x.breaks;
    y.cuts = x.breaks;
  }else{
    x.cuts <- x.breaks
    y.cuts <- y.breaks   
  }
  
  
  index.x <- cut(x, x.cuts, include.lowest = TRUE)
  index.y <- cut(y, y.cuts, include.lowest = TRUE)
  m <- tapply(x, list(index.x, index.y), FUN)
  if (identical(FUN, base::length)) 
    m[is.na(m)] <- 0
  if (missing(xlab)) 
    xlab <- deparse(substitute(xlab))
  if (missing(ylab)) 
    ylab <- deparse(substitute(ylab))
  if (show){
    image(x.cuts, y.cuts, m, col = col, xlab = xlab, ylab = ylab, 
          ...)
  } 
  midpoints <- function(x) (x[-1] + x[-length(x)])/2
  retval <- list()
  retval$counts <- m
  retval$counts_rel <- m/max(m)  
  retval$x.breaks = x.cuts
  retval$y.breaks = y.cuts
  retval$x = midpoints(x.cuts)
  retval$y = midpoints(y.cuts)
  retval$nobs = length(x)
  retval$bins = c(length(x.cuts)-1,length(y.cuts)-1)
  retval$call <- match.call()
  class(retval) <- "hist2d"
  retval
}


Loss_CAT217<-read.csv("~/.csv",row.names=1)
gainCutoff<-c(PAS217 = 1.1821, PAS332= 2.0451, SB1938=1.4362) #these values were determined from the gridmut/par distribution plots
lossCutoff<-c(PAS217 = 2.5532, PAS332= 1.6635, PAS482 = 6.6163, SB1938=1.1140)

#extract hit from phenotype category by Gridmut/par using calculated cutoffs
data217loss<-Loss_CAT217 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= lossCutoff["PAS217"] )

data482loss<-Loss_CAT482 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= lossCutoff["PAS482"] )

data332loss<-Loss_CAT332 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= lossCutoff["PAS332"] ) 

data1938loss<-Loss_CAT1938 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= lossCutoff["SB1938"] ) 

data217gain<-Gain_CAT217 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= gainCutoff["PAS217"] ) 

data332gain<-Gain_CAT332 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= gainCutoff["PAS332"] ) 

data1938gain<-Gain_CAT1938 %>% 
  rownames_to_column('Name') %>%
  filter(FoldChange >= gainCutoff["SB1938"] ) 

################################################################################################
# combine into gain and loss lists
genes_Loss<-list(PAS217=na.omit(data217loss$Name), 
                 PAS332=na.omit(data332loss$Name), 
                 PAS482=na.omit(data482loss$Name), 
                 SB1938=na.omit(data1938loss$Name))

genes_Gain<-list(PAS217=na.omit(data217gain$Name), 
                 PAS332=na.omit(data332gain$Name), 
                 SB1938=na.omit(data1938gain$Name))
# combine into single nucleation backgrounds, gain and loss lists
genes_Loss_noWT<-genes_Loss<-list(PAS332=na.omit(data332loss$Name), 
                                  PAS482=na.omit(data482loss$Name), 
                                  SB1938=na.omit(data1938loss$Name))

genes_Gain_noWT<-list(PAS332=na.omit(data332gain$Name), 
                 SB1938=na.omit(data1938gain$Name))
# combine into single nucleated loci, MAT v ECT, gain and loss lists
genes_Loss_singlenuc<-genes_Loss<-list(MAT=na.omit(c(data332loss$Name,data482loss$Name)), 
                                  ECT=na.omit(data1938loss$Name))

genes_Gain_singlenuc<-list(MAT=na.omit(data332gain$Name), 
                      ECT=na.omit(data1938gain$Name))
################################################################################################
# 2B upset plots by phenotype
################################################################################################
#plasma colors
show_col(viridis_pal(option = "plasma")(20))
cols<-viridis(n=20,option = "plasma")

pdf(file = "upsetPlots_plasmaLoss.pdf", width = 6, height = 6)
upset(fromList(genes_Loss), order.by = "freq",
      mainbar.y.label = "shared genes - Loss", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

pdf(file = "upsetPlots_plasmaGain.pdf", width = 4, height = 6)
upset(fromList(genes_Gain), order.by = "freq",
      mainbar.y.label = "shared genes - Gain", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

pdf(file = "upsetPlots_plasmaLoss_noWT.pdf", width = 6, height = 6)
upset(fromList(genes_Loss_noWT), order.by = "freq",
      mainbar.y.label = "shared genes - Loss", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

pdf(file = "upsetPlots_plasmaGain_noWT.pdf", width = 4, height = 6)
upset(fromList(genes_Gain_noWT), order.by = "freq",
      mainbar.y.label = "shared genes - Gain", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

pdf(file = "upsetPlots_plasmaLoss_singlenuc.pdf", width = 6, height = 6)
upset(fromList(genes_Loss_singlenuc), order.by = "freq",
      mainbar.y.label = "shared genes - Loss", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

pdf(file = "upsetPlots_plasmaGain_singlenuc.pdf", width = 4, height = 6)
upset(fromList(genes_Gain_singlenuc), order.by = "freq",
      mainbar.y.label = "shared genes - Gain", 
      sets.x.label = "enriched genes per strain",
      sets.bar.color = cols[1],
      main.bar.color = cols[3],
      matrix.color = cols[1],
      point.size = 3,
      text.scale = c(1.5, 1.5, 1.3, 1.1, 1.2, 1.8))
dev.off()

###################################################################################################
# Section 2C - Beeswarmplots with overlap annotations
###################################################################################################
# get the genes in overlaps
# make a venn diagram of all combos of intersections of genes above 2SD in the "on" cateogry
LossVenn<-venn(genes_Loss_noWT, show.plot=FALSE, intersections = TRUE)
# check that the counts for intersections match the counts in upset plots
lengths(attributes(LossVenn)$intersections)
# get the name of the intersections, should be strain names
names(lengths(attributes(LossVenn)$intersections))
# outupt the intersections aka what genes make up the counts in the upset or venn 
LossAttributes<-(attributes(LossVenn)$intersections) # this is a list
# convert the above list to a dataframe
OverlapsLoss<-data.frame(sapply(LossAttributes, "length<-", max(lengths(LossAttributes))))

# do as above for the Off category
GainVenn<-venn(genes_Gain_noWT, show.plot=FALSE, intersections = TRUE)
lengths(attributes(GainVenn)$intersections)
names(lengths(attributes(GainVenn)$intersections))
GainAttributes<-(attributes(GainVenn)$intersections) # this is a list
OverlapsGain<-data.frame(sapply(GainAttributes, "length<-", max(lengths(GainAttributes))))

#Following script trasnformed the above dataframe into a sinhgle comlumn indicating the number of backgrounds overlaped by mutant phenotype
LossSingle<-c(as.character(OverlapsLoss$'PAS332'),
              as.character(OverlapsLoss$'PAS482'),
              as.character(OverlapsLoss$'SB1938'))
LossSingle<-LossSingle[!is.na(LossSingle)]
LossDouble<-c(as.character(OverlapsLoss$'PAS332.PAS482'),
              as.character(OverlapsLoss$'PAS332.SB1938'),
              as.character(OverlapsLoss$'PAS482.SB1938'))
LossDouble<-LossDouble[!is.na(LossDouble)]
LossTriple<-as.character(OverlapsLoss$'PAS332.PAS482.SB1938')
LossTriple<-LossTriple[!is.na(LossTriple)]
OverlapsLossComb<-qpcR:::cbind.na(LossSingle,LossDouble,LossTriple)

#do so for gain of spreading phenotype
GainSingle<-c(as.character(OverlapsGain$'PAS332'),
              as.character(OverlapsGain$'SB1938'))
GainSingle<-GainSingle[!is.na(GainSingle)]
GainDouble<-as.character(OverlapsGain$'PAS332.SB1938')
OverlapsGainComb<-qpcR:::cbind.na(GainSingle,GainDouble,GainTriple)

#a function that bind the generated column of BG overlapped into the hit category
overlap<-function(x,phenotype){
  for (i in 1:nrow(x)) {
    x[i,'Overlap']<-colnames(phenotype)[which(phenotype == paste0(x[i,'Name']), arr.ind=T)[1,2]]}
}

data217loss_all<-Loss_CAT217
data217loss_all[,'Name']<-rownames(Loss_CAT217)

#process each phenotype category using the above function
for (i in 1:nrow(data217loss_all)) {
  data217loss_all[i,'Overlap']<-colnames(OverlapsLossComb)[which(OverlapsLossComb == paste0(data217loss_all[i,'Name']), arr.ind=T)[1,2]]}

CountPalette <- plasma(5)
data217loss_all[,"Background"]<-"PAS217"
for(i in 1:nrow(data217loss_all)){
  if (is.na(data217loss_all[i,'Overlap']) == TRUE) {
    data217loss_all[i,"OverlapNo"]<-"#DEDEDE" #color for non hits
  } else { if (data217loss_all[i,'Overlap']=="LossSingle") {
      data217loss_all[i,"OverlapNo"]<-CountPalette[1] #color for single background hit
  } else { if (data217loss_all[i,'Overlap']=="LossDouble") {
    data217loss_all[i,"OverlapNo"]<-CountPalette[2] #double
  } else { if (data217loss_all[i,'Overlap']=="LossTriple") {
    data217loss_all[i,"OverlapNo"]<-CountPalette[3] #triple
  } else { if (data217loss_all[i,'Overlap']=="LossQuadruple") {
    data217loss_all[i,"OverlapNo"]<-CountPalette[4] #quadruple
  }}}}}}
for(i in 1:nrow(data217loss_all)){
if (data217loss_all[i,'Name']== "Parent Iso1"||
    data217loss_all[i,'Name']== "Parent Iso2"||
    data217loss_all[i,'Name']== "Parent Iso3"||
    data217loss_all[i,'Name']== "Parent Iso4"||
    data217loss_all[i,'Name']== "Parent Iso5"||
    data217loss_all[i,'Name']== "Parent Iso6") {
  data217loss_all[i,"OverlapNo"]<-'#000000' #color for parent isolates
} else {next
    }}

#beeswarm plots
pdf("BeeswarmsLoss217.pdf", width = 3, height = 3, useDingbats = F)
ggplot(data= data217loss_all,mapping=aes(Background, FoldChange)) + 
  geom_beeswarm(aes(),col=data217loss_all$OverlapNo,size=0.75,cex=2.2,dodge.width=0.001) +
  geom_text_repel(data= subset(data217loss_all, Overlap==c('LossTriple','LossDouble')), #label the hits that are shared across 3 or 4 backgrounds
                  aes(label = Name),point.padding= unit(0.3, "lines"), box.padding = unit(0.3, "lines"),nudge_x = 0.1)+
  geom_text_repel(data= head(data217loss_all,n=10), #label top 10 hits
                 aes(label = Name),point.padding= unit(0.3, "lines"), box.padding = unit(0.3, "lines"),nudge_x = 0.1)+
  geom_hline(yintercept = sdcutoff_loss_217, col='red',size=0.25) + #sd cutoff
  #geom_text(aes(x= 1.3,y=sdcutoff, label=paste("\n",round(sdcutoff,digits = 4))),family= 'Helvetica', size = 3) +
  geom_hline(yintercept = as.numeric(percentilecutoff_loss_217), linetype = 2, col='blue',size=0.25) + #85th percentile cutoff
  #geom_text(aes(x= 1.3,y=percentilecutoff, label=paste("\n",round(percentilecutoff,digits = 4))),family= 'Helvetica', size = 3) +
  #scale_fill_manual(name = "Loss", values = c("#0D0887FF", "#7E03A8FF","#CC4678FF","#F89441FF"),labels = c("Single", "Double", "Triple", "Quadruple"),guide = guide_legend(fill = NULL,colour = NULL))+
  theme_classic() 
dev.off()