07_CTCF_Threshold_mm.Rmd

---
title: "CTCF threshold analysis"
author: "Mikhail Dozmorov"
date: "`r Sys.Date()`"
output:
  pdf_document:
    toc: no
  html_document:
    theme: cerulean
    toc: yes
---

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
# Set up the environment
library(knitr)
opts_chunk$set(cache.path='cache/', fig.path='img/', cache=F, tidy=T, fig.keep='high', echo=F, dpi=100, warnings=F, message=F, comment=NA, warning=F, results='as.is', fig.width = 10, fig.height = 6) #out.width=700, 
library(pander)
panderOptions('table.split.table', Inf)
set.seed(1)
```

# Libraries

```{r libraries}
library(tidyverse)
library(readxl)
library(writexl)
library(cowplot)
library(stringr)
library("ggsci")
library(scales)
# scales::show_col(pal_lancet("lanonc")(8))
mycols = pal_lancet("lanonc")(8)
# Color palette for the heatmap, https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf
col3 <- colorRampPalette(c('blue', 'white', 'red'))(20)
# col3 <- colorRampPalette(c('blue', 'gray', 'yellow'))(20)
# col3 <- colorRampPalette(c('green', 'black', 'red'))(20)
# col3 <- colorRamps::green2red(n = 20)
library(GenomicRanges)
```

# Settings

```{r settings}
dir_data <- "/Users/mdozmorov/Documents/Data/GoogleDrive/CTCF.dev"
# Project folder path
dir_project <- "/Users/mdozmorov/Documents/Work/GitHub/CTCF.dev"

# Input files
fileNameIn1 <- file.path(dir_project, "RData", "mm10.MA0139.1.bed")
fileNameIn2 <- file.path(dir_project, "data", "mm10-CTCF.bed")
# Output files
fileNameOut1 <- file.path(dir_project, "results", "Figure_mouse_pvalues_threshold.svg")
```

# Load data

```{r data}
# Read in data
mtx <- read_tsv(fileNameIn1, col_names = FALSE, col_types = c("cddcdcdd"))
gr_meme <- GRanges(seqnames = mtx$X1, ranges = IRanges(start = mtx$X2, end = mtx$X3), strand = mtx$X6)
gr_meme$pval <- mtx$X8
gr_meme$qval <- mtx$X9

mtx1 <- read_tsv(fileNameIn2, col_names = FALSE, col_types = c("cddccc"))
gr_encode <- GRanges(seqnames = mtx1$X1, ranges = IRanges(start = mtx1$X2, end = mtx1$X3))
```

## Summary statistics

```{r}
print("Summary of full MEME p-value distribution")
summary(gr_meme$pval)

print(paste("Length of full MEME CTCF set:", length(gr_meme)))
print(paste("Length of ENCODE SCREEN CTCF set:", length(gr_encode)))
```

## Scanning the threshold

```{r}
# thresholds <- seq(from = 1E-8, to = 1E-4, length.out = 10) %>% rev()
thresholds <- c(1.00E-4, 5.00E-4, 1.00E-5, 5.00E-5, 1.00E-6, 5.00E-6, 1.00E-7, 5.00E-7, 1.00E-8)
# formatC(thresholds, format = "e")

gr_meme_pos_list <- vector(mode = "numeric", length = length(thresholds))
gr_meme_neg_list <- vector(mode = "numeric", length = length(thresholds))
for (i in 1:length(thresholds)) {
  pval_cutoff <- thresholds[i]
  gr_meme_cutoff <- gr_meme[gr_meme$pval < pval_cutoff]
  gr_meme_overlaps <- countOverlaps(gr_meme_cutoff, gr_encode)
  gr_meme_pos_list[i] <- sum(gr_meme_overlaps > 0) / length(gr_meme_overlaps)
  gr_meme_neg_list[i] <- sum(gr_meme_overlaps == 0) / length(gr_meme_overlaps)
}

gr_summary <- data.frame(Threshold = thresholds, TP = gr_meme_pos_list, FP = gr_meme_neg_list)

mtx_to_plot <- pivot_longer(gr_summary, cols = c("TP", "FP"), names_to = "Type")
p <- ggplot(mtx_to_plot, aes(x = Threshold, y = value, color = Type)) +
  geom_line() +
  scale_x_continuous(breaks = thresholds, labels = scientific) +
  scale_y_continuous(breaks = seq(0, 1, 0.2))   + 
  scale_color_manual(values = mycols[1:2]) +
  coord_trans(x ="log10") +
  geom_vline(xintercept = 1E-6, linetype = "dashed", color = "black") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      # axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
      plot.title = element_text(hjust = 0.5),
      legend.text = element_text(size = 6)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  xlab("FIMO p-value threshold") +
  ylab("FP/TP rate")

# ggsave(filename = fileNameOut1, width = 4.5, height = 1.5)
dev.off()
svg(filename = fileNameOut1, width = 4.5, height = 1.5)
print(p)
dev.off()

```