Workflow Steps and Code Snippets

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( "--summit_granges"),
              type="character",
              help="Path to rds file of summits in granges format with directionschange and distancetoTSS as additional metadata columns"),
  make_option(c("--memedb_expressed"),
              type="character",
              help="Path to memedb file filtered for motifs where TFs are expressed in 4sU"))

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

suppressPackageStartupMessages(library(here, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ggplot2, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(memes, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(universalmotif, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(dplyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(plyranges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(GenomicRanges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm10, warn.conflicts=F, quietly=T))

#-------------------------------
## Import reference for sequence
#-------------------------------
mm.genome <- BSgenome.Mmusculus.UCSC.mm10

#-------------------------------
## read in prepared data
#-------------------------------

ChIPseq_summit_Granges <- readRDS(opt$summit_granges)

# Take 100bp windows around ChIP-seq summits
summit_flank_100bp <- ChIPseq_summit_Granges %>% 
  plyranges::anchor_center() %>% 
  plyranges::mutate(width = 100) 

# Take 100bp windows around ChIP-seq summits
summit_flank_1000bp <- ChIPseq_summit_Granges %>% 
  plyranges::anchor_center() %>% 
  plyranges::mutate(width = 1000) 

meme_db_expressed <- readRDS(opt$memedb_expressed)


# to_list() converts the database back from data.frame format to a standard `universalmotif` object.
options(meme_db = to_list(meme_db_expressed, extrainfo = FALSE))

# where is meme installed 
my_memepath="~/software/meme/bin/"
check_meme_install(meme_path=my_memepath)

summit_flank_100bp_seq <- summit_flank_100bp %>%
  get_sequence(mm.genome)

summit_flank_1000bp_seq <- summit_flank_1000bp %>%
  get_sequence(mm.genome)


#-------------------------------
## run fimo
#-------------------------------

fimo_results <-
  runFimo(summit_flank_1000bp_seq,
          meme_db_expressed,
          meme_path=my_memepath)

saveRDS (fimo_results, here("results/current/memes_bioc/fimo_1000bp/fimo.rds") )

print("Finished running fimo")
print("Analysis DONE")

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

option_list <- list(
  make_option(c("--ABC_all"), 
              type="character",
              help="path to abc results of dexlps condition"),
  make_option(c("--memedb_expressed"),
              type="character",
              help="Path to memedb file filtered for motifs where TFs are expressed in 4sU"),
  make_option(c("--output"),
              type="character",
              help="fimo results file")
  )

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

suppressPackageStartupMessages(library(here, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(ggplot2, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(memes, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(universalmotif, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(dplyr, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(plyranges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(GenomicRanges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm10, warn.conflicts=F, quietly=T))

#-------------------------------
## Import reference for sequence
#-------------------------------
mm.genome <- BSgenome.Mmusculus.UCSC.mm10

#-------------------------------
## read in prepared data
#-------------------------------
ABC_all <- read.delim(opt$ABC_all) %>% plyranges::as_granges(., seqnames=chr)

# no need to run fimo a bunch of times on the same enhancers regions, just because they are listed more than once (with different ABCscores)
ABC_unique <- unique(ABC_all)

meme_db_expressed <- readRDS(opt$memedb_expressed)

# to_list() converts the database back from data.frame format to a standard `universalmotif` object.
options(meme_db = to_list(meme_db_expressed, extrainfo = FALSE))

# where is meme installed 
my_memepath="~/software/meme/bin/"
check_meme_install(meme_path=my_memepath)

#-------------------------------
## get sequences
#-------------------------------

enhancer_seq <- ABC_unique %>%
  get_sequence(mm.genome)

#-------------------------------
## run fimo
#-------------------------------

# conda activate py_3
# perlbrew use perl-5.34.0
# nohup Rscript memes_runanalyses_ABCenhancerregions.r & (from within the script directory)

fimo_results <-
  runFimo(enhancer_seq,
          meme_db_expressed,
          meme_path=my_memepath)

print("Finished running fimo")

saveRDS (fimo_results, opt$output)

print("Done saving results")

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

option_list <- list(
  make_option(c( "--summit_granges"),
              type="character",
              help="Path to rds file of summits in granges format with directionschange and distancetoTSS as additional metadata columns"),
  make_option(c("--memedb_expressed"),
              type="character",
              help="Path to memedb file filtered for motifs where TFs are expressed in 4sU"))

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

suppressPackageStartupMessages(library(here, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(memes, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(universalmotif, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(GenomicRanges, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm10, warn.conflicts=F, quietly=T))
suppressPackageStartupMessages(library(plyranges, warn.conflicts=F, quietly=T))

#-------------------------------
## Import reference for sequence
#-------------------------------
mm.genome <- BSgenome.Mmusculus.UCSC.mm10

#-------------------------------
## read in prepared data
#-------------------------------

ChIPseq_summit_Granges <- readRDS(opt$summit_granges)

# Take 100bp windows around ChIP-seq summits
summit_flank_100bp <- ChIPseq_summit_Granges %>% 
  plyranges::anchor_center() %>% 
  plyranges::mutate(width = 100) 

# Take 100bp windows around ChIP-seq summits
summit_flank_1000bp <- ChIPseq_summit_Granges %>% 
  plyranges::anchor_center() %>% 
  plyranges::mutate(width = 1000) 

meme_db_expressed <- readRDS(opt$memedb_expressed)

# to_list() converts the database back from data.frame format to a standard `universalmotif` object.
options(meme_db = to_list(meme_db_expressed, extrainfo = FALSE))

# where is meme installed 
my_memepath="~/software/meme/bin/"
check_meme_install(meme_path=my_memepath)

#-------------------------------
# define inputs
#-------------------------------

summit_flank_100bp_seq <- summit_flank_100bp %>%
  get_sequence(mm.genome)

summit_flank_1000bp_seq <- summit_flank_1000bp %>%
  get_sequence(mm.genome)

summit_flank_seq_bydirchange <- summit_flank_100bp %>%
  # remove unchanged ones and only compare "up" vs "down"
  filter(directionchange !="ns")%>%
  # Get a list of chip peaks belonging to each set
  split(mcols(.)$directionchange) %>%
  # look up the DNA sequence of each peak within each group
  get_sequence(mm.genome)

#-------------------------------
## up vs downregulation
# run by directionchange to discover consensus motif separately
#-------------------------------

#-------------------------------
# STREME
#-------------------------------

print("Start running streme for 100bp")

stremeout_100bp_down <-  here("results/current/memes_bioc/streme_100bp_down/streme.xml")
if (!file.exists( stremeout_100bp_down )){
  runStreme(summit_flank_seq_bydirchange[["down"]], control="shuffle", objfun="de",
            meme_path="~/software/meme/bin/", silent=FALSE,
            outdir = dirname(stremeout_100bp_down))
}

stremeout_100bp_up <- here("results/current/memes_bioc/streme_100bp_up/streme.xml")
if (!file.exists( stremeout_100bp_up )){
  runStreme(summit_flank_seq_bydirchange[["up"]], control="shuffle", objfun="de",
            meme_path="~/software/meme/bin/",
            outdir = dirname(stremeout_100bp_up))
}

print("Finished running streme for up- and down-regions")


#-------------------------------
# DREME
#-------------------------------

dremeout_100bp_down <- here("results/current/memes_bioc/dreme_100bp_down/dreme.xml")
if (!file.exists(dremeout_100bp_down)){
  runDreme(summit_flank_seq_bydirchange[["down"]], "shuffle",
           meme_path="~/software/meme/bin/",
           outdir = dirname(dremeout_100bp_down))
  }

dremeout_100bp_up <- here("results/current/memes_bioc/dreme_100bp_up/dreme.xml")
if (!file.exists(dremeout_100bp_up)){
  runDreme(summit_flank_seq_bydirchange[["up"]], "shuffle",
           meme_path="~/software/meme/bin/",
           outdir = dirname(dremeout_100bp_up))
  }

print("Done running DREME")

#-------------------------------
## run ame - discriminative mode
#-------------------------------

# enriched in upregulated with "down" as control
ame_discr_up <- here("results/current/memes_bioc/ame_discr_up/ame.tsv")
if (!file.exists( ame_discr_up )){
  runAme(summit_flank_seq_bydirchange, control = "down",
         meme_path=my_memepath,
         outdir=dirname(ame_discr_up))
}
# enriched in downregulated with "up" as control
ame_discr_down <- here("results/current/memes_bioc/ame_discr_down/ame.tsv")
if (!file.exists(ame_discr_down)){
  runAme(summit_flank_seq_bydirchange, control = "up",
         meme_path=my_memepath,
         outdir=dirname(ame_discr_down))
}

print("Finished running ame in discriminative mode for direction of expressionchange")


#-------------------------------
## all summits
#-------------------------------

## run streme to discover consensus motif
#-------------------------------
#print("Starting streme 1000bp")
#stremeout_1000bp <- here("results/current/memes_bioc/streme_1000bp/streme.xml")
#if (!file.exists( stremeout_1000bp )){
#  runStreme(summit_flank_1000bp_seq, control="shuffle",
#            meme_path="~/software/meme/bin/",
#            outdir = dirname(stremeout_1000bp) )
#}

print("Starting streme 100bp for all summits")
stremeout_100bp <- here("results/current/memes_bioc/streme_100bp/streme.xml")
if (!file.exists( stremeout_100bp )){
  runStreme(summit_flank_100bp_seq, control="shuffle",
            meme_path="~/software/meme/bin/",
            outdir = dirname(stremeout_100bp) )
}
print("Finished running streme for 100bp summit regions")

# Option objfun="cd" does not seem to get passed on to streme
#print("Starting streme 1000bp with central enrichment")
#stremeout_cd_1000bp <- here("results/current/memes_bioc/streme_cd_1000bp/streme.xml")
#if (!file.exists( stremeout_cd_1000bp )){
#  runStreme(summit_flank_1000bp_seq[1:100], objfun="cd", control=NA,
#            meme_path="~/software/meme/bin/",
#            outdir = dirname(stremeout_cd_1000bp) )
#}
#print("Finished running streme for 1000bp summit regions")

print("Analysis DONE")

if (interactive()) {
  library(methods)
  Snakemake <- setClass(
    "Snakemake", 
    slots=c(
      input='list', 
      output='list',
      wildcards='list',
      threads='numeric'
    )
  )
  snakemake <- Snakemake(
      input=list(gtf="gencode.vM21.annotation.mRNA_ends_found.gtf.gz"),
      output=list(gtf="/fscratch/fanslerm/gencode.vM21.annotation.mRNA_ends_found.txcutr.w500.gtf",
                  fa="/fscratch/fanslerm/gencode.vM21.annotation.mRNA_ends_found.txcutr.w500.fa"),
      wildcards=list(width="500"),
      threads=1
  )
}

################################################################################
## Libraries and Parameters
################################################################################

library(txcutr)
library(BSgenome.Mmusculus.UCSC.mm10)
library(GenomicFeatures)

mm10 <- BSgenome.Mmusculus.UCSC.mm10

maxTxLength <- as.integer(snakemake@wildcards$width)

## set cores
BiocParallel::register(BiocParallel::MulticoreParam(snakemake@threads))

################################################################################
## Load Data, Truncate, and Export
################################################################################

txdb <- makeTxDbFromGFF(file=snakemake@input$gtf,
                        organism="Mus musculus",
                        taxonomyId=10090,
                        chrominfo=seqinfo(mm10))

txdb_result <- truncateTxome(txdb, maxTxLength)

exportGTF(txdb_result, snakemake@output$gtf)

exportFASTA(txdb_result, mm10, snakemake@output$fa)

if (interactive()) {
  library(methods)
  Snakemake <- setClass(
    "Snakemake", 
    slots=c(
      input='list', 
      output='list', 
      wildcards='list',
      threads='numeric'
    )
  )
  snakemake <- Snakemake(
    input=list(),
    output=list(gtf="txdb.mm10.ensGene.txcutr.w{width}.gtf",
                fa="txdb.mm10.ensGene.txcutr.w{width}.fa"),
    wildcards=list(width="500"),
    threads=1
  )
}

################################################################################
## Libraries and Parameters
################################################################################

library(txcutr)
library(BSgenome.Mmusculus.UCSC.mm10)
library(TxDb.Mmusculus.UCSC.mm10.ensGene)

mm10 <- BSgenome.Mmusculus.UCSC.mm10
txdb <- TxDb.Mmusculus.UCSC.mm10.ensGene

maxTxLength <- as.integer(snakemake@wildcards$width)

## set cores
BiocParallel::register(BiocParallel::MulticoreParam(snakemake@threads))

################################################################################
## Truncate and Export
################################################################################

txdb_result <- truncateTxome(txdb, maxTxLength)

exportGTF(txdb_result, snakemake@output$gtf)

exportFASTA(txdb_result, mm10, snakemake@output$fa)

myargs <- commandArgs(trailingOnly=TRUE)
bamfile <- myargs[1]
species <- myargs[2]

print("loading packages (ATACseqQC, ggplot, etc)...")
suppressPackageStartupMessages(library(ggplot2, quietly=TRUE))
suppressPackageStartupMessages(library(Rsamtools, quietly=TRUE))
suppressPackageStartupMessages(library(ATACseqQC, quietly=TRUE))
suppressPackageStartupMessages(library(ChIPpeakAnno, quietly=TRUE))
suppressPackageStartupMessages(library("GenomicAlignments", quietly=TRUE))

if (species == "mm") {
  suppressPackageStartupMessages(library(TxDb.Mmusculus.UCSC.mm10.knownGene, quietly=TRUE))
  suppressPackageStartupMessages(library(BSgenome.Mmusculus.UCSC.mm10, quietly=TRUE))
  txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
  bsgenome <- BSgenome.Mmusculus.UCSC.mm10
  genome <- Mmusculus
  print("species is 'mm' using mm10 for analysis")
  ### Note: Everything below is deprecated until I can figure out a way to port a 
  ### static/local package with snakemake
  # Note: phastCons60way was manually curated from GenomicAlignments, built, and installed as an R package
  # score was obtained according to: https://support.bioconductor.org/p/96226/
  # package was built and installed according to: https://www.bioconductor.org/packages/devel/bioc/vignettes/GenomicScores/inst/doc/GenomicScores.html
  # (section 5.1: Building an annotation package from a GScores object)
  #suppressWarnings(suppressPackageStartupMessages(library(GenomicScores, lib.loc="/users/dia6sx/snakeATAC/scripts/", quietly=TRUE)))
  #suppressWarnings(suppressPackageStartupMessages(library(phastCons60way.UCSC.mm10, lib.loc="/users/dia6sx/snakeATAC/scripts/", quietly=TRUE)))
} else if (species == "hs") {
  suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene, quietly=TRUE))
  suppressPackageStartupMessages(library(BSgenome.Hsapiens.UCSC.hg38, quietly=TRUE))
  txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  bsgenome <- BSgenome.Hsapiens.UCSC.hg38
  genome <- Hsapiens
  print("species is 'hs' using hg38 for analysis")
} else {
  print(paste("params ERROR: ATACseqQC is not configured to use species =", species))
  print("exiting...")
  quit(status=1)
}

doATACseqQC <- function(bamfile, txdb, bsgenome, genome) {
    # Fragment size distribution
    print(paste("generating output for ",strsplit(basename(bamfile),split='\\.')[[1]][1],"...",sep=""))
    print("calculating Fragment size distribution...")
    bamfile.labels <- gsub(".bam", "", basename(bamfile))
    loc_to_save_figures <- paste(dirname(dirname(bamfile)),"/qc/ATACseqQC",sep="")
    if (file.exists(loc_to_save_figures)) {
        print("Warning: old figures will be overwritten")
    } else {
        dir.create(loc_to_save_figures)
    }
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_fragment_size_distribution.png",sep="")
    png(png_file)
    fragSizeDist(bamfile, bamfile.labels)
    dev.off()

    # Adjust the read start sites
    print("adjusting read start sites...")
    ## bamfile tags to be read in
    possibleTag <- list("integer"=c("AM", "AS", "CM", "CP", "FI", "H0", "H1", "H2", 
                                    "HI", "IH", "MQ", "NH", "NM", "OP", "PQ", "SM",
                                    "TC", "UQ"), 
                    "character"=c("BC", "BQ", "BZ", "CB", "CC", "CO", "CQ", "CR",
                                "CS", "CT", "CY", "E2", "FS", "LB", "MC", "MD",
                                "MI", "OA", "OC", "OQ", "OX", "PG", "PT", "PU",
                                "Q2", "QT", "QX", "R2", "RG", "RX", "SA", "TS",
                                "U2"))
    bamTop100 <- scanBam(BamFile(bamfile, yieldSize = 100),
                     param = ScanBamParam(tag=unlist(possibleTag)))[[1]]$tag
    tags <- names(bamTop100)[lengths(bamTop100)>0]
    ## files will be output into outPath
    ## shift the coordinates of 5'ends of alignments in the bam file
    outPath <- paste(dirname(dirname(bamfile)),"/alignments_shifted", sep="")
    seqinformation <- seqinfo(txdb)
    gal <- readBamFile(bamfile, tag=tags, asMates=TRUE, bigFile=TRUE)
    shiftedBamfile <- file.path(outPath, paste(bamfile.labels,"_shifted.bam",sep=""))
    # check if shifted Bam file exists from previous run
    if (file.exists(shiftedBamfile)) {
        print("Shifted Bamfile found.")
        print("Loading in...")
        gal <- readBamFile(shiftedBamfile, tag=tags, asMates=TRUE, bigFile=TRUE)
        ## This step is mostly for formating so splitBam can
        ## take in bamfile. Implementing shift of 0 bp on positive strand
        ## and 0 bp on negative strand because shifted Bamfile should
        ## already have these shifts
        gal1 <- shiftGAlignmentsList(gal, positive = 0L, negative = 0L)
    } else {
        # shifted bam file does not exist check if
        # old shifted alignments directory exists
        # if so remove and create new one
        if (file.exists(outPath)){
            unlink(outPath,recursive=TRUE)
        }
        dir.create(outPath)
        print("*** creating shifted bam file ***")
        gal1 <- shiftGAlignmentsList(gal, outbam=shiftedBamfile)
    }

    # Promoter/Transcript body (PT) score
    print("calculating Promoter/Transcript body (PT) score...")
    txs <- transcripts(txdb)
    pt <- PTscore(gal1, txs)
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_ptscore.png",sep="")
    png(png_file)
    plot(pt$log2meanCoverage, pt$PT_score, 
        xlab="log2 mean coverage",
        ylab="Promoter vs Transcript",
        main=paste(bamfile.labels,"PT score"))
    dev.off()

    # Nucleosome Free Regions (NFR) score
    print("calculating Nucleosome Free Regions (NFR) score")
    nfr <- NFRscore(gal1, txs)
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_nfrscore.png",sep="")
    png(png_file)
    plot(nfr$log2meanCoverage, nfr$NFR_score, 
        xlab="log2 mean coverage",
        ylab="Nucleosome Free Regions score",
        main=paste(bamfile.labels,"\n","NFRscore for 200bp flanking TSSs",sep=""),
        xlim=c(-10, 0), ylim=c(-5, 5))
    dev.off()

    # Transcription Start Site (TSS) Enrichment Score
    print("calculating Transcription Start Site (TSS) Enrichment score")
    tsse <- TSSEscore(gal1, txs)
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_tss_enrichment_score.png",sep="")
    png(png_file)
    plot(100*(-9:10-.5), tsse$values, type="b", 
        xlab="distance to TSS",
        ylab="aggregate TSS score",
        main=paste(bamfile.labels,"\n","TSS Enrichment score",sep=""))
    dev.off()

    # Split reads, Heatmap and coverage curve for nucleosome positions
    print("splitting reads by fragment length...")
    genome <- genome
    outPath <- paste(dirname(dirname(bamfile)),"/alignments_split", sep="")
    TSS <- promoters(txs, upstream=0, downstream=1)
    TSS <- unique(TSS)
    ## estimate the library size for normalization
    librarySize <- estLibSize(bamfile)
    ## calculate the signals around TSSs.
    NTILE <- 101
    dws <- ups <- 1010
    splitBamfiles <- paste(outPath,"/",c("NucleosomeFree", 
                                             "mononucleosome",
                                             "dinucleosome",
                                             "trinucleosome"),".bam",sep="")
    # check if split Bam files exists from previous run
    if (all(file.exists(splitBamfiles))) {
        print("*** split bam files found! ***")
        print("Loading in...")
        sigs <- enrichedFragments(bamfiles=splitBamfiles,
                                    index=splitBamfiles, 
                                    TSS=TSS,
                                    librarySize=librarySize,
                                    TSS.filter=0.5,
                                    n.tile = NTILE,
                                    upstream = ups,
                                    downstream = dws)
    } else {
        # split bam files do not exist check if
        # old split alignments directory exists
        # if so remove and create new one
        if (file.exists(outPath)){
            unlink(outPath,recursive=TRUE)
        }
        print("*** creating split bam files ***")
        dir.create(outPath)
        ## split the reads into NucleosomeFree, mononucleosome, 
        ## dinucleosome and trinucleosome.
        ## and save the binned alignments into bam files.
        objs <- splitGAlignmentsByCut(gal1, txs=txs, genome=genome, outPath = outPath)
        #objs <- splitBam(bamfile, tags=tags, outPath=outPath,
            #        txs=txs, genome=genome,
            #       conservation=phastCons60way.UCSC.mm10,
            #      seqlev=paste0("chr", c(1:19, "X", "Y")))
        sigs <- enrichedFragments(gal=objs[c("NucleosomeFree", 
                                        "mononucleosome",
                                        "dinucleosome",
                                        "trinucleosome")], 
                                    TSS=TSS,
                                    librarySize=librarySize,
                                    TSS.filter=0.5,
                                    n.tile = NTILE,
                                    upstream = ups,
                                    downstream = dws)
    }
    ## log2 transformed signals
    sigs.log2 <- lapply(sigs, function(.ele) log2(.ele+1))
    ## plot heatmap
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_nucleosome_pos_heatmap.png",sep="")
    png(png_file)
    featureAlignedHeatmap(sigs.log2, reCenterPeaks(TSS, width=ups+dws),
                        zeroAt=.5, n.tile=NTILE)
    dev.off()
    ## get signals normalized for nucleosome-free and nucleosome-bound regions.
    out <- featureAlignedDistribution(sigs, 
                                    reCenterPeaks(TSS, width=ups+dws),
                                    zeroAt=.5, n.tile=NTILE, type="l", 
                                    ylab="Averaged coverage")
    ## rescale the nucleosome-free and nucleosome signals to 0~1
    range01 <- function(x){(x-min(x))/(max(x)-min(x))}
    out <- apply(out, 2, range01)
    png_file <- paste(loc_to_save_figures,"/",bamfile.labels,"_TSS_signal_distribution.png",sep="")
    png(png_file)
    matplot(out, type="l", xaxt="n", 
            xlab="Position (bp)", 
            ylab="Fraction of signal",
            main=paste(bamfile.labels,"\n","TSS signal distribution",sep=""))
    axis(1, at=seq(0, 100, by=10)+1, 
        labels=c("-1K", seq(-800, 800, by=200), "1K"), las=2)
    abline(v=seq(0, 100, by=10)+1, lty=2, col="gray")
    dev.off()

    print("QC Finished.")
    print("Generated QC figures can be found in qc folder under ATACseQC")
    print(paste("*** removing temp files in",outPath,"***"))
    unlink(outPath,recursive=TRUE)
    outPath <- paste(dirname(dirname(bamfile)),"/alignments_shifted", sep="")
    print(paste("*** removing temp files in",outPath,"***"))
    unlink(outPath,recursive=TRUE)
}

doATACseqQC(bamfile, txdb, bsgenome, genome)