shell:
  "TagBamWithReadSequenceExtended "
  "INPUT={input} "
  "OUTPUT={output.bam} "
  "SUMMARY={output.summary} "
  "BASE_QUALITY={params.base_qual} "
  "NUM_BASES_BELOW_QUALITY={params.bases_below_qual} "
  "BASE_RANGE={params.base_range} "
  "BARCODED_READ=1 "
  "DISCARD_READ=false "
  "TAG_NAME=XC "
  "2> {log}"

shell:
  "TagBamWithReadSequenceExtended "
  "INPUT={input} "
  "OUTPUT={output.bam} "
  "SUMMARY={output.summary} "
  "BASE_QUALITY={params.base_qual} "
  "NUM_BASES_BELOW_QUALITY={params.bases_below_qual} "
  "BASE_RANGE={params.base_range} "
  "BARCODED_READ=1 "
  "DISCARD_READ=true "
  "TAG_NAME=XM "
  "2> {log}"

shell:
  "FilterBAM "
  "INPUT={input} "
  "OUTPUT={output} "
  "TAG_REJECT=XQ "
  "2> {log}"

shell:
  "TrimStartingSequence "
  "INPUT={input} "
  "OUTPUT={output.bam} "
  "OUTPUT_SUMMARY={output.summary} "
  "SEQUENCE={params.adapter_sequence} "
  "MISMATCHES={params.mismatches} "
  "NUM_BASES={params.num_bases} "
  "2> {log}"

shell:
  "PolyATrimmer "
  "INPUT={input} "
  "OUTPUT={output.bam} "
  "OUTPUT_SUMMARY={output.summary} "
  "MISMATCHES={params.mismatches} "
  "NUM_BASES={params.num_bases} "
  "2> {log}"

shell:
  "picard SamToFastq "
  "INPUT={input} "
  "FASTQ={output} "
  "2> {log}"

shell:
  "picard SortSam "
  "INPUT={input} "
  "OUTPUT={output} "
  "SORT_ORDER=queryname "
  "2> {log}"

shell:
  "picard MergeBamAlignment "
  "ALIGNED_BAM={input.aligned} "
  "UNMAPPED_BAM={input.unaligned} "
  "REFERENCE_SEQUENCE={input.reference} "
  "OUTPUT={output} "
  "INCLUDE_SECONDARY_ALIGNMENTS=false "
  "PAIRED_RUN=false "
  "2> {log}"

shell:
  "TagReadWithGeneExon "
  "INPUT={input.bam} "
  "OUTPUT={output.bam} "
  "ANNOTATIONS_FILE={input.annot} "
  "TAG=GE "
  "CREATE_INDEX=true "
  "2> {log}"

shell:
  "BAMTagHistogram "
  "INPUT={input} "
  "OUTPUT={output} "
  "TAG=XC "
  "READ_QUALITY={params.read_quality} "
  "2> {log}"

script:
  "../scripts/processing/align_report.Rmd"

shell:
  "python scripts/processing/extract_dge.py -i {input.umi_obs} {params.whitelist_arg}"
  "-o {output.dge} --sample {params.reads} --seed {params.seed} "
  "--tpt_threshold {params.tpt_threshold} 2> {log}"

script:
  "../scripts/analyses/advanced_dge_downsampling.R"

script:
  "../scripts/analyses/reads_on_target.R"

script:
  "../scripts/analyses/tapseq_vs_cropseq.Rmd"

script:
  "../scripts/analyses/filter_mmix_umi_obs.R"

shell:
  "python scripts/processing/extract_dge.py -i {input.umi_obs} -o {output.dge} "
  "{params.whitelist_arg} --sample {params.reads} --seed {params.seed} --tpt_threshold "
  "{params.tpt_threshold} 2> {log}"

script:
  "../scripts/analyses/screen_data_qc.Rmd"

shell:
  "python scripts/analyses/collapse_perturbations.py -i {input.pert_status} -t {input.grna_targets} "
  "-o {output} 2> {log}"

script:
  "../scripts/analyses/differential_expression.R"

script:
  "../scripts/analyses/compare_covariates.Rmd"

script:
  "../scripts/analyses/process_de_results.R"

script:
  "../scripts/analyses/map_enhancers.Rmd"

shell:
  "wget -O {output} {params.url}"

