fig_01_genome_wide_marks_chr_x.Rmd

---
title: "Figure 1. Genome-wide and chromosome X histone marks analysis."
output:
  workflowr::wflow_html:
    code_folding: show
editor_options:
  chunk_output_type: console
params:
  datadir: "./data/"
  chromhmm: "./data/bed/E008_15_coreMarks_hg38lift_noblock_fullnames.bed"
---

# Summary

Supplementary code for panel 1 figures.

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
library(ggplot2)
library(tidyverse)
library(wigglescout)
library(ggpubr)
library(ggrastr)
library(knitr)
library(treemap)
library(ggridges)
library(karyoploteR)
library(cowplot)
library(rtracklayer)
library(ggrepel)
library(scales)

source("./code/embed_functions.R")
source("./code/globals.R")

my_svg <- function(file, width, height) {
  library(svglite)
  svglite(file = file, width = width, height = height, bg = "white")
}

dir.create("./figures_data")
dir.create("./output")

knitr::opts_chunk$set(dev = c('png', 'my_svg'), fig.ext = c("png", "svg"), fig.width = 6, fig.height = 6)
knitr::opts_chunk$set(class.source='fold-show')
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
```


## Helper functions

```{r helper-functions, class.source = NULL}
#' Calculate INRC from a mapped read counts table, and append such values
#' to it.
#'
#' @param counts Counts table. Corresponding file is provided as part of the 
#'   included metadata.
#' @param selector Counts column used. Final_mapped represents the final number
#'   of reads after deduplication and blacklisting.

#' @return A table including INRC and INRC norm to naive reference
calculate_inrc <- function(counts, selector = "final_mapped") {
  counts$condition <- paste(counts$celltype, counts$treatment, sep="_")
  
  inputs <- counts[counts$ip == "Input", c("library", selector)]
  colnames(inputs) <- c("library", "input_reads")
  
  non_inputs <- counts[counts$ip != "Input",]
  counts <- merge(non_inputs, inputs, by.x="input", by.y="library")
  counts$inrc <- counts[, selector] / counts[, "input_reads"]
  
  references <- counts[grepl("_Ni_pooled", counts$library), c("ip", "inrc")]
  colnames(references) <- c("ip", "ref_inrc")
  
  counts <- merge(counts, references, by="ip")
  counts$norm_to_naive <- counts$inrc / counts$ref_inrc
  
  id_vars <- c("ip", "treatment", "celltype", "condition", "replicate", "norm_to_naive")
  inrc <- counts[, c(id_vars)]
  
  inrc$condition <- factor(
    inrc$condition,
    levels = c(
      "Naive_Untreated",
      "Primed_Untreated",
      "Naive_EZH2i",
      "Primed_EZH2i"
    )
  )
  
  inrc
}

#' Barplot INRC pooled vs replicates per condition
#'
#' @param inrc Table with the INRC values
#' @param ip Which IP to plot
#' @param colors Corresponding colors
inrc_barplot <- function(inrc, ip, colors, font = 16) {
  inrc <- inrc[inrc$ip == ip, ]
  
  # So paired test takes right replicates
  inrc <- inrc[order(inrc$condition, inrc$replicate), ]
  
  max_v <- max(abs(inrc$norm_to_naive))
  
  aesthetics <- aes(x = .data[["condition"]],
                    y = .data[["norm_to_naive"]],
                    color = .data[["condition"]])
  
  my_comp <- list(c("Naive_Untreated", "Primed_Untreated"),
                  c("Naive_Untreated", "Naive_EZH2i"),
                  c("Primed_Untreated", "Primed_EZH2i"),
                  c("Naive_EZH2i", "Primed_EZH2i"))
  
  stats_method <- "t.test"
  
  ggplot(inrc[inrc$replicate != 'pooled',], aesthetics) +
    geom_point() +
    stat_compare_means(
      method = stats_method,
      paired = FALSE,
      comparisons = my_comp,
      label = "p.format"
    ) +
    geom_bar(
      data = inrc[inrc$replicate == 'pooled',],
      stat = 'identity',
      alpha = 0.6,
      aes(fill = condition)
    ) +
    scale_fill_manual(values = colors) +
    scale_color_manual(values = colors) +
    labs(
      x = "",
      y = 'INRC fraction vs Naïve',
      title = paste(ip, "MINUTE-ChIP"),
      caption = paste(stats_method, "signif. test, paired")
    ) +
    theme_classic(base_size = font) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    ylim(0, 1.5)
}


