Workflow Steps and Code Snippets

suppressPackageStartupMessages(library(optparse, warn.conflicts=F, quietly=T))

option_list <- list(
  make_option(c( "--ABC_DexLPS_all"),
              type="character",
              help="Path to file with all ABC predictions passing 0.02 (in DexLPS condition)"),
  make_option(c("--ABC_LPS_all"),
              type="character",
              help="Path to file with all ABC predictions passing 0.02 (in LPS condition)"),
  make_option(c("--contrast_DexVSDexLPS"),
              type="character",
              help="Path to annotated tsv file of DeSeq2 contrast of DexLPS vs LPS"),
  make_option(c("--chipseq_summits"),
              type="character",
              help="Path to summit file of IDR peaks"),
  make_option(c("--igv"),
              type="character",
              help="Path to IGV snapshot of genomic locus (created manually)")
)

opt <- parse_args(OptionParser(option_list=option_list))

options(stringsAsFactors = FALSE)

suppressPackageStartupMessages(library(dplyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(org.Mm.eg.db, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(DESeq2, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ggplot2, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ComplexHeatmap, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(TxDb.Mmusculus.UCSC.mm10.knownGene, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(topGO, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ggExtra, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(plyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(here, warn.conflicts=F, quietly=T))

#set defaults for ggplot2 figures
theme_update(panel.background = element_rect(fill = "transparent", colour = NA),
             plot.background = element_rect(fill = "transparent", colour = NA),
             legend.background = element_rect(fill = "transparent", colour = NA),
             legend.key = element_rect(fill = "transparent", colour = NA),
             text=element_text(size=8, family = "ArialMT", colour="black"),
             title=element_text(size=10, family="ArialMT", colour="black"),
             panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             axis.text = element_text(size=8, family="ArialMT", colour="black"),
             axis.line = element_line(colour="black"),
             axis.ticks = element_line(colour="black"))

#----------------------------------------------------------------------
##--NOTE: example code for additional plots can be found in abc_visualizations.Rmd
#----------------------------------------------------------------------

#----------------------------------------------------------------------
##----------- load data
#----------------------------------------------------------------------

chipseq_summits <- rtracklayer::import.bed(opt$chipseq_summits)

contrast_DexLPSvLPS <- read.delim(opt$contrast_DexVSDexLPS)
# rename first column for merge later
colnames(contrast_DexLPSvLPS)[1] <- "EnsemblID"

ABC_DexLPS_all <- read.delim(opt$ABC_DexLPS_all)
ABC_LPS_all <- read.delim(opt$ABC_LPS_all)
# remove version number from ensembl IDs
ABC_DexLPS_all$TargetGene = gsub("\\.[0-9]*$","",ABC_DexLPS_all$TargetGene)
ABC_LPS_all$TargetGene = gsub("\\.[0-9]*$","",ABC_LPS_all$TargetGene)


#----------------------------------------------------------------------
##----------- data exploration
#----------------------------------------------------------------------

# What do the distributions of ABC scores look like

ggplot()+
  geom_density(data=ABC_DexLPS_all %>% filter(!class=="promoter"), aes(x=ABC.Score, colour="DexLPS"))+
  geom_density(data=ABC_LPS_all %>% filter(!class=="promoter"), aes(x=ABC.Score, colour="LPS"))+
  scale_colour_manual("",
                      breaks = c("DexLPS", "LPS"),
                      values = c("#339966", "#0066CC"))

# How do the ABC scores of promoters compare to those of enhancers?

ggplot()+
geom_density(data=ABC_DexLPS_all %>% filter(class=="promoter"), aes(x=ABC.Score, colour="promoter"))+
geom_density(data=ABC_DexLPS_all %>% filter(class=="genic"), aes(x=ABC.Score, colour="genic"))+
geom_density(data=ABC_DexLPS_all %>% filter(class=="intergenic"), aes(x=ABC.Score, colour="intergenic"))+
scale_colour_manual("",
                    breaks = c("promoter", "genic", "intergenic"),
                    values = c("red", "green", "blue"))

# What distance is used to define sth as promoter?

ggplot()+
  geom_density(data=ABC_DexLPS_all, aes(x=distance, colour=class))+
  scale_x_continuous(trans = "log10")+
  scale_colour_manual("",
                      breaks = c("promoter", "genic", "intergenic"),
                      values = c("red", "green", "blue"))

ggplot()+
  geom_density(data=ABC_DexLPS_all %>% filter(class=="promoter"), aes(x=distance, colour=isSelfPromoter))+
  scale_x_continuous(trans = "log10")+
  scale_colour_manual("isSelfPromoter",
                      breaks = c("True", "False"),
                      values = c("red", "blue"))+
  labs(title="All promoter regions, split by whether they're promoters for their target gene")+
  theme(
    plot.title = element_text(size=12)
  )
#----------------------------------------------------------------------
# ------ ------  Number of enhancers per gene
#----------------------------------------------------------------------
# Are most genes regulated by a single enhancer region or by multiple? How many does the ABC model assign per gene?

plot_enhancers_per_gene <- function(){
  DexLPS_TargetGeneCounts <- plyr::count(ABC_DexLPS_all %>% filter(!class=="promoter"), vars="TargetGene")
  LPS_TargetGeneCounts <- plyr::count(ABC_LPS_all %>% filter(!class=="promoter"), vars="TargetGene")

  gg <- ggplot()+
    geom_histogram(data=DexLPS_TargetGeneCounts, aes(x=freq, fill="DexLPS", colour="DexLPS"), binwidth = 1,  alpha=0.2)+
    geom_histogram(data=LPS_TargetGeneCounts, aes(x=freq, fill="LPS", colour="LPS"), binwidth = 1, alpha=0.2)+
    scale_fill_manual("",
                      breaks = c("DexLPS", "LPS"),
                      values = c("#339966", "#0066CC")) +
    scale_colour_manual("",
                        breaks = c("DexLPS", "LPS"),
                        values = c("#339966", "#0066CC")) +
    guides(color = "none")+
    theme(legend.position = c(0.8,0.8))+
    labs(title=" ",
         x="# of enhancers per gene",
         y="counts")

  return(gg)
}

gg_enh_per_gene <- plot_enhancers_per_gene()
gg_enh_per_gene

print("We find an average of")
plyr::count(ABC_DexLPS_all %>% filter(!class=="promoter"), vars="TargetGene") %>% pull(freq) %>% mean()
print("enhancers per gene for the DexLPS condition and")
plyr::count(ABC_LPS_all %>% filter(!class=="promoter"), vars="TargetGene") %>% pull(freq) %>% mean()
print("in LPS.")

#----------------------------------------------------------------------
# ------ ------  Number of genes per enhancer
#----------------------------------------------------------------------
# ABC allows for multiple assignments (an individual enhancer assigned to more than one promoter)  

plot_genes_per_enhancer <- function(){

  DexLPS_EnhancerCounts <- plyr::count(ABC_DexLPS_all %>% filter(!class=="promoter"),vars="name")
  LPS_EnhancerCounts <- plyr::count(ABC_LPS_all %>% filter(!class=="promoter"),vars="name")

  gg_full <- ggplot()+
    geom_histogram(data=DexLPS_EnhancerCounts, aes(x=freq, fill="DexLPS",colour="DexLPS"), binwidth = 1, alpha=0.2)+
    geom_histogram(data=LPS_EnhancerCounts, aes(x=freq, fill="LPS",colour="LPS"), binwidth = 1, alpha=0.2)+
    scale_fill_manual("",
                      breaks = c("DexLPS", "LPS"),
                      values = c("#339966", "#0066CC")) +
    scale_colour_manual("",
                        breaks = c("DexLPS", "LPS"),
                        values = c("#339966", "#0066CC")) +
    guides(color = FALSE)+
    labs(title="Number of genes an individual enhancer is assigned to",
         x="# of genes per enhancer",
         y="counts")
  plot(gg_full)

  gg_zoom <- ggplot()+
    geom_histogram(data=DexLPS_EnhancerCounts, aes(x=freq, fill="DexLPS",colour="DexLPS"), binwidth = 1, alpha=0.2)+
    geom_histogram(data=LPS_EnhancerCounts, aes(x=freq, fill="LPS",colour="LPS"), binwidth = 1, alpha=0.2)+
    scale_fill_manual("",
                      breaks = c("DexLPS", "LPS"),
                      values = c("#339966", "#0066CC")) +
    scale_colour_manual("",
                        breaks = c("DexLPS", "LPS"),
                        values = c("#339966", "#0066CC")) +
    xlim(c(0,8))+
    guides(color = FALSE)+
    theme(legend.position = c(0.8,0.8))+
    labs(title=" ",
         x="# of genes per enhancer",
         y="counts")
  plot(gg_zoom)

  print("We find an average of")
  DexLPS_EnhancerCounts %>% pull(freq) %>% mean() %>% print()
  print("genes per enhancer in the DexLPS condition and")
  LPS_EnhancerCounts %>% pull(freq) %>% mean() %>% print()
  print("in LPS.")

  return(gg_zoom)
}

gg_genes_per_enhancer_zoom <- plot_genes_per_enhancer()


#----------------------------------------------------------------------
# ------ Comparing ABC Scores (including promoter regions) between the conditions
#----------------------------------------------------------------------

# merge regions from both conditions but also keep those that are only in one
# 
# add info no whether the target genes are DE or now and what direction they're changed

LPS_GR<- plyranges::as_granges(ABC_LPS_all, seqnames=chr)
DexLPS_GR <- plyranges::as_granges(ABC_DexLPS_all, seqnames=chr)

#----------------------------------------------------------------------
merge_LPSandDexLPS_ABCdata <- function(ABC_LPS_all_GR, ABC_DexLPS_all_GR){
  ABC_LPS_all_unique <- unique(ABC_LPS_all_GR)
  ABC_DexLPS_all_unique <- unique(ABC_DexLPS_all_GR)

  # and then figure out which ones overlap with one another
  overlap_ABC_all_unique <- ChIPpeakAnno::findOverlapsOfPeaks (ABC_LPS_all_unique,
                                                               ABC_DexLPS_all_unique,
                                                               connectedPeaks = "keepAll",
                                                               ignore.strand = TRUE)
  # we coerce the info of overlapping regions into a dataframe
  overlap_ABC_all_unique_df <- as.data.frame(overlap_ABC_all_unique$overlappingPeaks)

  # and clean up column names
  colnames(overlap_ABC_all_unique_df)<- gsub("ABC_LPS_all_unique...ABC_DexLPS_all_unique.","", colnames(overlap_ABC_all_unique_df))
  colnames(overlap_ABC_all_unique_df)<- gsub(".1", "", colnames(overlap_ABC_all_unique_df))
  colnames(overlap_ABC_all_unique_df)<- paste( c(rep("LPS_",27),rep("DexLPS_",27),rep("",2)), colnames(overlap_ABC_all_unique_df), sep="")

  # we use this overlap to define a new variable that holds the info of which regions of the conditions form a pair
  overlap_ABC_all_unique_df$pairID <- paste("pair", seq(1, nrow(overlap_ABC_all_unique_df)), sep="" )

  #-----------------for DexLPS-------------------
  # merge the pairID onto the original dataframes
  ABC_DexLPS_all_wpairID <- merge(ABC_DexLPS_all,
                                  overlap_ABC_all_unique_df[,c("DexLPS_name","pairID")],
                                  by.x="name",
                                  by.y="DexLPS_name",
                                  all.x=TRUE)
  # for those that don't correspond to a pair, we assign condition and name as pairID
  ABC_DexLPS_all_wpairID <- ABC_DexLPS_all_wpairID %>% 
    mutate(pairID = case_when(is.na(pairID) ~ paste("DexLPS",name, sep="_"),
                              TRUE ~ pairID) )
  # combination of pairID and gene they're assigned to forms the assignment variable.
  # we will use this later to plot values from the 2 conditions, that have the same assignment
  ABC_DexLPS_all_wpairID$assignment <- paste(ABC_DexLPS_all_wpairID$pairID,
                                             ABC_DexLPS_all_wpairID$TargetGene,
                                             sep="_")

  #-----------------for LPS-------------------
  # merge the pairID onto the original dataframes
  ABC_LPS_all_wpairID <- merge(ABC_LPS_all,
                               overlap_ABC_all_unique_df[,c("LPS_name","pairID")],
                               by.x="name",
                               by.y="LPS_name",
                               all.x=TRUE)
  # for those that don't correspond to a pair, we assign condition and name as pairID
  ABC_LPS_all_wpairID <- ABC_LPS_all_wpairID %>% 
    mutate(pairID = case_when(is.na(pairID) ~ paste("LPS",name, sep="_"),
                              TRUE ~ pairID) )

  ABC_LPS_all_wpairID$assignment <- paste(ABC_LPS_all_wpairID$pairID,
                                          ABC_LPS_all_wpairID$TargetGene,
                                          sep="_")

  merged_conditions <- merge(ABC_DexLPS_all_wpairID,
                                 ABC_LPS_all_wpairID,
                                 by="assignment", 
                                 all=TRUE,
                                 suffixes = c(".DexLPS", ".LPS"))
  merged_conditions <- merged_conditions %>%
    mutate(TargetGene = case_when(TargetGene.DexLPS==TargetGene.LPS ~ TargetGene.DexLPS,
                                  is.na(TargetGene.DexLPS) ~ TargetGene.LPS,
                                  is.na(TargetGene.LPS) ~ TargetGene.DexLPS,
                                  TRUE ~ NA_character_)
    )

  return(merged_conditions)
}

merged_conditions_ABC <- merge_LPSandDexLPS_ABCdata(LPS_GR,DexLPS_GR)

#----------------------------------------------------------------------
# ------ merge on gene expression results
#----------------------------------------------------------------------

# merge the RNAseq results to use as colour later
merged_conditions_ABC <- merge(merged_conditions_ABC, contrast_DexLPSvLPS, 
                               by.x="TargetGene", by.y="EnsemblID")
merged_conditions_ABC <- merged_conditions_ABC %>%
  mutate(change=case_when(padj<0.05 & log2FoldChange>0.58 ~ "up",
                          padj<0.05 & log2FoldChange<(-0.58) ~ "down",
                          TRUE ~ "no change"))

# replace the na with 0
merged_conditions_ABCnomissing <- merged_conditions_ABC %>% 
  mutate(ABC.Score.DexLPS = case_when(is.na(ABC.Score.DexLPS) ~ 0,
                                      TRUE ~ ABC.Score.DexLPS) ) %>%
  mutate(ABC.Score.LPS = case_when(is.na(ABC.Score.LPS) ~ 0,
                                   TRUE ~ ABC.Score.LPS) )

cor_ABCscores <- cor(merged_conditions_ABCnomissing$ABC.Score.DexLPS,
                     merged_conditions_ABCnomissing$ABC.Score.LPS)

gg_merged_conditions_ABCnomissing <- ggplot(merged_conditions_ABCnomissing)+
  geom_point(aes(x=ABC.Score.DexLPS, y=ABC.Score.LPS, fill=log2FoldChange), alpha=0.2, size=1, stroke = 0, shape=21)+
  scale_fill_gradient(low="blue", high="red")+
  ylim(c(0,0.8))+
  xlim(c(0,0.8))+
  geom_text(x=0.2, y=0.75, label=paste("r = ",format(cor_ABCscores, digits = 2)))+
  labs(fill="log2FC", x="ABC score Dex+LPS", y="ABC score LPS")+
  theme(legend.position = "top")

gg_merged_conditions_ABCnomissing

#----------------------------------------------------------------------
# ------ plot foldchange in ABC score vs foldchange in expression
#----------------------------------------------------------------------

cor_ABCdiff_log2FC <- with(merged_conditions_ABCnomissing %>% filter(change!="no change"),
                           cor((ABC.Score.DexLPS - ABC.Score.LPS),log2FoldChange))

gg_ABCdiff_log2FC <- ggplot(merged_conditions_ABCnomissing %>% filter(change!="no change"), aes(x=(ABC.Score.DexLPS - ABC.Score.LPS), y=log2FoldChange)) +
  geom_hex(bins=70)+
  viridis::scale_fill_viridis()+
  geom_smooth(method = "lm", colour="black")+
  geom_text(x=0.6, y=5, label=paste("r = ",format(cor_ABCdiff_log2FC, digits = 2)))+
  labs(x="ABC score Dex+LPS - ABC score LPS", y="expression log2FC(Dex+LPS / LPS)")+
  theme(legend.position = "top")



#----------------------------------------------------------------------
# ------ ABC score by peak presence
#----------------------------------------------------------------------

plot_ABC_with_marginals <- function(whatchange){
  gg <- ggplot(data= merged_conditions_ABCnomissing %>% filter(change==whatchange) , 
               aes(x=ABC.Score.DexLPS, y=ABC.Score.LPS ) )+
    geom_point( aes( fill=whatchange) , size=1, alpha=1,stroke = 0, shape=21) +
    geom_abline(slope = 1)+
    ylim(c(0,0.8))+
    xlim(c(0,0.8))+
    labs(x="ABC score Dex+LPS", y="ABC score LPS")+
    scale_fill_manual(
      breaks=c("up","no change","down"),
      values=c("red","black","blue") )+
    theme(legend.position = "none")

  gg_marg <- ggExtra::ggMarginal(gg, type="histogram", size=2,
                                 xparams = list(bins=85), yparams = list(bins=85))

  return(gg_marg)
}

ggplot_ABC_with_marginals_up <- plot_ABC_with_marginals("up")
ggplot_ABC_with_marginals_up
ggplot_ABC_with_marginals_down <-plot_ABC_with_marginals("down")
ggplot_ABC_with_marginals_down

#----------------------------------------------------------------------
# ------ show marginals as bar plot
#----------------------------------------------------------------------
gg_marginals_dexlps <- 
  ggplot()+
  geom_histogram(data=merged_conditions_ABCnomissing %>% filter(change=="up"), aes(x=ABC.Score.DexLPS, fill="up"), alpha=0.2, bins=85)+
  geom_histogram(data=merged_conditions_ABCnomissing %>% filter(change=="down"), aes(x=ABC.Score.DexLPS, fill="down"), alpha=0.2, bins=85)+
  scale_fill_manual(
    breaks=c("up","down"),
    values=c("red","blue") )+
  xlim(c(NA,0.3))+
  labs(fill="gene expression change",x="ABC score DexLPS", y="ABC score LPS")+
  theme(legend.position = c(0.7,0.7))
gg_marginals_lps <- 
  ggplot()+
  geom_histogram(data=merged_conditions_ABCnomissing %>% filter(change=="up"), aes(x=ABC.Score.LPS, fill="up"), alpha=0.2, bins=85)+
  geom_histogram(data=merged_conditions_ABCnomissing %>% filter(change=="down"), aes(x=ABC.Score.LPS, fill="down"), alpha=0.2, bins=85)+
  scale_fill_manual(
    breaks=c("up","down"),
    values=c("red","blue") )+ 
  xlim(c(NA,0.3))+
  labs(fill="gene expression change",x="ABC score DexLPS", y="ABC score LPS")+
  theme(legend.position = c(0.7,0.7))

#----------------------------------------------------------------------
# ------ ABC score by peak presence
#----------------------------------------------------------------------
# Is there a difference in ABC score between the enhancers overlapping with a GR summit vs the ones that don't?


DexLPS_GR_summits_overlap <- IRanges::findOverlaps(DexLPS_GR,chipseq_summits) #(query,subject)
DexLPS_GR$haspeakID <- FALSE
DexLPS_GR$haspeakID[queryHits(DexLPS_GR_summits_overlap)] <- TRUE


plot_ABCscore_ofDEgenes_byGRoverlap <- function(){
  DE_ENSEMBL <- contrast_DexLPSvLPS %>%  filter(padj<0.05 & abs(log2FoldChange)>0.58 ) %>% pull(EnsemblID)

  DexLPS_enhancers_DEsubset <- 
    as.data.frame(DexLPS_GR) %>% 
    filter(TargetGene %in% DE_ENSEMBL) %>%
    filter(!class=="promoter")

  DexLPS_enhancers_DEsubset %>% filter(haspeakID==TRUE) %>% dplyr::pull(ABC.Score) %>% mean() %>% print()
  DexLPS_enhancers_DEsubset %>% filter(haspeakID==FALSE) %>% dplyr::pull(ABC.Score) %>% mean() %>% print()

  # show distribution for the two groups

  gg1 <- ggplot()+
    geom_density(data=DexLPS_enhancers_DEsubset %>% filter(haspeakID==TRUE), 
                 aes(x=ABC.Score, colour="haspeakID", fill="haspeakID"), alpha=0.3)+
    geom_density(data=DexLPS_enhancers_DEsubset %>% filter(haspeakID==FALSE), 
                 aes(x=ABC.Score, colour="nopeakID", fill="nopeakID"), alpha=0.3)+
    scale_x_continuous(trans = "log2")+
    scale_colour_manual("",
                        labels = c("has GR peak", "no GR peak"),
                        breaks = c("haspeakID", "nopeakID"),
                        values = c("darkmagenta", "cadetblue"))+
    scale_fill_manual("",
                      labels = c("has GR peak", "no GR peak"),
                      breaks = c("haspeakID", "nopeakID"),
                      values = c("darkmagenta", "cadetblue"))+
    guides(colour="none")+
    #xlim(c(-3,0))+
    labs(title=" ",
         x="ABC score Dex+LPS",
         y="density")+
    theme(
      axis.text = element_text(size=12, family="ArialMT", colour="black"),
      axis.title = element_text(size=12, family="ArialMT", colour="black"),
      axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
      legend.position = c(0.8,0.8)
    )
  plot(gg1)

  gg2 <- DexLPS_enhancers_DEsubset %>%
    ggplot( aes(x=hic_contact_pl_scaled_adj, y=log2(activity_base))) +
    geom_hex()+
    facet_wrap(~haspeakID, labeller="label_both")+
    viridis::scale_fill_viridis()+
    labs(title = " ", x="Hi-C contact", y="log2(base activity)")

  plot(gg2)
  all_plots=list()
  all_plots[[1]] <- gg1
  all_plots[[2]] <- gg2
  return(all_plots)
}

gg_ABCscore_ofDEgenes_byGRoverlap <- plot_ABCscore_ofDEgenes_byGRoverlap()


#----------------------------------------------------------------------
# ------ load IGV plot
#----------------------------------------------------------------------
igv <- png::readPNG( opt$igv )

gg_igv <- ggplot() + 
  ggpubr::background_image(igv) +
  # so it doesn't get squished
  #coord_fixed()+
  # This ensures that the image leaves some space at the edges
  theme(plot.margin = margin(t=0, l=0, r=0, b=0, unit = "cm"),
        axis.line = element_blank())


#----------------------------------------------------------------------
# ------ put figure panel together
#----------------------------------------------------------------------

# save giant files separately
ggsave(here("./results/current/Figures/Figure_abcresults_3B.bmp"), gg_merged_conditions_ABCnomissing,
       width=75, height=100, units="mm",
       bg="white")
ggsave(here("./results/current/Figures/Figure_abcresults_3B.png"), gg_merged_conditions_ABCnomissing,
       width=75, height=100, units="mm",
       bg="white")

gg_placeholder <- ggplot() + 
  theme(plot.margin = margin(t=0, l=0, r=0, b=0, unit = "cm"),
        axis.line = element_blank())

# first row
gg_r1c1 <-  ggpubr::ggarrange(gg_enh_per_gene, gg_genes_per_enhancer_zoom,
                            labels = c("A", NA),
                            ncol = 1, nrow = 2, heights = c(1,1))

gg_r1c2 <-  ggpubr::ggarrange(gg_placeholder, ggplot_ABC_with_marginals_up, ggplot_ABC_with_marginals_down, 
                            labels = c("B", "C", NA),
                            ncol = 3, nrow=1, widths=c(1,1,1))

gg_r1 <- ggpubr::ggarrange(gg_r1c1, gg_r1c2, 
                           labels = c(NA, NA),
                           ncol = 2, nrow=1, widths=c(1,3))

# second row
gg_r2 <-  ggpubr::ggarrange(gg_ABCdiff_log2FC, gg_ABCscore_ofDEgenes_byGRoverlap[[1]], gg_ABCscore_ofDEgenes_byGRoverlap[[2]], 
                            labels = c("D","E", "F"),
                            ncol = 3, nrow=1, widths=c(1,1,1.5))


gg_tophalf <- ggpubr::ggarrange(gg_r1, gg_r2, 
                                nrow=2, heights =c(1,1))
gg_tophalf

full_panel <- ggpubr::ggarrange(gg_tophalf, gg_igv, 
                                labels=c(NA,"G"),
                                nrow=2, heights = c(1.8,1))

#full_panel
ggsave(here("./results/current/Figures/Figure_abcresults.png"), full_panel,
       width=300, height=300, units="mm",
       bg="white")
ggsave(here("./results/current/Figures/Figure_abcresults.pdf"), full_panel,
       width=300, height=300, units="mm",
       bg="white")

suppressPackageStartupMessages(library(optparse, warn.conflicts=F, quietly=T))

option_list <- list(
  make_option(c("--log2fcthresh"),
              type="numeric",
              help="Log2FC threshold used in addition to adj.pval to define significant genes"),
  make_option(c("--chipseq_summits"),
              type="character",
              help="Path to summit file of IDR peaks"),
  make_option(c("--genekey"),
              type="character",
              help="Path to biomart genekey that mapps ensembl geeneIDs to MGI symbols"),
  make_option(c("--contrast_DexVSDexLPS"),
              type="character",
              help="Path to annotated tsv file of DeSeq2 contrast of DexLPS vs LPS"),
  make_option(c("--meme_db_path"),
              type="character",
              help="Path to JASPAR motif db file"),
  make_option(c( "--rna_nascent_fpkm"),
              type="character",
              help="FPKM matrix of 4sU experiment"),
  make_option(c("-o", "--outdir"),
              type="character",
              help="Path to output directory"))

opt <- parse_args(OptionParser(option_list=option_list))

# set output for logfile to retrieve stats for plot later
sink(file=paste0(opt$outdir,"figure_proxanno_prepdata.out"))

suppressPackageStartupMessages(library(memes, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(universalmotif, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(biomaRt, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(TxDb.Mmusculus.UCSC.mm10.knownGene, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm10, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(dplyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ChIPseeker, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(stringr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(plyranges, warn.conflicts=F, quietly=T))

#-------------------------------
## Import references
#-------------------------------
# for gene annotation
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
# for the sequence
# either use masked or unmasked (the mask does NOT seem to be for repeats though)
mm.genome <- BSgenome.Mmusculus.UCSC.mm10

#-------------------------------
### Determine expressed genes
#-------------------------------
rna_nascent <- read.table(opt$rna_nascent_fpkm, header=TRUE)
print("Determine expressed genes using 4sU data")
# what does the read count distribution look like?
# compute median per gene and plot it as histogram
rna_nascent$median_genecounts <- apply(rna_nascent[,-1], 1, FUN=median)
# lots of medians that are below 1
hist(log10(rna_nascent$median_genecounts))
# filter based on expression
expressed_genes <- rna_nascent %>% 
  dplyr::filter(median_genecounts > 0) 

#-------------------------------
### Load genekey to annotate ensembl to mgi
#-------------------------------
geneKey <- read.delim(opt$genekey)

#merge gene annotations to results table
expressed_genes <- merge(expressed_genes,
                         geneKey, 
                         by.x="Geneid", 
                         by.y="ensembl_gene_id")

#-------------------------------
## Getting the sequences
#-------------------------------
# import summit of ChIPseq peaks
ChIPseq_summits <- read.table(opt$chipseq_summits)
ChIPseq_ranges <- GRanges(seqnames = ChIPseq_summits[,c("V1")],
                          ranges = IRanges(start=ChIPseq_summits[,c("V2")],
                                           end=ChIPseq_summits[,c("V3")]-1)) # to make up for 0 vs 1 encoding
ChIPseq_ranges$id <- c(1:length(ChIPseq_ranges))


# NOTE: if we only want to annotate to genes that are expressed, we could use ChIPpeakAnno and a filtered annoDB object instead
# annotate it to genes
print("Annotate ChIPseq summit to closest gene (using genomic reference)")
summitAnno <- annotatePeak(ChIPseq_ranges, 
                           tssRegion=c(-3000, 3000), 
                           TxDb=txdb, annoDb = "org.Mm.eg.db")

summitAnno_df <- summitAnno %>% as.data.frame()

# see which ones are DE genes and add that info to GRanges as column "directionchange"
print("Add info on which genes are DE to the annotated summits")
DE_4sU <- read.delim(opt$contrast_DexVSDexLPS)

summitAnno_df <- left_join(summitAnno_df, 
                           DE_4sU[,c("mgi_symbol","padj","log2FoldChange")],
                           by = c("SYMBOL" = "mgi_symbol"))

summitAnno_df <- summitAnno_df %>% mutate(change = case_when(padj<0.05 & log2FoldChange > opt$log2fcthresh ~ "up",
                                                             padj<0.05& log2FoldChange < -opt$log2fcthresh ~"down",
                                                             TRUE ~ "ns")
)

# save info on gene annotation
ChIPseq_ranges$mgi_symbol[match(summitAnno_df$id, ChIPseq_ranges$id )] <- summitAnno_df$SYMBOL

# assign directionchange as metadata column
ChIPseq_ranges$directionchange[match(summitAnno_df$id, ChIPseq_ranges$id )] <- summitAnno_df$change

# ass distance to TSS
ChIPseq_ranges$distanceToTSS[match(summitAnno_df$id, ChIPseq_ranges$id )] <- summitAnno_df$distanceToTSS

#-------------------------------
## Prefilter motifdb to motifs that are expressed in celltype
#-------------------------------
print("Prefilter meme_db for those motifs expressed in our 4sU data")
meme_db <- read_meme(opt$meme_db_path) %>% 
  to_df()

meme_db_expressed <- meme_db %>% 
  # the altname slot of meme_db contains the gene symbol (this is database-specific)
  # avoid mismatches cased by casing and keep motif if at least one part of composite is expressed
  tidyr::separate(altname, into=c("tf1", "tf2"), sep="::",remove=FALSE) %>%
  filter( str_to_upper(tf1) %in% str_to_upper(expressed_genes$mgi_symbol) | str_to_upper(tf2) %in% str_to_upper(expressed_genes$mgi_symbol)) %>%
  # we don't need the split TF info downstream
  dplyr::select(!c("tf1","tf2"))


print("Number of motifs pre-filtering: ")
nrow(meme_db)
print("Number of motifs post-filter: ")
nrow(meme_db_expressed)

#-------------------------------
## OPTIONAL: only run with motifs of interest
#-------------------------------
meme_motifsOI <- 
  meme_db_expressed %>% 
  filter(
    grepl("STAT", str_to_upper(altname)) |
      grepl("NR3C", str_to_upper(altname))
  )

#-------------------------------
## FIGURES on peak gene annotation
#-------------------------------

# filter peaks for those annotated to genes that are expressed
summitAnno_expr <- subset(summitAnno,
                          summitAnno@anno$SYMBOL %in% expressed_genes$mgi_symbol)

# filter the df version in the same fashion
summitAnno_df_expr <- summitAnno_df %>% filter(SYMBOL %in% expressed_genes$mgi_symbol)

#---------------------------------
# --- some stats
#---------------------------------

distbygene <-  summitAnno_df_expr  %>% 
  group_by(SYMBOL, change) %>%
  summarise(min_dist=min(abs(distanceToTSS)), 
            mean_dist=mean(abs(distanceToTSS))) %>%
  ungroup() %>%
  mutate(logmindist=log2(min_dist+1))

# why 30kb cutoff
distbygene_allDE <- distbygene  %>% 
  filter(!change=="ns") %>%
  mutate(change=factor(change,levels=c("down","up")))

print("We need to justify why we picked a cutoff of 30kb.")
print("From a genecentric perspective, we want to include the peak regions that most likely have a regulating function on the gene.")

print("With a cutoff of 30kb, how many genes DONT have at least one peak within that range?")
tbl <- table(distbygene_allDE$min_dist > 30000)
tbl[2]/(tbl[1]+tbl[2])

print("How many genes do we lose of both sets by using that cutoff?")
print("In the upregulated fraction:")
table( (distbygene %>% filter(change=="up"))$min_dist > 30000)
print("In the downregulated fraction:")
table( (distbygene %>% filter(change=="down"))$min_dist > 30000)

print("Min and mean dist for the genes with log2FC >", opt$log2fcthresh)
distbygene %>%
  filter(change=="up") %>%
  summarise_all(mean) %>%
  print()
print("Min and mean dist for the genes with log2FC <", opt$log2fcthresh )
distbygene %>%
  filter(change=="down") %>%
  summarise_all(mean) %>%
  print()

print("How many peaks per UPregulated gene:")
summitAnno_df_expr  %>% 
  filter(change=="up")%>%
  filter(abs(distanceToTSS)<30000) %>%
  group_by(SYMBOL) %>%
  summarise(count=n()) %>%
  pull(count) %>% 
  mean()
print("How many peaks per DOWNregulated gene:")
summitAnno_df_expr  %>% 
  filter(change=="down")%>%
  filter(abs(distanceToTSS)<30000) %>%
  group_by(SYMBOL) %>%
  summarise(count=n()) %>%
  pull(count) %>% 
  mean()

print("The distances between the peaks mapping to the same gene.")
print("returns NA if only one 1 peak is annotated to the gene - those are excluded")
print("For upregulated genes:")
summitAnno_df_expr  %>% 
  filter(change=="up")%>%
  filter(abs(distanceToTSS)<30000) %>%
  group_by(SYMBOL) %>%
  summarise(meanpeakdist = mean(dist(distanceToTSS))) %>% # 
  filter(!is.na(meanpeakdist)) %>%
  pull(meanpeakdist) %>%
  mean()
print("For downregulated genes:")
summitAnno_df_expr  %>% 
  filter(change=="down")%>%
  filter(abs(distanceToTSS)<30000) %>%
  group_by(SYMBOL) %>%
  summarise(meanpeakdist = mean(dist(distanceToTSS))) %>% 
  filter(!is.na(meanpeakdist)) %>%
  pull(meanpeakdist) %>%
  mean()

#--------------------------------------
#- permutations
#--------------------------------------

# Difference in means
groupdiff <- diff(tapply(distbygene_allDE$min_dist, distbygene_allDE$change, mean))
print("The mean minimum distance is smaller for the upregulated set")
print(paste("The group difference is: ",groupdiff))
print("Permutation test to see if this difference between the groups is meaningful")

#Permutation test
permutation.test <- function(group, outcome, n, reference){
  distribution=c()
  result=0
  for(i in 1:n){
    distribution[i]=diff(by(outcome, sample(group, length(group), FALSE), mean))
  }
  result=sum(abs(distribution) >= abs(groupdiff))/(n)
  return(list(result, distribution, groupdiff))
}

permtest_res <- permutation.test(distbygene_allDE$change, distbygene_allDE$min_dist, 100000, groupdiff)

#--------------------------------------
#- export objects
#--------------------------------------

#---------------------------------
# --- export results of the permutation test
saveRDS(permtest_res,
        file=paste0(opt$outdir,"permtest_res.rds"))

#---------------------------------
# --- export up and downregulated summitfraction for deeptools

table(ChIPseq_ranges$directionchange)
dir.create(paste0(opt$outdir,"peaks_annot2DEgenes_30kb_log2FC0.58/"))
export.bed(ChIPseq_ranges %>% filter(directionchange == "up") %>% filter(abs(distanceToTSS)<30000) ,
           con=paste0(opt$outdir,"peaks_annot2DEgenes_30kb_log2FC0.58/UP_summit_unmerged.bed"))
export.bed(ChIPseq_ranges %>% filter(directionchange == "down") %>% filter(abs(distanceToTSS)<30000) ,
           con=paste0(opt$outdir,"peaks_annot2DEgenes_30kb_log2FC0.58/DOWN_summit_unmerged.bed"))

#-------------------------------
## export objects to run memes afterwards

saveRDS(ChIPseq_ranges,
        file=paste0(opt$outdir,"../memes_bioc/ChIPseq_summit_Granges.rds"))

saveRDS(meme_db_expressed,
        file=paste0(opt$outdir,"../memes_bioc/meme_db_4sUexpressed.rds"))

#-------------------------------
## export objects for ggplot figures
saveRDS(summitAnno,
        file=paste0(opt$outdir,"summitAnno.rds"))
saveRDS(summitAnno_expr,
        file=paste0(opt$outdir,"summitAnno_expr.rds"))
saveRDS(summitAnno_df_expr,
        file=paste0(opt$outdir,"summitAnno_df_expr.rds"))

sink()

suppressPackageStartupMessages(library(optparse, warn.conflicts=F, quietly=T))

option_list <- list(
  make_option(c("--ctss_pool1"), 
              type="character",
              help="Path for ctss file from FANTOM5 of pool 1"),
  make_option(c("--ctss_pool2"), 
              type="character",
              help="Path for ctss file from FANTOM5 of pool 2"),
  make_option(c("--liftoverchain"), 
              type="character",
              help="Path to liftover chain file for mm9 to mm10"),
  make_option(c("--gencode_mm9_geneanno"), 
              type="character",
              help="Path to genomic reference file for mm9"),
  make_option(c("--gencode_mm10_geneanno"), 
              type="character",
              help="GENCODE genomic reference for assembly mm10, prefiltered for gene entries"),
  make_option(c("--outdir"), 
              type="character",
              help="Output directory")
)
opt <- parse_args(OptionParser(option_list=option_list))

suppressPackageStartupMessages(library(dplyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm9, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ChIPseeker, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(CAGEr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(rtracklayer, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ggplot2, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(GenomicRanges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(here, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(org.Mm.eg.db, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(TxDb.Mmusculus.UCSC.mm10.knownGene, warn.conflicts=F, quietly=T))

outdir <- here(opt$outdir)

# Workflow is based on the CAGEr vignette
# https://www.bioconductor.org/packages/release/bioc/vignettes/CAGEr/inst/doc/CAGEexp.html


#----------------------   Import CAGE samples from BMDMs
#------------------------------------------------------------------------------------

# get URL for public samples, download and read necessariy columns into ctss file (see snakemake workflow)
# We want to import bone-marrow derived macrophage samples through CAGEr.
# After looking at list of available datasets we can decide what samples best fit our needs.
# Let's see what samples are available through FANTOM5
#data(FANTOM5mouseSamples)
#head(FANTOM5mouseSamples)
# The FANTOM5 dataframe holds descriptions of the samples and the url where they can be retrieved.
# There's an easy way to import samples that match a certain term into a CAGEset object.
#mac_samples <- FANTOM5mouseSamples[grep("macrophage, bone marrow derived",
#                                        FANTOM5mouseSamples[,"description"]),]

print("NOTE: reference genome for the public CAGE files is mm9!")

ce <- CAGEr::CAGEexp( genomeName = "BSgenome.Mmusculus.UCSC.mm9",
                      inputFiles = c(opt$ctss_pool1 ,opt$ctss_pool2),
                      inputFilesType = "ctss",
                      sampleLabels   = c("pool1","pool2")
)
# To actually read in the data into the object we use getCTSS() function, that will add an experiment called tagCountMatrix to the CAGEexp object.
ce <- CAGEr::getCTSS(ce)
ce

#------------------------------------------------------------------------------------
#----------------------   QC
#------------------------------------------------------------------------------------

ncbim37_anno <- rtracklayer::import.gff(opt$gencode_mm9_geneanno)

ce <- annotateCTSS(ce, ncbim37_anno)
colData(ce)[,c("librarySizes", "promoter", "exon", "intron", "unknown")]
plotAnnot(ce, "counts")

corr.m <- plotCorrelation2( ce, samples = "all"
                            , tagCountThreshold = 1, applyThresholdBoth = FALSE
                            , method = "pearson")

#------------------------------------------------------------------------------------
#----------------------   Get read clusters
#------------------------------------------------------------------------------------

print("Merging samples")
#Now we can merge them
ce <- mergeSamples(ce, mergeIndex = c(1,1), 
                   mergedSampleLabels = c("BMDM"))

# redo annotation since this gets reset during merging
ce <- annotateCTSS(ce, ncbim37_anno)

print("The total library size is:")
print(librarySizes(ce))

# Check if data follows a power law distribution
plotReverseCumulatives(ce, fitInRange = c(5, 3000), onePlot = TRUE)

print("Normalizing reads")
# Since we don't really care about making comparisons between different population we could prob just skip the normalization
# The fit range is chosen from the plot. We take the alpha from the ref distribution and set T to a million to get the tag count per million (TPM)
ce <- normalizeTagCount(ce, method = "powerLaw", fitInRange = c(5, 3000), alpha = 1.15, T = 1*10^6)
#mac_CAGEset@tagCountMatrix

print("Cluster the tags")
# After normalization we can cluster the tags.
# Clustering, only seems to work with the CAGEset object (due to some problems with the IRanges column)
# From the CAGEr vignette:
# "Transcription start sites are found in the promoter region of a gene and reflect the transcriptional activity of that promoter (Figure 5). TSSs in the close proximity of each other give rise to a functionally equivalent set of transcripts and are likely regulated by the same promoter elements. Thus, TSSs can be spatially clustered into larger transcriptional units, called tag clusters (TCs) that correspond to individual promoters. CAGEr supports three methods for spatial clustering of TSSs along the genome, two ab initio methods driven by the data itself, as well as assigning TSSs to predefined genomic regions:"

ce <- clusterCTSS(ce,
           threshold=1,
           thresholdIsTpm = TRUE,
           nrPassThreshold = 1,
           method="distclu",
           maxDist=20,
           removeSingletons = TRUE,
           keepSingletonsAbove = 3)

# Let's have a look what the result looks like
head(tagClustersGR(ce, sample = "BMDM"))

# calculate cumulative distribution for every tag cluster in each of the samples
ce <- cumulativeCTSSdistribution(ce, clusters = "tagClusters", useMulticore = T)
# determine the positions of selected quantiles
ce <- quantilePositions(ce, clusters = "tagClusters", qLow = 0.1, qUp = 0.9)

# How many tagclusters do we have in total?
length(tagClustersGR(ce, sample = "BMDM"))

# histogram of interquantile width
plotInterquantileWidth(ce, clusters = "tagClusters", tpmThreshold = 3, qLow = 0.1, qUp = 0.9)

print("Retrieving clusters as GenomicRanges")
clusters_gr <- tagClustersGR(ce, sample="BMDM")

#------------------------------------------------------------------------------------
#----------------------   Liftover coordinates to mm10
#------------------------------------------------------------------------------------

print("Liftover to mm10")

# * Now we can lift over the intervals to mm10
# * Annotate them with peakanno
# * pick the most highly expressed one for each gene

liftover <- function(peaks_gr_mm9){ #input is a GenomicRanges object in mm9 coordinates
  #lift peak locations from mm9 to mm10
  chain <- rtracklayer::import.chain(opt$liftoverchain)
  on.exit( close( file(opt$liftoverchain)) )

  peaks_gr_mm10 <- rtracklayer::liftOver(peaks_gr_mm9, chain)
  peaks_gr_mm10 <- GenomicRanges::GRanges(unlist(peaks_gr_mm10))

  return(peaks_gr_mm10)
}

mac_cage_mm10 <- liftover( clusters_gr )
ggplot(as.data.frame(mac_cage_mm10), aes(x=width)) +
  geom_histogram(bins = 100)

# Liftover coordinates of dominant_ctss
dominant_ctss <- liftover( 
  GRanges(
    seqnames = seqnames(clusters_gr), 
    ranges = IRanges(start = clusters_gr$dominant_ctss,
                     end = clusters_gr$dominant_ctss),
    score=clusters_gr$score)
)

#------------------------------------------------------------------------------------
#-----------------   Annotate TSS clusters to reference gene coordinates
#------------------------------------------------------------------------------------
print("Annotate TSS coordinates")
# Use coordinates of the dominant ctss downstream

mac_cage_anno <- ChIPseeker::annotatePeak(dominant_ctss, 
                                          tssRegion=c(-1000, 1000), #more stringent than default
                                          level = "gene",
                                          TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene, 
                                          annoDb = "org.Mm.eg.db")

# For those that are reasonably close to a TSS,
# check for each gene, which position has the highest score.

mac_cage_maxscore <- as.data.frame(mac_cage_anno) %>% 
  filter(abs(distanceToTSS)<=30000) %>%
  mutate(SYMBOL=as.factor(SYMBOL)) %>%
  filter(SYMBOL!="") %>%
  group_by(SYMBOL)%>%
  filter(score == max(score))%>%
  filter(distanceToTSS == min(distanceToTSS )) # for tied score, use shorter distance

nrow(mac_cage_maxscore)

ggplot(mac_cage_maxscore, aes(x=distanceToTSS)) +
  geom_histogram(bins = 100)

#------------------------------------------------------------------------------------
#----------   retrieve gene coordinates and promoterregion from reference
#------------------------------------------------------------------------------------

gencode_mm10_geneanno <- rtracklayer::import.gff(opt$gencode_mm10_geneanno)
genecoords <- as.data.frame(gencode_mm10_geneanno) %>%
  dplyr::select("seqnames","start","end","strand","gene_id") %>%
  mutate(score=0) %>%
  dplyr::mutate(gene_id=gsub("\\.[0-9]*$","",gene_id)) %>%
  dplyr::filter(!"gene_id"=="")%>%
  dplyr::select("seqnames","start","end","gene_id","score","strand")


gencode_mm10_promoterregions <- promoters(gencode_mm10_geneanno)
gencode_mm10_promoterregions <- as.data.frame(gencode_mm10_promoterregions) %>%
  dplyr::select("seqnames","start","end","strand","gene_id") %>%
  mutate(score=0) %>%
  dplyr::mutate(gene_id=gsub("\\.[0-9]*$","",gene_id)) %>%
  dplyr::filter(!"gene_id"=="")%>%
  dplyr::select("seqnames","start","end","gene_id","score","strand")

#------------------------------------------------------------------------------------
#-----------------   export files
#------------------------------------------------------------------------------------

write.table(genecoords, 
            file = paste0(outdir,"reference_genecoords.bed"),
            sep="\t", 
            col.names = FALSE,
            quote=FALSE,
            row.names = FALSE)

write.table(gencode_mm10_promoterregions, 
            file = paste0(outdir,"reference_promoterregions.bed"),
            sep="\t", 
            col.names = FALSE,
            quote=FALSE,
            row.names = FALSE)


rtracklayer::export.bed(as.data.frame(mac_cage_mm10),
                        paste0(outdir, "mac_cage_tssclusterregions.bed"),
                        format="bed")

rtracklayer::export.bed(as.data.frame(dominant_ctss),
                        paste0(outdir, "mac_cage_dominant_ctss.bed"),
                        format="bed")

rtracklayer::export.bed(mac_cage_maxscore,
                        paste0(outdir,"mac_cage_maxscore.bed"),
                        format="bed")