shell:
  "wget -O {output.bam} {params.url} ; picard BuildBamIndex I={output.bam}"

shell:
  "rm -r data/k562_chromatin_data/HiC/5kb_resolution_intrachromosomal; "
  "wget -c {params.url} -O - | tar -xz K562/5kb_resolution_intrachromosomal/chr*/MAPQG0/; "
  "mv K562/5kb_resolution_intrachromosomal "
  "data/k562_chromatin_data/HiC/5kb_resolution_intrachromosomal; "
  "rmdir K562"

script:
  "../scripts/analyses/chromatin_analyses_screen.Rmd"

script:
  "../scripts/analyses/hic_analysis.Rmd"

script:
  "../scripts/analyses/abc_analysis_screen.Rmd"

shell:
  "Rscript scripts/analyses/extract_cropseq_vectors.R -b {input.bam} -d {input.dge_summary} "
  "-p {params.vector_pattern} -o {output}"

script:
  "../scripts/analyses/screen_grna_synthesis_errors.Rmd"

script:
  "../scripts/chromatin_annotated_etps/enhancer_screen_etps.R"

script:
  "../scripts/chromatin_annotated_etps/genome_wide_etps.R"

script:
  "../scripts/chromatin_annotated_etps/gasperini_screen_etps.R"

script:
  "../scripts/chromatin_annotated_etps/gasperini_pilot_etps.R"

script:
  "../scripts/chromatin_annotated_etps/fulco_tiling_etps.R"

shell:
  "wget -O {output.gtf}.gz {params.gtf_url} 2> {log} && gzip -d {output.gtf}.gz ; "
  "wget -O {output.fasta}.gz {params.fasta_url} 2>> {log} && gzip -d {output.fasta}.gz"

shell:
  "Rscript scripts/processing/create_tapseq_annot.R -t {input} -b {params.BSgenome} "
  "-f {output.fasta} -o {output.gtf} 2> {log}"

shell:
  "python scripts/processing/cropseq_vector_reference.py -i {input} -o {params.output_bn} "
  "--prefix {params.vector_prefix}"

shell:
  "touch {output.fasta} {output.gtf}"

shell:
  "cat {input.genome_gtf} > {output.cropseq_ref_gtf} && "
  "cat {input.vectors_gtf} >> {output.cropseq_ref_gtf} ; "
  "cat {input.genome_fasta} > {output.cropseq_ref_fasta} && "
  "cat {input.vectors_fasta} >> {output.cropseq_ref_fasta} ; "

shell:
  "picard CreateSequenceDictionary "
  "REFERENCE={input} "
  "OUTPUT={output} "
  "2> {log}"

shell:
  "ConvertToRefFlat "
  "ANNOTATIONS_FILE={input.annot} "
  "SEQUENCE_DICTIONARY={input.seq_dict} "
  "OUTPUT={output} "
  "2> {log}"

shell:
  "mkdir -p {output} && "  # create genomeDir directory
  "STAR --runThreadN {threads} "
  "--runMode genomeGenerate "
  "--genomeDir {output} "
  "--genomeFastaFiles {input.fasta} "
  "--sjdbGTFfile {input.gtf} "
  "--sjdbOverhang {params.sjdb_overhang} "
  "--outFileNamePrefix {params.outprefix} ; "
  "mv data/alignment_references/{wildcards.align_ref}/star.Log.out {log}"

shell:
  "cmd='fastq-dump --split-files --gzip --outdir {params.outdir} {wildcards.run}' ; "
  "echo $cmd > {log} ; $cmd >> {log}"

shell:
  "mv -v {params.outdir}/{params.access}_1.fastq.gz {output.fwd} > {log} & "
  "mv -v {params.outdir}/{params.access}_2.fastq.gz {output.rev} >> {log}"

shell:
  "bash scripts/download_fastq/merge_fastq.sh {params.outdir} {output.fwd} {output.rev} "
  "> {log}"

shell:
  "bash scripts/download_fastq/merge_fastq_i7.sh {params.outdir} {output.fwd} {output.rev} "
  "> {log}"