#' Summarize stats per chromosome on a scaled bigWig file
#'
#' @param bwfile BigWig file to summarize
#' @param chromosomes Array of chromosome names to include.
#'
#' @return A data frame with stats per chromosome: mean, chr size, #reads
#'   (estimated as (score * chr size) / fraglen), %reads.
scaled_reads_per_chromosome <- function(bwfile, chromosomes, fraglen = 150) {
  granges <- unlist(summary(BigWigFile(bwfile)))
  df <- data.frame(granges[seqnames(granges) %in% chromosomes, ])
  rownames(df) <- df$seqnames
  
  # Calculate scaled number of reads as mean x chromosome length / read length
  df$nreads <- (df$score * df$width) / fraglen
  
  # Perc of total
  df$perc <- (df$nreads / sum(df$nreads)) * 100
  
  # Perc size
  df$size <- df$width / sum(df$width)

  df$group <- basename(bwfile)
  df[chromosomes, ]
}

chromosomes <- paste0("chr", c(1:22, "X"))

# Fix some parameters on treemap function to remove some clutter from nb.
chr_treeplot <- partial(
  treemap,
  index = "seqnames",
  vSize = "nreads",
  vColor = "score",
  type = "value",
  mapping = c(0, 3),
  range = c(0, 3),
  fontsize.labels = 16,
  fontsize.legend = 16,
  fontsize.title = 20
)


ridges_chromosome_plot <- function(values, column, color, main_seqs = chromosomes, scale = 1.7) {
  value_name <- column
  df <- values[values$seqnames %in% main_seqs, c("seqnames", value_name)]
  colnames(df) <- c("seqnames", "value")
  df$value <- as.numeric(df$value)
  df$seqnames <- factor(df$seqnames, levels = rev(main_seqs))

  df_summary <- df %>% group_by(seqnames) %>%
    summarise(value=median(value, na.rm = T))
  x_nudge <- quantile(df$value, 0.02, na.rm = T)

  ggplot(df, aes(x = value, y = seqnames, fill = seqnames)) +
    geom_density_ridges(
      rel_min_height = 0.001,
      scale = 1.7,
      calc_ecdf = TRUE,
      quantile_lines = TRUE, quantiles = 2,
    ) +
    theme_default(base_size = 12) +
    labs(y = "", x = "log2FC") +
    scale_fill_manual(values = c(color, rep("#bbbbbb", 22))) + theme(legend.position = "none") +
    geom_vline(xintercept = 0, linetype = "dashed", size = 0.2) +
    geom_text(data=df_summary,
              aes(label=sprintf("%1.2f", value)),
              position=position_nudge(y=0.35, x = x_nudge), colour="black", size=3)
}

get_long_format_heatmap_data <- function(df, mark) {
  columns <- grep("mean_cov", colnames(df), value = T)
  main_seqs <- paste0("chr", c(1:22, "X"))
  df <- df[df$seqnames %in% main_seqs, c("seqnames", columns)]
  
  summary_mat <- df %>% group_by(seqnames) %>% summarise_at(columns, mean, na.rm = TRUE)
  
  to_plot <- summary_mat %>% 
    select("seqnames", contains(mark) & contains("mean_cov") & !contains("rep"))
  
  # Reorder chromosomes and conditions
  conditions <- c("Ni", "Ni_EZH2i", "Pr", "Pr_EZH2i")
  to_plot$seqnames <- gsub("chr", "", to_plot$seqnames)
  to_plot$seqnames <- factor(to_plot$seqnames, levels = c(1:22, "X"))
  colnames(to_plot) <- c("seqnames", conditions)
  to_plot_melt <- pivot_longer(to_plot, !seqnames)
  to_plot_melt$name <- factor(to_plot_melt$name, levels = rev(conditions))
  to_plot_melt
}

colors_list <- c("Naive_EZH2i"="#82c5c6",
                 "Naive_Untreated"="#278b8b",
                 "Primed_EZH2i"="#f49797",
                 "Primed_Untreated"="#f44b34")

```

## Data input

```{r table-data-input}
genes <- read.table("./data/meta/Kumar_2020_master_gene_table_rnaseq_shrunk_annotated.tsv",
  header = T, sep = "\t",
  colClasses = c(rep("character", 5), rep("factor", 3), rep("numeric", 86)))

bins <- read.table("./data/meta/Kumar_2020_master_bins_10kb_table_final_raw.tsv",
  header = T, sep = "\t",
  colClasses = c(rep("character", 5), rep("numeric", 112)))

