EDA_PWMScan.Rmd

---
title: "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(rtracklayer)
library(stringr)
library(cowplot)
library(stringr)
library("ggsci")
library(scales)
# scales::show_col(pal_lancet("lanonc")(8))
mycols = pal_lancet("lanonc")(8)
library(plyranges)
library(pheatmap)
library(patchwork)
library(GenometriCorr)
```

# Settings

```{r settings}
# Project folder path
dir_data <- "/Users/mdozmorov/Documents/Work/GitHub/CTCF.dev/"
# Genome assembly
genome_assembly <- "hg38"
genome_assembly <- "mm10"

# Input files
fileNameIn1 <- file.path(dir_data, "")
# Output files
fileNameOut1 <- file.path(dir_data, paste0("results/Table_stats_", genome_assembly, ".csv"))
fileNameOut2 <- file.path(dir_data, paste0("results/Figure_stats_", genome_assembly, ".svg"))
fileNameOut3 <- file.path(dir_data, paste0("results/Figure_jaccard_", genome_assembly, ".svg"))
```

# Load data

```{r data}
# Read in data
files <- list.files(path = file.path(dir_data, "data/PWMscan/"), pattern = genome_assembly, full.names = TRUE)
# Object to store GRanges 
bed_list_all <- list() # All chromosomes
bed_list_standard <- list() # Standard chromosomes
# Additional columns
extraCols_names <- c("character", "character")
names(extraCols_names) <- c("motif", "pval")
# Process each file
for(i in 1:length(files)) {
  print(basename(files[i]))
  # Read data
  bed_file <- import(files[i], format = "bed", genome = genome_assembly, colnames = c("chrom", "start", "end", "name", "score", "strand", "motif", "pval"), extraCols = extraCols_names)
  # Save all data
  bed_list_all <- c(bed_list_all, list(bed_file))
  # Keep standard chromosomes
  bed_file_standard <- keepStandardChromosomes(bed_file)
  # Save filtered data
  bed_list_standard <- c(bed_list_standard, list(bed_file_standard))
}
# Add names
names(bed_list_all) <- basename(files) %>% str_replace(., ".bed", "") %>% str_replace(., paste0(genome_assembly, ".PWMScan."), "")
names(bed_list_standard) <- basename(files) %>% str_replace(., ".bed", "") %>% str_replace(., paste0(genome_assembly, ".PWMScan."), "")
```

# Summaries

## Genomewide Number

```{r fig.height=3, fig.width=10}
stats_number_pos <- c()
stats_number_neg <- c()
for(i in 1:length(bed_list_standard)) {
  stats_number_pos <- c(stats_number_pos, length(bed_list_standard[[i]][strand(bed_list_standard[[i]]) == "+" ]))
  stats_number_neg <- c(stats_number_neg, length(bed_list_standard[[i]][strand(bed_list_standard[[i]]) == "-" ]))
}
mtx_to_plot <- data.frame(Number = c(stats_number_pos, stats_number_neg), 
                          List = c(names(bed_list_standard), names(bed_list_standard)), 
                          Strand = c(rep("+", length(stats_number_pos)), rep("-", length(stats_number_neg))))
# Wide format for saving
mtx_to_save <- pivot_wider(mtx_to_plot, names_from = "Strand", values_from = "Number")
mtx_to_save <- mtx_to_save %>% mutate(Total = `+` + `-`, .after = List) %>% arrange(desc(Total))
colnames(mtx_to_save) <- c("Database", "Total CTCF motifs", "Positive strand", "Negative strand")
write_csv(mtx_to_save, fileNameOut1)

stats_number_total <- aggregate(mtx_to_plot$Number, list(mtx_to_plot$List), sum)
colnames(stats_number_total) <- c("List", "Total")
mtx_to_plot$List <- factor(mtx_to_plot$List, levels = stats_number_total$List[order(stats_number_total$Total)])
ggplot(mtx_to_plot, aes(x = List, y = Number, fill = Strand)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  theme_bw(base_size = 15) +
  # get rid of the grid 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        text = element_text(size = 15),
        legend.position = c(0.8, 0.30)) +
  scale_fill_manual(values = mycols[1:2]) + # change colors 
  xlab("Database") +
  ylab("Number of motifs") +
  ggtitle(paste(genome_assembly, "genome assembly"))
ggsave(fileNameOut2, width = 10, height = 3)
```

## Chromosome-specific Number

```{r fig.height=16, fig.width=10}
bed_to_number_per_chr_plot <- function(i = 1) {
  bed_list_standard_selected = bed_list_standard[[i]]
  # Summarize counts per chromosome and strand
  mtx_to_plot <- bed_list_standard_selected %>% group_by(seqnames, strand) %>% summarise(Number = n()) %>% as.data.frame()
  # Reorder chromosomes
  mtx_to_plot$seqnames <- factor(mtx_to_plot$seqnames, levels = (unique(mtx_to_plot$seqnames)))
  p <- ggplot(mtx_to_plot, aes(x = seqnames, y = Number, fill = strand)) +
    geom_bar(stat = "identity") +
    # coord_flip() +
    theme_bw(base_size = 15) +
    # get rid of the grid 
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          text = element_text(size = 15),
          legend.position = c(0.8, 0.80)) +
    theme(axis.text.x = element_text(angle = 30, vjust = 0.5, hjust=1)) +
    scale_fill_manual(values = mycols[1:2]) + # change colors 
    xlab("Chromosome") +
    ylab("Number of motifs") +
    ggtitle(paste(genome_assembly, names(bed_list_standard)[i]))
  return(p)
}