shell:
  "GatherMolecularBarcodeDistributionByGene "
  "INPUT={input} "
  "OUTPUT={output} "
  "NUM_CORE_BARCODES={params.ncells} "
  "EDIT_DISTANCE={params.edit_distance} "
  "READ_MQ={params.read_mq} "
  "MIN_BC_READ_THRESHOLD={params.min_umi_reads} "
  "RARE_UMI_FILTER_THRESHOLD={params.rare_umi_filter} "
  "2> {log}"

shell:
  "python scripts/processing/extract_dge.py -i {input.umi_obs} -o {output.dge} "
  "{params.whitelist_arg} --tpt_threshold {params.tpt_threshold} 2> {log}"

shell:
  "Rscript scripts/processing/perturbation_status.R --infile {input} --outfile {output} "
  "--vector_patter {params.vector_prefix} --min_txs {params.min_txs} --trim 2> {log}"

script:
  "../scripts/processing/dge_report.Rmd"

script:
  "../scripts/processing/dge_report_no_perts.Rmd"

/* Define a margin before every image element */
img {
  margin-top: 3em;
}

# load chromatin annotated ETPs generated by scripts/chromatin_annotated_etps/enhancer_screen_etps.R
pairs <- read.csv(here(snakemake@input), stringsAsFactors = FALSE)

# plot abc score vs distance to TSS
ggplot(pairs, aes(x = dist_to_tss, y = abc_score_norm, color = prop_grna_hits)) +
  facet_wrap(~enh_chr, ncol = 1) +
  geom_point(data = filter(pairs, significant == 0), color = "gray69", alpha = 0.25) +
  geom_point(data = filter(pairs, significant == 1)) +
  labs(x = "Distance to TSS", y = "ABC score", "ABC score vs distance to TSS") +
  scale_color_gradient(low = "blue", high = "orangered", na.value = "gray69") +
  scale_x_log10() +
  scale_y_log10() +
  theme_bw()

# plot abc score vs pvalue
ggplot(pairs, aes(x = abc_score_norm, y = -log10(pvalue), color = as.factor(significant))) +
  facet_wrap(~enh_chr, ncol = 2) +
  geom_point() +
  labs(x = "ABC score", color = "Significant") +
  scale_color_manual(values = c("1" = "steelblue", "0" = "gray69")) +
  theme_bw()

# function to calculate area under roc curve
calc_auroc <- function(category, prediction) {

  # order by prediction
  dat <- data.frame(category, prediction) %>%
    arrange(desc(prediction))

  # compute auc
  roc_obj <- roc(dat$category, dat$prediction)
  as.numeric(auc(roc_obj))

}

# only retain data on pairs that have at least one significant hit
sig_pairs <- pairs %>%
  group_by(sample, gene) %>%
  filter(sum(significant) > 0) %>%
  select(-c(2, 5:14, 23:29, 32))

# convert to long format and set pretty predictor names
sig_pairs_long <- sig_pairs %>%
  gather(key = "predictor", value = "predictor_value",
         -c(sample, perturbation, gene, significant)) %>%
  mutate(predictor = case_when(
    predictor == "int_freq" ~ "Hi-C int. freq.",
    predictor == "dist_to_tss" ~ "Distance to TSS",
    predictor == "activity_score" ~ "Activity score",
    predictor == "abc_score_norm" ~ "ABC score",
    TRUE ~ predictor))

# calculate auroc for all genes with significant enhancers
auroc <- sig_pairs_long %>%
  group_by(sample, predictor) %>%
  summarize(auc = calc_auroc(category = significant, prediction = predictor_value))

# plot auroc per gene and predictor
ggplot(auroc, aes(x = predictor, y = auc)) +
  facet_wrap(~sample, ncol = 2) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(y = "AUC", title = "AUROC") +
  coord_flip() +
  theme_bw() +
  theme(axis.title.x = element_blank())

# calculate mean and median auroc
auroc_stats <- auroc %>%
  group_by(predictor) %>%
  summarize(avg_auroc = mean(auc),
            med_auroc = median(auc))

# add predictor type for plot
auroc_stats <- auroc_stats %>%
  mutate(type = case_when(
    predictor %in% c("H3K27ac", "H3K27me3", "H3K4me1", "H3K4me3") ~ "Histone PTM",
    predictor %in% c("Hi-C int. freq.", "Distance to TSS") ~ "Proximity",
    predictor %in% c("Activity score", "ABC score") ~ "ABC model",
    predictor == "Dnase.seq" ~ "DNA accessibility",
    predictor == "POLR2A" ~ "Transcription"))