counts_file <- file.path(params$datadir, "meta", "Kumar_2020_stats_summary.csv")
counts <- read.table(counts_file, sep="\t", header = T, na.strings = "NA", stringsAsFactors = F)
```

# Global quantification

## Read counts

### H2AUb

```{r h2aub-inrc-barplot, fig.height = 8, fig.width = 8}
inrc <- calculate_inrc(counts) 
inrc_barplot(inrc, "H2Aub", colors_list)
```

You can download data values here: `r embed_last_plot_data(name="H2AUb_global.tsv")`.

### H3K27m3

```{r h3k27m3-inrc-barplot, fig.height = 8, fig.width = 8}
inrc_barplot(inrc, "H3K27m3", colors_list)
```

You can download data values here: `r embed_last_plot_data(name="H3K27m3_global.tsv")`.

### H3K4m3

```{r h3k4m3-inrc-barplot, fig.height = 8, fig.width = 8}
inrc_barplot(inrc, "H3K4m3", colors_list)
```

You can download data values here: `r embed_last_plot_data(name="H3K4m3_global.tsv")`.

## ChromHMM

Here global average per ChromHMM categories are shown. 

### H3K27m3

```{r chromhmm-h3k27m3, fig.width = 13, fig.height = 14}
labels <- gsub("_pooled.hg38.scaled.bw", "", basename(bwfiles$k27))
labels <- gsub("_H9", "", labels)

chromhmm <- params$chromhmm

plot_bw_loci_summary_heatmap(bwfiles$k27, chromhmm, labels = labels, remove_top=0.001)
```

You can download data values here: `r embed_last_plot_data(name="H3K27m3_chromhmm.tsv")`.

### H3K4m3

```{r chromhmm-h3k4m3, fig.width = 13, fig.height = 14}
labels <- gsub("_pooled.hg38.scaled.bw", "", basename(bwfiles$k4))
labels <- gsub("_H9", "", labels)

plot_bw_loci_summary_heatmap(bwfiles$k4, chromhmm, labels = labels, remove_top=0.001)
```

You can download data values here: `r embed_last_plot_data(name="H3K4m3_chromhmm.tsv")`.

### H2AUb

```{r chromhmm-h2aub, fig.width = 13, fig.height = 14}
labels <- gsub("_pooled.hg38.scaled.bw", "", basename(bwfiles$ub))
labels <- gsub("_H9", "", labels)

plot_bw_loci_summary_heatmap(bwfiles$ub, chromhmm, labels = labels, remove_top=0.001)
```

You can download data values here: `r embed_last_plot_data(name="H2Aub_chromhmm.tsv")`.

## Genome-wide 10kb bins

### H3K27m3 Naïve vs Primed

```{r h3k27m3-ni-vs-pr-10kb-bins, warning = F}
points_color <- "#112233"
points_shape <- "."
raster_dpi <- 300

ggplot(bins, aes(x=log2(H3K27m3_Ni_mean_cov), y=log2(H3K27m3_Pr_mean_cov))) +
    rasterise(geom_point(size = 1, alpha = 0.2, color = points_color, shape = points_shape), dpi = raster_dpi) +
    geom_density2d(binwidth = 0.1) +
    geom_abline(slope = 1, linetype = "dashed") +
    theme_default() +
    labs(x = "Log2 H3K27m3 Naïve - FPGC",
         y = "Log2 H3K27m3 Primed - FPGC", title = "H3K27m3", subtitle = "10kb bins") + 
    coord_cartesian(xlim = c(-8, 8), ylim = c(-8, 8))
```

### H2AUb Naïve vs Primed

```{r h2aub-ni-vs-pr-10kb-bins, warning = F}
ggplot(bins, aes(x=log2(H2Aub_Ni_mean_cov), y=log2(H2Aub_Pr_mean_cov))) +
    rasterise(geom_point(size = 1, alpha = 0.2, color = "#2f1547", shape = points_shape), dpi = raster_dpi) +
    geom_abline(slope = 1, linetype = "dashed") +
    geom_density2d(binwidth = 0.1) +
    theme_default() +
    labs(x = "Log2 H2Aub Naïve - FPGC",
         y = "Log2 H2Aub Primed - FPGC", title = "H2Aub", subtitle = "10kb bins") + 
    coord_cartesian(xlim = c(-8, 8), ylim = c(-8, 8))