p_list <- list()
for(i in 1: length(bed_list_standard)) {
  print(names(bed_list_standard)[i])
  p_list <- c(p_list, list(bed_to_number_per_chr_plot(i = i)))
}
wrap_plots(p_list, nrow = length(p_list))
```

## Jaccard / GenometriCorr

```{r jaccard, warning=FALSE, fig.width=6.5, fig.height=6}
save_pheatmap_svg <- function(x, filename, width=4.5, height=3, units = "in", res = 300) {
  stopifnot(!missing(x))
  stopifnot(!missing(filename))
  svg(filename, width=width, height=height)
  grid::grid.newpage()
  grid::grid.draw(x$gtable)
  dev.off()
}
# Jaccard calculations
jaccard <- function(gr_a, gr_b) {
  intersects <- GenomicRanges::intersect(gr_a, gr_b, ignore.strand = TRUE)
  intersection <- sum(width(intersects))
  union <- sum(width(GenomicRanges::union(gr_a, gr_b, ignore.strand = TRUE)))
  DataFrame(intersection, union, 
            jaccard = intersection/union,
             n_intersections = length(intersects))
}
# GenometriCorr of smaller query vs. larger reference
# Returns either Jaccard or ECDF area correlation
GCorr <- function(gr1, gr2, return_value = "relative.distances.ecdf.area.correlation") {
  # Select smaller query and larger reference
  # GenometriCorr itselv
  gr1_vs_gr2 <- GenometriCorrelation(query, reference, awhole.only = TRUE, showProgressBar = FALSE, permut.number = 0)
  # What to return
  if (return_value == "relative.distances.ecdf.area.correlation") {
    return(gr1_vs_gr2$awhole$relative.distances.ecdf.area.correlation)
  }
  if (return_value == "jaccard.measure") {
    return(gr1_vs_gr2$awhole$jaccard.measure)
  }
}

# Correlation matrix, empty
mtx_to_plot <- matrix(data = 0, nrow = length(bed_list_standard), ncol = length(bed_list_standard))
# Fill it in
for (i in 1:length(bed_list_standard)) {
  for (j in 1:length(bed_list_standard)) {
    # If diagonal, set to zero
    if (i == j) mtx_to_plot[i, j] <- 0
    # Process only one half, the other is symmetric
    if (i > j) {
      gr1 <- bed_list_standard[[i]]
      gr2 <- bed_list_standard[[j]]
      # Jaccard
      query <- gr1
      reference <- gr2
      mtx_to_plot[i, j] <- mtx_to_plot[j, i] <- jaccard(query, reference)[["jaccard"]]
      # GenometriCorr
      # if (length(gr1) >= length(gr2)) {
      #   reference <- gr1
      #   query <- gr2
      # } else {
      #   reference <- gr2
      #   query <- gr1
      # }
      # mtx_to_plot[i, j] <- mtx_to_plot[j, i] <- GCorr(query, reference)
    }
  }
}

# Trim row/colnames
rownames(mtx_to_plot) <- colnames(mtx_to_plot) <- str_trunc(names(bed_list_standard), width = 25) 
# Adjust clustering method
# if(genome_assembly == "hg38") {
#   clust_method <- "euclidean"
# }
# if(genome_assembly == "mm10") {
#   clust_method <- "correlation"
# }
clust_method <- "euclidean"
# Save the plot
# png("man/figures/excluderanges_hg38_jaccard.png", width = 1000, height = 900, res = 200)
p <- pheatmap(data.matrix(mtx_to_plot), cluster_cols = T, cluster_rows = T,
         clustering_method = "ward.D2",# "ward.D", 
         clustering_distance_rows = clust_method, 
         clustering_distance_cols = clust_method, 
         # annotation_row = mydf, annotation_colors = list(Group = c(PR = mycols[1], CR = mycols[2])),
         treeheight_row = 40,
         treeheight_col = 0, 
         display_numbers = TRUE)
save_pheatmap_svg(p, filename = fileNameOut3, width = 5.5, height = 4.5, units = "in", res = 300)
# dev.off()
```


```{r eval=FALSE}
# MDS
mds <- cmdscale(d = (1 - data.matrix(mtx_to_plot))) %>% as_tibble()
colnames(mds) <- c("Dim.1", "Dim.2")
mds <- mds %>% mutate(Groups = sample_annotation$Tumor, Samples = sample_annotation$Sample.Name)
mds$Groups <- factor(mds$Groups, levels = c("PR", "CR"))
ggplot(mds, aes(x = Dim.1, y = Dim.2)) + # , color = Groups, label = Samples
  geom_point(size = 2) +
  # scale_color_lancet() +
  scale_color_manual(values = mycols[1:2]) +
  # scale_color_brewer(palette = "Spectral") +
  theme_bw() +
  xlab("Coordinate 1") + ylab("Coordinate 2") 
#  geom_label_repel(show.legend = FALSE) #, force=1, box.padding=0.5, label.padding = 0.2, direction = c("y"), segment.color = 'grey50', max.overlaps = 15)
ggsave(file.path(dir_project, "manuscript/figures/figure_replicates/FigureS1A.svg"), width = 4, height = 3)
```