# plot average auroc per predictor
ggplot(auroc_stats, aes(x = fct_reorder(predictor, desc(avg_auroc)), y = avg_auroc, fill = type)) +
  geom_bar(stat = "identity") +
  labs(x = "Predictor", y = "Average AUROC", fill = "Type") +
  scale_fill_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank())

# only retain data on pairs that have at least one significant hit within 5 Mb of the TSS
sig_pairs_1mb <- pairs %>%
  group_by(sample, gene) %>%
  filter(dist_to_tss <= 1e6) %>%
  filter(sum(significant) > 0) %>%
  select(-c(2, 5:14, 23:29, 32))

# compute number of enhancers per gene
enh_per_gene <- summarize(sig_pairs_1mb, enhancers = n())

# plot number of enhancers per gene
ggplot(enh_per_gene, aes(x = gene, y = enhancers)) +
  facet_wrap(~sample, scales = "free_x") +
  geom_bar(stat = "identity") +
  labs(x = "Gene", y = "Enhancers", title = "Enhancers per gene within 1Mb") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# convert to long format and set pretty predictor names
sig_pairs_1mb_long <- sig_pairs_1mb %>%
  gather(key = "predictor", value = "predictor_value",
         -c(sample, perturbation, gene, significant)) %>%
  mutate(predictor = case_when(
    predictor == "int_freq" ~ "Hi-C int. freq.",
    predictor == "dist_to_tss" ~ "Distance to TSS",
    predictor == "activity_score" ~ "Activity score",
    predictor == "abc_score_norm" ~ "ABC score",
    TRUE ~ predictor))

# calculate auroc for all genes with significant enhancers
auroc_1mb <- sig_pairs_1mb_long %>%
  group_by(sample, predictor) %>%
  summarize(auc = calc_auroc(category = significant, prediction = predictor_value))

# plot auroc per gene and predictor
ggplot(auroc_1mb, aes(x = predictor, y = auc)) +
  facet_wrap(~sample, ncol = 2) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(y = "AUC", title = "AUROC 1Mb") +
  coord_flip() +
  theme_bw() +
  theme(axis.title.x = element_blank())

# calculate mean and median auroc
auroc_stats_1mb <- auroc_1mb %>%
  group_by(predictor) %>%
  summarize(avg_auroc = mean(auc),
            med_auroc = median(auc))

# add predictor type for plot
auroc_stats_1mb <- auroc_stats_1mb %>%
  mutate(type = case_when(
    predictor %in% c("H3K27ac", "H3K27me3", "H3K4me1", "H3K4me3") ~ "Histone PTM",
    predictor %in% c("Hi-C int. freq.", "Distance to TSS") ~ "Proximity",
    predictor %in% c("Activity score", "ABC score") ~ "ABC model",
    predictor == "Dnase.seq" ~ "DNA accessibility",
    predictor == "POLR2A" ~ "Transcription"))

# plot average auroc per predictor
ggplot(auroc_stats_1mb, aes(x = fct_reorder(predictor, desc(avg_auroc)), y = avg_auroc,
                            fill = type)) +
  geom_bar(stat = "identity") +
  labs(x = "Predictor", y = "Average AUROC", fill = "Type") +
  scale_fill_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank())

# function to calculate roc curve
calc_roc <- function(category, prediction) {

  # order by prediction
  dat <- data.frame(category, prediction) %>%
    arrange(desc(prediction))

  # compute roc
  roc_obj <- roc(dat$category, dat$prediction)

  # return tpr and fpr
  data.frame(TPR = rev(roc_obj$sensitivities), FPR = rev(1 - roc_obj$specificities))

}

# calculate roc for all genes with significant enhancers
roc_1mb <- sig_pairs_1mb_long %>%
  group_by(sample, predictor) %>%
  do(calc_roc(category = .$significant, prediction = .$predictor_value))

# order according to average auroc
order <- auroc_stats_1mb %>%
  arrange(desc(avg_auroc)) %>%
  pull(predictor)

roc_1mb <- roc_1mb %>%
  ungroup() %>%
  mutate(predictor = factor(predictor, levels = order))