```

### H3K4m3 Naïve vs Primed

```{r h3k4m3-ni-vs-pr-10kb-bins, warning = F}
ggplot(bins, aes(x=log2(H3K4m3_Ni_mean_cov), y=log2(H3K4m3_Pr_mean_cov))) +
    rasterise(geom_point(size = 1, alpha = 0.2, color = "#614925", shape = points_shape), dpi = raster_dpi) +
    geom_abline(slope = 1, linetype = "dashed") +
    geom_density2d(binwidth = 0.1) +
    theme_default() +
    labs(x = "Log2 H3K4m3 Naïve - FPGC",
         y = "Log2 H3K4m3 Primed - FPGC", title = "H3K4m3", subtitle = "10kb bins") + 
    coord_cartesian(xlim = c(-8, 8), ylim = c(-8, 8))
```

# Chromosome X

## 10kb bin per chromosome plot

```{r chromosomes-jitter, warning=F, message=F}
columns <- grep("mean_cov", colnames(bins), value = T)
main_seqs <- paste0("chr", c(1:22, "X"))
df <- bins[bins$seqnames %in% main_seqs, c("seqnames", columns)]

summary_mat <- df %>% group_by(seqnames) %>% summarise_at(columns, mean, na.rm = TRUE)

to_plot <- summary_mat %>% select("seqnames", contains("mean_cov") & !contains("IN") & !contains("rep"))
# Reorder chromosomes
to_plot$seqnames <- gsub("chr", "", to_plot$seqnames)
to_plot$seqnames <- factor(to_plot$seqnames, levels = c(1:22, "X"))
to_plot_melt <- pivot_longer(to_plot, !seqnames)

to_plot_melt$name <- gsub("_mean_cov", "", to_plot_melt$name)
to_plot_melt$ip <- str_split_fixed(to_plot_melt$name, "_", 2)[, 1]
to_plot_melt$condition <- str_split_fixed(to_plot_melt$name, "_", 2)[, 2]

ggplot(to_plot_melt, aes(color = ip, x = condition, y = value, label = seqnames)) +
  geom_boxplot(color = "gray", alpha = 0.9) +
  geom_jitter(position = "dodge") +
  geom_jitter(data = to_plot_melt %>% filter(seqnames == "X"), color = "black", size = 3.5, position = "dodge") +
  geom_label_repel(data = to_plot_melt %>% filter(
    seqnames == "X" & ip == "H3K27m3" & condition == "Ni"), color = "black", box.padding = 1.5) +
  theme_default(base_size=12) + 
  facet_wrap(. ~ ip, nrow = 1) +
  geom_hline(yintercept = 1, linetype = "dashed", alpha = 0.8) +
  scale_color_manual(
    values = c("H2Aub" = gl_mark_colors$H2Aub,
               "H3K27m3" = gl_mark_colors$H3K27m3,
               "H3K4m3" = gl_mark_colors$H3K4m3)) +
  labs(y="FPGC", title =
         "Mean 10kb bin RPGC per chromosome and histone mark",
       subtitle = "Chromosome X in black")

```

## Whole chromosome-view of histone marks autosome vs chromosome X

These figures are made using the package karyoploteR: https://academic.oup.com/bioinformatics/article/33/19/3088/3857734

### H3K27m3

```{r h3k27m3-karyoplot, warning = F, message = F}
kp <-
  plotKaryotype(
    genome = "hg38",
    plot.type = 1,
    main = "H3K27m3 Naïve vs Primed",
    chromosomes = c("chr7", "chrX")
  )

d1 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$k27[[1]]),
  data.panel = 1,
  col = "#092ba8",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

d2 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$k27[[3]]),
  data.panel = 2,
  col =
    "#5d9ddd",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

# Store these values on a table
df1 <- cbind(data.frame(d1$latest.plot$computed.values$windows, d1$latest.plot$computed.values$density))
df2 <- cbind(data.frame(d2$latest.plot$computed.values$windows, d2$latest.plot$computed.values$density))
df_karyo <- left_join(df1, df2)
colnames(df_karyo) <- c("seqnames", "start", "end", "width", "strand", "H3K27m3_Pr", "H3K27m3_Ni")
final_df <- df_karyo
```

### H3K4m3
```{r h3k4m3-karyoplot, warning = F, message = F}
kp <-
  plotKaryotype(
    genome = "hg38",
    plot.type = 1,
    main = "H3K4m3 Naïve vs Primed",
    chromosomes = c("chr7", "chrX")
  )

d1 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$k4[[3]]),
  data.panel = 2,
  col =
    "#ffab45",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

d2 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$k4[[1]]),
  data.panel = 1,
  col = "#e76e3b",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

# Store these values on a table
df1 <- cbind(data.frame(d1$latest.plot$computed.values$windows, d1$latest.plot$computed.values$density))
df2 <- cbind(data.frame(d2$latest.plot$computed.values$windows, d2$latest.plot$computed.values$density))
df_karyo <- left_join(df1, df2)
colnames(df_karyo) <- c("seqnames", "start", "end", "width", "strand", "H3K4m3_Pr", "H3K4m3_Ni")
final_df <- left_join(final_df, df_karyo)
```


### H2Aub
```{r h2aub-karyoplot, warning = F, message = F}
kp <-
  plotKaryotype(
    genome = "hg38",
    plot.type = 1,
    main = "H2Aub Naïve vs Primed",
    chromosomes = c("chr7", "chrX")
  )

d1 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$ub[[1]]),
  data.panel = 1,
  col = "#400c84",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

d2 <- kpPlotDensity(
  kp,
  rtracklayer::import(bwfiles$ub[[3]]),
  data.panel = 2,
  col =
    "#a07af0",
  chromosomes = c("chr7", "chrX"),
  window.size = 500000
)

# Store these values on a table
df1 <- cbind(data.frame(d1$latest.plot$computed.values$windows, d1$latest.plot$computed.values$density))
df2 <- cbind(data.frame(d2$latest.plot$computed.values$windows, d2$latest.plot$computed.values$density))
df_karyo <- left_join(df1, df2)
colnames(df_karyo) <- c("seqnames", "start", "end", "width", "strand", "H2Aub_Pr", "H2Aub_Ni")
final_df <- left_join(final_df, df_karyo)
```

## H3K27m3 naïve per chromosome treemap

H3K27me3 is highly abundant on X chromosome on naïve cells.

If we take a look at coverage per chromosome for both Naïve and Primed cells:

```{r treemap-h3k27m3-naive, fig.width=8, fig.height=7, warning = FALSE}
values <- scaled_reads_per_chromosome(bwfiles$k27[[1]], chromosomes = chromosomes)

chr_treeplot(
  values,
  palette = c("#ffffff", gl_condition_colors[["Naive_Untreated"]]),
  fontcolor.labels = "#555555",
  border.col = c("white"),
  title = "H3K27m3 - Naïve"
)
```

Values can be downloaded here: `r embed_df(values[, c("seqnames", "nreads", "score")])`.

In this and subsequent plots, each rectangle's size is proportional to the
number of reads mapped to its corresponding chromosome. Color intensity
represents mean coverage per chromosome, and rectangles are ordered according
to size. Top-left is the highest value.

## Per chromosome distribution of logFC values

These figures are made using the package ggridges: https://wilkelab.org/ggridges/

### Primed vs Naive RNA-seq

```{r rnaseq-per-chr-pr-vs-ni-ridgeplot, warning = F, message = F}
ridges_chromosome_plot(genes, "RNASeq_DS_Pr_vs_Ni_log2FoldChange", "#F08080") +
  labs(title = "RNA Seq log2FC distribution Primed vs Naïve DESeq2") +
  coord_cartesian(xlim=c(-7, 7))
```

`r embed_df(genes[, c("seqnames", "RNASeq_DS_Pr_vs_Ni_log2FoldChange")], name = "rnaseq_ridgeplot.tss", text = "Values")` used in this plot.


### Primed vs Naive H3K27m3

```{r h3k27m3-per-chr-pr-vs-ni-ridgeplot, warning = F, message = F}
ridges_chromosome_plot(genes, "H3K27m3_DS_Pr_vs_Ni_log2FoldChange", "#3e5aa8") +
  labs(title = "H3K27m3 Primed vs Naïve DESeq2") + coord_cartesian(xlim=c(-7, 7))
```

`r embed_df(genes[, c("seqnames", "H3K27m3_DS_Pr_vs_Ni_log2FoldChange")], name = "rnaseq_ridgeplot.tss", text = "Values")` used in this plot.

### EZH2i vs Naive RNA-seq data

```{r rnaseq-per-chr-ezh2i-vs-ni-ridgeplot, warning = F, message = F}
ridges_chromosome_plot(genes, "RNASeq_DS_EZH2i_vs_Ni_log2FoldChange", "#F08080") +
  labs(title = "RNA Seq EZH2i vs Naïve DESeq2") +
  coord_cartesian(xlim=c(-5, 5))
```

`r embed_df(genes[, c("seqnames", "RNASeq_DS_EZH2i_vs_Ni_log2FoldChange")], name = "rnaseq_ridgeplot.tss", text = "Values")` used in this plot.