# create labels to add auroc to plots
auroc_labels <- auroc_1mb %>%
  left_join(select(auroc_stats_1mb, predictor, avg_auroc), by = "predictor") %>%
  mutate(predictor = factor(predictor, levels = order)) %>%
  mutate(label = paste0("AUC: ", round(auc, digits = 2))) %>%
  mutate(x = 0.8, y = if_else(sample == "chr11", true = 0.25, false = 0.1))

# plot roc curves for every predictor
ggplot(roc_1mb, aes(x = FPR, y = TPR, color = sample)) +
  facet_wrap(~predictor, ncol = 2) +
  geom_line(size = 1) +
  geom_abline(lty = "dashed", color = "darkgray") +
  geom_text(data = auroc_labels, aes(x, y, label = label), show.legend = FALSE) +
  labs(color = "Region") +
  scale_color_manual(values = c("chr11" = "indianred3", "chr8" = "steelblue")) +
  theme_bw() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 13))

library(tidyverse)
library(here)

# define functions ---------------------------------------------------------------------------------

# calculate transcripts-per-transcript for umi observations to filter for chimeric reads
calc_tpt <- function(umi_obs) {

  # calculate total number of reads for every cbc - umi combination
  total_reads <- umi_obs %>%
    group_by(Cell.Barcode, Molecular_Barcode) %>%
    summarize(total_reads = sum(Num_Obs))

  # add total reads to umi_obs and calculate TPT
  umi_obs %>%
    left_join(total_reads, by = c("Cell.Barcode", "Molecular_Barcode")) %>%
    mutate(tpt = Num_Obs / total_reads) %>%
    select(-total_reads)

}

# filter according to a specified tpt threshold
filter_tpt <- function(umi_obs_tpt, tpt_threshold) {
  umi_obs_tpt %>%
    filter(tpt >= tpt_threshold) %>%
    select(-tpt)
}

# downsample umi_obs for a specifc number of genic reads
downsample_reads <- function(umi_obs, reads, level = c("cell", "sample")) {

  level <- match.arg(level)

  # "uncount" umi_observations so that each read is a row by repeating each line number-of-reads
  # times (i.e. num_obs)
  reps <- rep(seq_len(nrow(umi_obs)), times = umi_obs$Num_Obs)
  umi_obs_reads <- umi_obs[reps, -4]

  # randomly sample reads per cell or the whole sample, depending on level parameter
  if (level == "cell") {
    sampled_reads <- umi_obs_reads %>%
      group_by(Cell.Barcode) %>%
      sample_n(size = reads, replace = FALSE)
  }else{
    sampled_reads <- sample_n(umi_obs_reads, size = reads, replace = FALSE)
  }

  # count observations of each sample-cell-umi combination
  sampled_umi_obs <- sampled_reads %>%
    group_by(Cell.Barcode, Gene, Molecular_Barcode) %>%
    summarize(Num_Obs = n()) %>%
    ungroup()

}

# extract dge data for each sample by calculating number of umis per gene and  cell
extract_dge <- function(umi_obs){
  umi_obs %>%
    count(Cell.Barcode, Gene) %>%
    rename(txs = n)
}

# convert dge from long format to wide matrix
create_dge_matrix <- function(dge_long) {
  dge_long %>%
    spread(key = Cell.Barcode, value = txs, fill = 0) %>%
    as.data.frame(stringsAsFactors == FALSE) %>%
    rename(GENE = Gene)
}


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

# get parameters
sample_reads <- as.numeric(snakemake@wildcards$reads)
sampling <- match.arg(snakemake@wildcards$sampling, choices = c("fixed", "avg"))
panel <- match.arg(snakemake@wildcards$panel, choices = c("onGenes", "onTarget"))

# set sed for reproducible sampling
set.seed(snakemake@params$seed)

# load input data
umi_obs <- read.table(here(snakemake@input$umi_obs), header = TRUE, stringsAsFactors = FALSE,
                      sep = "\t")
target_genes <- read.csv(here(snakemake@input$target_genes), stringsAsFactors = FALSE)
cbc_whitelist <- readLines(here(snakemake@input$whitelist))

# filter for cell barcodes on CBC whitelist
umi_obs_filt <- filter(umi_obs, Cell.Barcode %in% cbc_whitelist)

# remove any vector transcripts if vector_pattern is provided
if (!is.null(snakemake@params$vector_pattern)) {

  umi_obs_filt <- filter(umi_obs_filt, !grepl(Gene, pattern = snakemake@params$vector_pattern))

}

# only retain target genes if panel is set to onTarget
if (panel == "onTarget") {

  umi_obs_filt <- filter(umi_obs_filt, Gene %in% target_genes$gene)

}


# downsample genic reads ---------------------------------------------------------------------------

if (sampling == "fixed") {

  # calculate number of reads per cell
  reads_per_cell <- umi_obs_filt %>%
    group_by(Cell.Barcode) %>%
    summarize(reads = sum(Num_Obs))

  # filter for cells that have at least the desired sampling size of reads
  cells_filt <- reads_per_cell %>%
    filter(reads >= sample_reads) %>%
    pull(Cell.Barcode)

  # retain umi observations ony for these cells
  umi_obs_filt <- filter(umi_obs_filt, Cell.Barcode %in% cells_filt)

  # downsample umi observations
  ds_umi_obs <- downsample_reads(umi_obs_filt, reads = sample_reads, level = "cell")

}else{

  # calculate number of reads to sample from umi_observation
  ncells <- n_distinct(umi_obs_filt$Cell.Barcode)
  sample_reads <- ncells * sample_reads

  # abort if scaled sampling size is larger than the number of available reads
  if (sample_reads > sum(umi_obs_filt$Num_Obs)) {
    stop("Not enough input reads for desired sampling size per sample!")
  }

  # downsample umi observations
  ds_umi_obs <- downsample_reads(umi_obs_filt, reads = sample_reads, level = "sample")

}

# extract dge --------------------------------------------------------------------------------------

# calculate tpt and filter for chimeric reads using the specified tpt threshold
umi_obs_tpt <- calc_tpt(ds_umi_obs)
umi_obs_tpt <- filter_tpt(umi_obs_tpt, tpt_threshold = snakemake@params$tpt_threshold)

# extract dge and create dge matrix
dge <- extract_dge(umi_obs_tpt)
dge_mat <- create_dge_matrix(dge)

# output file
if (tools::file_ext(snakemake@output[[1]]) == "gz") {
  outfile <- gzfile(here(snakemake@output[[1]]))
}else{
  outfile <- here(snakemake@output[[1]])
}

# save dge matrix to output file
write.table(dge_mat, file = outfile, sep = "\t", row.names = FALSE, quote = FALSE)

/* Define a margin before every image element */
img {
  margin-top: 3em;
}

# load processed de output (includes confidence level and distance to TSS)
results <- here(snakemake@input$processed_results) %>%
  read.csv(stringsAsFactors = FALSE) %>%
  mutate(sample = paste0("Chr", sub("iScreen1", "", sample)),
         significant = if_else(pval_adj_allTests < 0.05, true = "sig", false = "non_sig"))

# encode chip-seq files
chip_files <- here(snakemake@input$encode_chip)
names(chip_files) <- sub("_.*", "", basename(chip_files))

# function to calculate average assay signal per enhancer
avg_assay_signal <- function(assay, enh_coords) {
  overlaps <- findOverlaps(enh_coords, assay)
  overlaps <- split(overlaps, from(overlaps))
  names(overlaps) <- names(enh_coords[as.numeric(names(overlaps))])
  sapply(overlaps, function(x) {
    sum(mcols(assay[to(x)])$score * width(assay[to(x)])) / sum(width(assay[to(x)]))
  })
}

# calculate average chromatin signal at enhancers
enh_chrom_signal <- chip_data %>%
  sapply(FUN = avg_assay_signal, enh_coords = enh_coords) %>%
  as.data.frame() %>%
  rownames_to_column(var = "perturbation")

# transform to data.frame and merge with cis enhancer - gene pairs
cis_enh_perts <- cis_enh_perts %>%
  select(sample, perturbation, enh_chr, gene, pvalue, prop_grna_hits, dist_to_tss, significant) %>%
  left_join(enh_chrom_signal, by = "perturbation")

# function to randomly draw n non-significant pairs for a significant association. any enhancers
# from drawn non-significant pairs are excluded from subsequent draws
draw_control_pairs <- function(sig, non_sig, n = 3) {

  # get distance to tss bins for which a control pair needs to be drawn (repeat n times)
  bins <- select(sig, enh_chr, dTSS_bin)
  bins <- bins[rep(seq_len(nrow(bins)), each = n), ]

  # iteratively select a random negative control for each bin and remove any pairs involving that
  # enhancer from pool of non-significant pairs for any subsequenct draws
  random_ctrls <- data.frame()
  for (i in seq_len(nrow(bins))) {

    # get pairs on same chromosome and same distance to TSS bin
    nonSig_bin <- filter(non_sig, enh_chr == bins[i, "enh_chr"],
                         dTSS_bin == bins[i, "dTSS_bin"])

    # randomly draw a control pair (if possible)
    if (nrow(nonSig_bin)) {

      ctrl <- sample_n(nonSig_bin, size = 1)
      random_ctrls <- rbind(random_ctrls, ctrl)

      # remove any pairs involving that enhancer from non_sig
      non_sig <- filter(non_sig, perturbation != ctrl$perturbation)

    }else{
      warning("Not enough unique enhancers to draw control!\n", call. = FALSE)
    }
  }

  return(random_ctrls)

}

# create distance to TSS bins
dtss_range <- range(cis_enh_perts$dist_to_tss)
bins <- seq(from = dtss_range[1], to = dtss_range[2], by = 1e4)

# bin distance to tss
cis_enh_perts <- cis_enh_perts %>%
  mutate(dTSS_bin = cut(dist_to_tss, bins))

# extract significant enhancer - gene associations
sig_pairs <- filter(cis_enh_perts, significant == "sig")

# get other tested enhancer - gene pairs but whose enhancers appear only in non-sigificant tests
nonSig_pairs <- cis_enh_perts %>%
  filter(significant == "non_sig", !perturbation %in% sig_pairs$perturbation,
         perturbation != "GATA1")  # no GATA1, because it's the only on chrX

# randomly draw 3 control pairs for each significant enhancer - gene association
set.seed(20190617)
nonSig_ctrls <- draw_control_pairs(sig_pairs, non_sig = nonSig_pairs, n = 3)

# combine significant and control pairs
all_pairs <- bind_rows(sig_pairs, nonSig_ctrls)

# plot distance to TSS for significant pairs and sampled non-significant control pairs
ggplot(all_pairs, aes(x = dist_to_tss, fill = significant)) +
  facet_wrap(~sample) +
  geom_histogram(bins = 30) +
  labs(x = "Distance to TSS", y = "pairs",
       title = "Distance to TSS of significant vs non-significant control pairs") +
  theme_bw()

# figure caption
cap <- "Distance to TSS distribution of significant enhancers and sampled non-significant controls"

# extract unique enhancers from gene-enhancer pairs and convert to long format
enh_plot <- all_pairs %>%
  select(-c(sample, gene, pvalue, prop_grna_hits, dist_to_tss, dTSS_bin)) %>%
  distinct() %>%
  gather(key = "assay", value = "fc_over_ctrl", -c(perturbation, enh_chr, significant))

# pairwise comparisons using wilcoxon test and format adjusted p-value for plotting
pw_tests <- compare_means(fc_over_ctrl ~ significant, group.by = "assay", data = enh_plot,
                          method = "wilcox.test", p.adjust.method = "holm") %>%
  mutate(p.adj.plot = if_else(p.adj < 0.001, true = "< 0.001",
                              false = paste("=", round(p.adj, digits = 2))))

# plot difference in enhancer relevant chromatin marks between significant and non-significant hits
ggplot(enh_plot, aes(x = significant, y = fc_over_ctrl + 1, fill = significant)) +
  facet_wrap(~assay, ncol = 2) +
  geom_boxplot(notch = TRUE) +
  stat_pvalue_manual(pw_tests, label = "P adj. {p.adj.plot}", y.position = log10(145),
                     inherit.aes = FALSE) +
  labs(y = "Chromatin FC over control") +
  scale_x_discrete(labels = c("Non significant", "Significant")) +
  scale_y_log10(limits = c(NA, 180)) +
  scale_fill_manual(values = c(non_sig = "darkgray", sig = "steelblue")) +
  theme_bw() + 
  theme(axis.title.x = element_blank(), legend.position = "none", text = element_text(size = 15),
        panel.grid = element_blank())

# figure caption
cap <- "Chromatin marks of significant enhancers vs non-significant control enhancers"