bkATAC_1allAnalysis.Rmd

---
title: | 
  | Script for analysing bulk ATAC-seq of reprogramming intermediate
author: |
  | Fernando J. Rossello
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  html_document:
    toc: true
    toc_depth: 2
    toc_float: 
      collapsed: false
fontsize: 12pt
pagetitle: "bkATAC"
---

# Preamble
## Things to note
- All input files can be downloaded from http://hrpi.ddnetbio.com/
- All input files are assumed to be in the data folder.
- For reproducibility purposes, the randomisation seed used throughout the manuscript was saved as an object and it should be loaded in all relevant sections (object `seed` in `data/seed.Rdata`loaded in each section's code).
- All genome browser tracks can be explored at http://hrpi.ddnetbio.com at the _bulk ATAC-seq iPSC reprog. section_.

## Load required libraries.

```{r session_setup, include = TRUE, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
library(dplyr)
library(edgeR)
library(irlba)
library(sjstats)
library(Mfuzz)
library(patchwork)
library(data.table)
library(ggplot2)
library(ggrepel)
library(knitr)
library(annotatr)
library(plyr)
```
## Define colour palettes and themes.

```{r color_labels_themes}

colorMedia = c("black","darkorange","blue","darkolivegreen2","red2","pink2")
names(colorMedia) = c("D0, 3, 7","Primed","t2iLGoY","5iLAF","NHSM","RSeT")
shapeTime = c(0,1,2,15,16,17,18)
names(shapeTime) = c("D0","D3","D7","D13","D21","P3","P10")

gplot.theme <- theme(panel.border=element_blank(), 
        panel.background=element_blank(),
        axis.line=element_line(),
        axis.line.y=element_line("black", size=1),
        axis.line.x=element_line("black", size=1),
        axis.text=(element_text(colour="black", size = 6, face="bold")),
        axis.title=(element_text(colour="black", size = 8, face="bold")),
        axis.ticks = (element_line("black", size=0.5)),
        legend.title = element_text(colour="black", size = 9, face="bold"),
        legend.text = element_text(colour="black", size= 8, face="bold"),
        legend.key.size=unit(0.1,"cm"),
        # legend.margin = unit(0, "cm"),
        legend.position="none",
        plot.margin=unit(c(0,2,0,0),"mm"),
        strip.text.x = element_text(size = 6, colour = "black", face = "bold", margin = margin(0,0,0,0, "cm"))
  )
```

# Exploratory Analysis - PCA.
## Import read counts of consensus peak set.
```{r import_read_counts}

peak.counts <-
  read.table(file = "data/bkATAC_consensus_peak_set_counts.txt",
             header = T,
             stringsAsFactors = F)
```

### Reorder samples and save peak information.
```{r reorder_samples}
order <- colnames(peak.counts)[7:ncol(peak.counts)]
order <- order[c(1, 2, 11:14, 3, 5, 7, 9, 19, 21, 15, 17, 4, 6, 8, 10, 20, 22, 16, 18)]
order <- factor(order, levels = order)

peak.info <- colnames(peak.counts)[1:6]
```

### Reorder peak counts table.

```{r prepare_counts_table}

peak.counts <- peak.counts[, c(peak.info, levels(order))]

# Change Geneid column name to PeakID

colnames(peak.counts)[1] <- "PeakID"
rownames(peak.counts) <- peak.counts$PeakID
```

## Within and between library normalization.

### Quantile-normalise log2(FPKM + 1)

```{r quantile_log2_fpkm}

peaks.fpkm <-
  rpkm(peak.counts[, 7:ncol(peak.counts)],
       log = F, gene.length = peak.counts$Length)

# Keep peaks with more than 5 FPKMs in at least 2 samples for exploratory analysis.

keep <- rowSums(peaks.fpkm > 5) >= 2

peaks.fpkm <- peaks.fpkm[keep, ]
peaks.fpkm <- log2(peaks.fpkm + 1)
peaks.fpkm <- normalizeQuantiles(peaks.fpkm)
```

## PCA 

```{r pca}

# Prepare annotation table.
pca.plot.annotation <- data.table(sample_id = colnames(peaks.fpkm))

pca.plot.annotation$sample_labels <- c(
  paste(rep("FM_", 3), rep("D", 3), rep(c(0, 3, 7), each = 2), rep(c("_R1", "_R2")), sep = ""),
  paste(rep("PM_", 2), rep("D", 2), rep(c(13, 21), each = 2), rep(c("_R1", "_R2")), sep = ""),
  paste(rep("PM_", 2), rep("P", 2), rep(c(3, 10), each = 2), rep(c("_R1", "_R2")), sep = ""),
  paste(rep("NM_", 2), rep("D", 2), rep(c(13, 21), each = 2), rep(c("_R1", "_R2")), sep = ""),
  paste(rep("NM_", 2), rep("P", 2), rep(c(3, 10), each = 2), rep(c("_R1", "_R2")), sep = "")
)

pca.plot.annotation$media = "D0, 3, 7"
pca.plot.annotation[grep("PM", sample_labels)]$media = "Primed"
pca.plot.annotation[grep("NM", sample_labels)]$media = "t2iLGoY"
pca.plot.annotation$media = factor(pca.plot.annotation$media, 
                           levels = c("D0, 3, 7", "Primed", "t2iLGoY"))
pca.plot.annotation$timept = tstrsplit(pca.plot.annotation$sample_labels, "_")[[2]]
pca.plot.annotation$timept = factor(pca.plot.annotation$timept, 
                            levels = unique(pca.plot.annotation$timept))

base::load("data/seed.Rdata")
.Random.seed <- seed
pca <- prcomp_irlba(t(peaks.fpkm))

pca.plot.annotation$PC1 = pca$x[,1]
pca.plot.annotation$PC2 = pca$x[,2]
pca.plot.annotation$PC3 = pca$x[,3]

ggplot(pca.plot.annotation, aes(PC1, PC2, color = media, label = timept)) +
  geom_point(size = 3) + geom_text_repel(size = 5) +
  scale_color_manual(values = colorMedia[1:3]) +
  scale_shape_manual(values = c(16, 15)) + theme_classic(base_size = 24)
  
ggplot(pca.plot.annotation, aes(PC1, PC3, color = media, label = timept)) +
  geom_point(size = 3) + geom_text_repel(size = 5) +
  scale_color_manual(values = colorMedia[1:3]) +
  scale_shape_manual(values = c(16, 15)) + theme_classic(base_size = 24)


```

# Fuzzy clustering.

## Aggregate sample counts and calculate FPKM.

```{r aggregate_sample_counts_fpkm}
# Prepare conditions labels.

conditions.labels <-
  c(
  "D0",
  "D3",
  "D7",
  "D13.*Primed",
  "D21.*Primed",
  "P3.*Primed",
  "P10.*Primed",
  "D13.*Smith",
  "D21.*Smith",
  "P3.*Smith",
  "P10.*Smith"
  )

# Aggregate peak counts.
aggregated.peak.counts <-
  sapply(conditions.labels, function(x)
    rowSums(peak.counts[, grep(x,
                               colnames(peak.counts),
                               value = TRUE)]))

colnames(aggregated.peak.counts) <-
  colnames(aggregated.peak.counts) %>% gsub("\\.\\*", "_", .)

```
## Quantile-normalise log2(RPKM + 1)

```{r aggregated_counts_fpkm_log2}

aggregated.peak.fpkm <-
  rpkm(aggregated.peak.counts,
       log = F,
       gene.length = peak.counts$Length)

# Keep peaks with more than 10 FPKMs in any condition.

aggregated.keep <- rowSums(aggregated.peak.fpkm > 10) >= 1

aggregated.peak.fpkm <- aggregated.peak.fpkm[aggregated.keep, ]
aggregated.peak.fpkm <- log2(aggregated.peak.fpkm + 1)
aggregated.peak.fpkm <- normalizeQuantiles(aggregated.peak.fpkm)

```

## Discard low coefficient of variation peaks.

```{r cv_20_peaks}
cv.20 <- apply(aggregated.peak.fpkm, 1, function(x) cv(x) * 100) > 20
cv.20 %>% table()

aggregated.peak.fpkm.high.cv <- aggregated.peak.fpkm[cv.20, ]

```

## Run fuzzy clustering.

```{r fuzzy_clustering}

exprSet <- ExpressionSet(assayData = aggregated.peak.fpkm.high.cv)

s.exprSet <- standardise(exprSet)

m <- mestimate(s.exprSet)
# Create seed as save it as an object.
# set.seed(1234)
# seed <- .Random.seed
# save(seed, file = "data/seed.Rdata")

load("data/seed.Rdata")
.Random.seed <- seed

cl.s.exprSet <- mfuzz(s.exprSet, c = 8, m = m)


# Retrieve peak cluster membership

peak.cluster.membership <- as.data.frame(cl.s.exprSet$cluster) %>% 
  dplyr::mutate(PeakID = rownames(.))
colnames(peak.cluster.membership)[1] <- "Cluster"

```
### Identify and subset high affinity peaks.

```{r core_peaks}

core.exprSet <- acore(s.exprSet,
                      cl = cl.s.exprSet,
                      min.acore = 0.8)

# Peak ids of high affinity peaks
core.exprSet.peaks.id <- sapply(c(1:8), function(x) as.vector(core.exprSet[[x]]$NAME)
)
core.exprSet.peaks.id <- as.character(unlist(core.exprSet.peaks.id))

```

## Create a table of high affinity cluster  peaks.
It contains z-scores and peak details.
Please note that for the benefit of biological flow, clusters have been relabeled as follows: 1 = "C 1", 2 = "C 2", 3 = "C 7", 4 = "C 4", 5 = "C 3", 6 = "C 6", 7 = "C 5", 8 = "C 8". This has also been indicated as `Cluster_Labels` accross the manuscript's source data and supplementary tables.
.

```{r high_affinity_peaks_z_score}

# Peaks z-scores

peak.z.scores <- as.data.frame(s.exprSet@assayData$exprs) %>% 
  dplyr::mutate(., PeakID = rownames(.))

# Table of high affinity peaks annotated.
high.affinity.peaks.annotated <- dplyr::left_join(peak.z.scores,
                                                  peak.cluster.membership,
                                                  by = "PeakID") %>%
  dplyr::left_join(peak.counts[, 1:6], by = "PeakID") %>%
  dplyr::filter(PeakID %in% core.exprSet.peaks.id)

high.affinity.peaks.annotated <-
  high.affinity.peaks.annotated %>% mutate(
  Cluster_Label = case_when(
  Cluster == "1" ~ "C 1",
  Cluster == "2" ~ "C 2",
  Cluster == "3" ~ "C 7",
  Cluster == "4" ~ "C 4",
  Cluster == "5" ~ "C 3",
  Cluster == "6" ~ "C 6",
  Cluster == "7" ~ "C 5",
  TRUE ~ "C 8"
  )
  )
  
```

## Cluster dynamics.

### High affinity peaks cluster counts.

```{r membership_counts}

core.exprSet.cluster.counts <- lapply(core.exprSet, nrow) %>%
  unlist() %>% as.data.frame() %>% 
  mutate(Cluster = c(1:8))

colnames(core.exprSet.cluster.counts)[1] <- "members"

core.exprSet.cluster.counts <-
  core.exprSet.cluster.counts %>% mutate(
  Cluster_Label = case_when(
  Cluster == "1" ~ "C 1",
  Cluster == "2" ~ "C 2",
  Cluster == "3" ~ "C 7",
  Cluster == "4" ~ "C 4",
  Cluster == "5" ~ "C 3",
  Cluster == "6" ~ "C 6",
  Cluster == "7" ~ "C 5",
  TRUE ~ "C 8"
  )
  )
  
core.exprSet.cluster.counts$Cluster_Label <-
  factor(x = core.exprSet.cluster.counts$Cluster_Label,
  levels = paste("C", 1:8, sep = " "))

knitr::kable(
  core.exprSet.cluster.counts %>% dplyr::select(Cluster_Label, members) %>% arrange(Cluster_Label)
  )
  
```

### Cluster dynamics of high affinity peaks.

```{r  cluster_dynamics_table}

high.affinity.peaks.melt <-
  melt(dplyr::select(
  high.affinity.peaks.annotated,
  -c("Chr", "Start", "End", "Strand", "Length", "Cluster_Label")
  ),
  id = c("PeakID", "Cluster"))

# Define stage and media

high.affinity.peaks.melt <-
  high.affinity.peaks.melt %>% mutate(stage = case_when(
  grepl("D13", variable) ~ "D13",
  grepl("D21", variable) ~ "D21",
  grepl("P3", variable) ~ "P3",
  grepl("P10", variable) ~ "P10",
  TRUE ~ as.character(variable)
  )) %>%
  mutate(media = case_when(
  grepl("D0|D3|D7", variable) ~ "FM",
  grepl("Smith", variable) ~ "NM",
  grepl("Primed", variable) ~ "PM"
  ))
  
# Order Stages                                                          
high.affinity.peaks.melt$stage <-
  factor(high.affinity.peaks.melt$stage,
  levels = c("D0", "D3", "D7", "D13", "D21", "P3", "P10"))
  
# Order Clusters
high.affinity.peaks.melt$Cluster <-
  factor(high.affinity.peaks.melt$Cluster, levels = c(1, 2, 5, 4, 7, 6, 3, 8))



cluster.dynamics.plot.values <-
  high.affinity.peaks.melt %>% group_by(media, stage, Cluster) %>% dplyr::summarise(
  mean = mean(value),
  sd = sd(value)
  )
  
# Define Cluster labels (please see comment above)

cluster.dynamics.plot.values <-
  cluster.dynamics.plot.values %>% mutate(
  Cluster_Label = case_when(
  Cluster == "1" ~ "C 1",
  Cluster == "2" ~ "C 2",
  Cluster == "3" ~ "C 7",
  Cluster == "4" ~ "C 4",
  Cluster == "5" ~ "C 3",
  Cluster == "6" ~ "C 6",
  Cluster == "7" ~ "C 5",
  TRUE ~ "C 8"
  )
  )

cluster.dynamics.plot.values$Cluster_Label  <-
  factor(
  x = cluster.dynamics.plot.values$Cluster_Label,
  levels = c("C 1", "C 2", "C 3", "C 4", "C 5", "C 6", "C 7", "C 8")
  )


```

Ribbons equal mean +/- SD (ggplot code hidden, see Rmd for details).

```{r cluster_dynamics_plot, fig.align="center", fig.height=9, fig.width=4.5, warning=FALSE, echo=FALSE}

# Write facet labels as described 
# here:https://stackoverflow.com/questions/3472980/how-to-change-facet-labels

cluster.labels.atac.seq <- list(
  "C 1" = "C 1 (n = 12024)",
  "C 2" = "C 2 (n = 7779)",
  "C 3" = "C 3 (n = 5077)",
  "C 4" = "C 4 (n = 3334)",
  "C 5" = "C 5 (n = 9117)",
  "C 6" = "C 6 (n = 10129)",
  "C 7" = "C 7 (n = 4885)",
  "C 8" = "C 8 (n = 7739)"
  )
  
cluster_labeller_atac <- function(variable,value){
  return(cluster.labels.atac.seq[value])
}

ggplot(data = NULL, aes(stage, mean, group = 1)) +
  geom_line(
  data = dplyr::filter(cluster.dynamics.plot.values, media %in% c("FM", "PM")),
  color = "#FF8C00",
  size = 0.5
  ) +
  geom_line(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "PM"),
  stage %in% c("D0", "D3", "D7")
  ),
  color = "black",
  size = 0.5
  ) +
  geom_ribbon(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "PM"),
  stage %in% c("D0", "D3", "D7")
  ),
  aes(
  x = stage,
  ymax = mean + sd,
  ymin = mean - sd
  ),
  fill = "black",
  alpha = 0.5
  ) +
  geom_ribbon(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "PM"),
  stage %in% c("D7", "D13", "D21", "P3", "P10")
  ),
  aes(
  x = stage,
  ymax = mean + sd,
  ymin = mean - sd
  ),
  fill = "#FF8C00",
  alpha = 0.5
  ) +
  geom_line(
  data = dplyr::filter(cluster.dynamics.plot.values, media %in% c("FM", "NM")),
  color = "#2302FE",
  size = 0.5
  ) +
  geom_line(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "NM"),
  stage %in% c("D7", "D13", "D21", "P3", "P10")
  ),
  color = "#2302FE",
  size = 0.5
  ) +
  geom_line(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "NM"),
  stage %in% c("D0", "D3", "D7")
  ),
  color = "black",
  size = 0.5
  ) +
  geom_ribbon(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "NM"),
  stage %in% c("D0", "D3", "D7")
  ),
  aes(
  x = stage,
  ymax = mean + sd,
  ymin = mean - sd
  ),
  fill = "black",
  alpha = 0.5
  ) +
  geom_ribbon(
  data = dplyr::filter(
  cluster.dynamics.plot.values,
  media %in% c("FM", "NM"),
  stage %in% c("D7", "D13", "D21", "P3", "P10")
  ),
  aes(
  x = stage,
  ymax = mean + sd,
  ymin = mean - sd
  ),
  fill = "#2302FE",
  alpha = 0.5
  ) +
  ylab("Chromatin accessibiliyt (Z-scaling)") +
  xlab("Reprogramming stages") +
  facet_wrap( ~ Cluster_Label, nrow = 4, ncol = 2, labeller = cluster_labeller_atac) +
  gplot.theme + theme(legend.position = "none")

```

# Basic cluster annotation.

```{r basic_hg19_annotation, message=FALSE, warning=FALSE}

hg19.basic.genes <- c('hg19_basicgenes', 'hg19_genes_intergenic')

hg19.basic.genes.annotation <- build_annotations(genome = "hg19", annotations = hg19.basic.genes)

high.affinity.peaks.annotated.granges <-
  makeGRangesFromDataFrame(high.affinity.peaks.annotated,
  keep.extra.columns = T)
  
high.affinity.peaks.annotated.granges <-
  annotate_regions(
  regions = high.affinity.peaks.annotated.granges,
  annotations = hg19.basic.genes.annotation,
  ignore.strand = T,
  quiet = FALSE
  )

```

## Basic annotation table.

```{r basic_annotation_table}

# Annotation type data frame

high.affinity.peaks.annotated.granges.df <-
  data.frame(high.affinity.peaks.annotated.granges) %>%
  dplyr::select(PeakID, Cluster_Label, annot.type) %>%
  dplyr::distinct(PeakID, annot.type, .keep_all = T)#


basic.annotation.type.freqs <-
  high.affinity.peaks.annotated.granges.df %>%
  dplyr::group_by(Cluster_Label, annot.type) %>%
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(total = sum(n), rel.freq = n / total)


basic.annotation.type.freqs$Cluster_Label <- factor(x = basic.annotation.type.freqs$Cluster_Label,
                                                    levels = paste("C", 1:8, sep = " "))
knitr::kable(
  basic.annotation.type.freqs %>% dplyr::select(Cluster_Label, annot.type, rel.freq) %>% arrange(Cluster_Label)
  )

```

## Basic annotation plot.

```{r basic_annotation_plot}

ggplot(data = basic.annotation.type.freqs, aes(x = Cluster_Label, y = rel.freq, fill = annot.type)) +
  geom_bar(position = "stack", stat = "identity") +
  coord_flip() +
  scale_x_discrete(limits = rev(levels(
  basic.annotation.type.freqs$Cluster_Label
  ))) +
  scale_fill_brewer(palette = "Paired") +
  xlab("ATAC-seq clusters") +
  ylab("Proportion of Cluster specific ATAC-seq peaks class") +
  labs(fill = "Class") +
  gplot.theme + theme(
  legend.position = "right",
  legend.title = element_text(size = 10),
  legend.text = element_text(size = 9),
  axis.title = element_text(size = 10),
  axis.text = element_text(size = 9)
  )
```


# Integration of chromatin accessibility and gene expression profiles.
We show global gene expression profiles of genes associated with peaks in the clusters identified above.
Also, we show profiles of selected genes associated with primed and naive reprogramming and somatic regulatory elements. 

Based on detailed cluster annotation (clusters annotated with Homer, see methods) and peaks asociated with gene transcription starting sites TTSs (see methods). In cases were a gene has more than one associated peak, the closest peak to the TSS is selected.

For simplicity, objects with the annotated clusters and z-scores of expressed genes are loaded (`data/bkATAC_annotated_clusters.Rdata` and `data/bkATAC_gene_expression_bulk_RNA_seq.Rdata`).

## Load annotated clusters and process.
```{r annotated_clusters}
# Load annotated cluster and subset peaks closest to TSSs

load(file = "data/bkATAC_annotated_clusters.Rdata") # Creates `annotated.clusters` object.

# Filter out peaks with no Ensembl ID associated.

annotated.clusters.table <- ldply(annotated.clusters, rbind) %>% filter(!Nearest.Ensembl == "" ) %>% .[, -1] %>% dplyr::select(-c("Peak.Score", "Focus.Ratio.Region.Size"))

# In cases were a gene has more than one associated peak, the closest peak to the TSS is selected.
annotated.clusters.table.nearest.TSS <-
  as.data.frame(annotated.clusters.table) %>%
  dplyr::group_by(Nearest.Ensembl) %>%
  dplyr::slice(which.min(abs(Distance.to.TSS)))

# Add cluster information and chromatin accessiblity z-scores.  
annotated.clusters.table.nearest.TSS <-
  dplyr::left_join(
  annotated.clusters.table.nearest.TSS,
  dplyr::select(
  high.affinity.peaks.annotated,
  -c("Chr", "Start", "End", "Strand", "Length")
  ),
  by = "PeakID"
  )

```

## Genes associated to ATAC-seq peak clusters.

### Expression table.

```{r z_scores_rna_seq}
base::load(file = "data/bkATAC_gene_z_scores_bkRNA_seq.Rdata") # Load RNA-seq gene z-scores object name :`bkRNA.seq.z.scores`

# Simplify TSS table to include only cluster information and gene ENSEMBL ids.

table.nearest.TSS.cluster <- dplyr::select(annotated.clusters.table.nearest.TSS, c("Nearest.Ensembl", "Cluster", "Cluster_Label"))

colnames(table.nearest.TSS.cluster)[1] <- "Ensembl_id"

# Add cluster information to expressed genes. Please note that not all expressed genes have a cluster associated to.

bkRNA.seq.z.scores.annotated <-
  dplyr::inner_join(bkRNA.seq.z.scores,
  table.nearest.TSS.cluster,
  by = "Ensembl_id")
  
```

## Table of number of expressed genes associated to ATAC-seq peak clusters.

```{r table_peak-associated_genes}

bkRNA.seq.z.scores.cluster.counts <-
  bkRNA.seq.z.scores.annotated %>%
  group_by(Cluster_Label) %>%
  dplyr::count()

knitr::kable(bkRNA.seq.z.scores.cluster.counts)
```

## Expression profile of genes associated to ATAC-seq peak clusters.

```{r expression_genes_cluster, warning=FALSE, fig.align="center", fig.height=9, fig.width=4.5}

bkRNA.seq.z.scores.annotated.melt <-
  melt(
  bkRNA.seq.z.scores.annotated %>% dplyr::select(-c("geneName", "Cluster")),
  id = c("Ensembl_id", "Cluster_Label")
  )
  
bkRNA.seq.z.scores.annotated.melt <-
  bkRNA.seq.z.scores.annotated.melt %>%
  mutate(media = ifelse(grepl("Fibroblast", variable), "FM",
                        ifelse(grepl("Primed", variable), "PM", "NM")))  %>%
  mutate(stage = gsub(".*_", "", variable))

bkRNA.seq.z.scores.annotated.melt$stage <-
  factor(
  bkRNA.seq.z.scores.annotated.melt$stage,
  levels = c("D00", "D03", "D07", "D13", "D21", "P03", "P10"),
  labels = c("D0", "D3", "D7", "D13", "D21", "P3", "P10")
  )

# Labeller defined as above.  
cluster.labels.rna.seq <- list("C 1" = "C 1 (n = 1823)", 
                                 "C 2" = "C 2 (n = 1347)",
                                 "C 3" = "C 3 (n = 892)",
                                 "C 4" = "C 4 (n = 339)",
                                 "C 5" = "C 5 (n = 1239)",
                                 "C 6" = "C 6 (n = 1027)", 
                                 "C 7" = "C 7 (n = 790)",
                                 "C 8" = "C 8 (n = 1547)")

cluster_labeller_rna <- function(variable,value){
  return(cluster.labels.rna.seq[value])
}

ggplot(
  data = bkRNA.seq.z.scores.annotated.melt %>%
  dplyr::filter(media %in% c("FM", "PM")),
  aes(stage, value, group = Ensembl_id, colour = media)
  ) +
  stat_summary(aes(group = 1),
  fun.y = mean,
  geom = "line",
  size = 1) +
  stat_summary(
  data = bkRNA.seq.z.scores.annotated.melt %>%
  dplyr::filter(media %in% c("FM", "NM")),
  aes(group = 1),
  fun.y = mean,
  geom = "line",
  size = 1
  ) +
  scale_color_manual(values = c("FM" = "black", "NM" = "#0000FF", "PM" = "#FF8C00"),
  name = "Culture Media") +
  facet_wrap( ~ Cluster_Label,
  nrow = 4,
  ncol = 2,
  labeller = cluster_labeller_rna) +
  ylab("Gene expression (z-score)") +
  xlab("Reprogramming stages") +
  gplot.theme

```

## Selected genes chromatin accessibility and expression profiles.
The following genes were selected: _SNAI1_, _TWIST2_, _ZEB2_ , _BHLHE40_, _F11R_, _ZIC5_, _OCLN_, _NLRP7_.

### Chromatin accessiblity profiles of selected genes.

```{r chromatin_accessibility_selected_genes}

selected.genes <- data.frame(geneName = c("SNAI1", "TWIST2", "ZEB2", "BHLHE40", "F11R", "ZIC5", "OCLN", "NLRP7"))

selected.genes <- dplyr::inner_join(selected.genes,
dplyr::select(bkRNA.seq.z.scores.annotated,
c("Ensembl_id", "geneName")),
by = "geneName")

# Chromatin accessibility

chromatin.accessibiliy.selected.genes <-
  annotated.clusters.table.nearest.TSS %>% ungroup() %>% dplyr::select(
  "Gene.Name",
  "D0",
  "D3",
  "D7",
  "D13_Primed",
  "D21_Primed",
  "P3_Primed",
  "P10_Primed",
  "D13_Smith",
  "D21_Smith",
  "P3_Smith",
  "P10_Smith"
  ) %>% filter(Gene.Name %in% selected.genes$geneName)
  
chromatin.accessibiliy.selected.genes.melt <-
  melt(chromatin.accessibiliy.selected.genes, id = c("Gene.Name"))

chromatin.accessibiliy.selected.genes.melt <-
  chromatin.accessibiliy.selected.genes.melt %>%
  mutate(stage = case_when(
  grepl("D13", variable) ~ "D13",
  grepl("D21", variable) ~ "D21",
  grepl("P3", variable) ~ "P3",
  grepl("P10", variable) ~ "P10",
  TRUE ~ as.character(variable)
  )) %>%
  mutate(media = case_when(
  grepl("D0|D3|D7", variable) ~ "FM",
  grepl("Smith", variable) ~ "NM",
  grepl("Primed", variable) ~ "PM"
  ))
  

chromatin.accessibiliy.selected.genes.melt$stage <-
  factor(
  chromatin.accessibiliy.selected.genes.melt$stage,
  levels = c("D0", "D3", "D7", "D13", "D21", "P3", "P10")
  )


chromatin.accessibiliy.selected.genes.melt$Gene.Name <-
  factor(
  chromatin.accessibiliy.selected.genes.melt$Gene.Name,
  levels = c(
  "SNAI1",
  "TWIST2",
  "ZEB2",
  "BHLHE40",
  "F11R",
  "ZIC5",
  "OCLN",
  "NLRP7"
  )
  )
  

plot.chromatin.accessibiliy.selected.genes <- ggplot(
  data = chromatin.accessibiliy.selected.genes.melt %>%
  filter(
  media %in% c("FM", "PM")),
  aes(stage, value, group = Gene.Name, colour = media)
  ) +
  geom_line(size = 1.5) +
  geom_line(
  data = chromatin.accessibiliy.selected.genes.melt %>%
  filter(
  media %in% c("FM", "NM")
  ),
  aes(stage, value, group = Gene.Name, colour = media), size = 1.5
  ) +
  scale_color_manual(values = c("FM" = "black", "NM" = "#0000FF", "PM" = "#FF8C00"),
  name = "Culture Media") +
  facet_wrap( ~ Gene.Name, nrow = 8, ncol = 1) +
  ylab("Chromatin Accessiblity (Z-score)") +
  xlab("Reprogramming stages") +
  gplot.theme

```

### Gene expression profiles of selected genes.

```{r gene_expression_profile, warning=FALSE}

gene.expression.selected.genes.melt <-
  dplyr::inner_join(selected.genes,
    bkRNA.seq.z.scores.annotated.melt,
    by = "Ensembl_id")

gene.expression.selected.genes.melt$geneName <-
  factor(
  gene.expression.selected.genes.melt$geneName,
  levels = c(
  "SNAI1",
  "TWIST2",
  "ZEB2",
  "BHLHE40",
  "F11R",
  "ZIC5",
  "OCLN",
  "NLRP7"
  )
  )

gene.expression.selected.genes.plots <-
  ggplot(data = gene.expression.selected.genes.melt %>%
  filter(media %in% c("FM", "PM")),
  aes(stage, value, group = geneName, colour = media)) +
  geom_line(size = 1.5) +
  geom_line(
  data = gene.expression.selected.genes.melt %>%
  filter(media %in% c("FM", "NM")),
  aes(stage, value, group = geneName, colour = media),
  size = 1.5
  ) +
  scale_color_manual(values = c("FM" = "black", "NM" = "#0000FF", "PM" = "#FF8C00"),
  name = "Culture Media") +
  facet_wrap( ~ geneName, nrow = 8, ncol = 1) +
  ylab("Gene expression (Z-score)") +
  xlab("Reprogramming stages") +
  gplot.theme

```

```{r accessibliyt_expression_profiles_selected_genes, fig.height=12.5, fig.width=7.5, echo = FALSE, warning=FALSE}

plot.chromatin.accessibiliy.selected.genes | gene.expression.selected.genes.plots

```

# Motif enrichment analysis.

## Top 25 most enriched motifs accross clusters (based on p-value).
For simplicity, objects with enriched motifs accross clusters are loaded (`data/bkATAC_motif_enrichment.Rdata` and `data/bkATAC_gata2.ap2.gamma.motif.enrichment.Rdata`)

In circumstances where motifs have identical symbols, the motif with the highest -(log10(p-value)) was kept; e.g, cases such as EWS:ERG-Fusion and ERG. 

```{r motif_enrichment, fig.height=5, fig.width=12.5}

 # Load `bkATAC_motif_enrichment` object containing information of top 25 most enriched motifs accross clusters.
base::load("data/bkATAC_motif_enrichment.Rdata")

# Add cluster labels as above

bkATAC_motif_enrichment <- bkATAC_motif_enrichment %>% dplyr::mutate(
  Cluster_Label = case_when(
  cluster == "Cluster 1" ~ "C 1",
  cluster == "Cluster 2" ~ "C 2",
  cluster == "Cluster 3" ~ "C 7",
  cluster == "Cluster 4" ~ "C 4",
  cluster == "Cluster 5" ~ "C 3",
  cluster == "Cluster 6" ~ "C 6",
  cluster == "Cluster 7" ~ "C 5",
  TRUE ~ "C 8"
  )
  )

bkATAC_motif_enrichment  <-  data.table(bkATAC_motif_enrichment)
bkATAC_motif_enrichment$motifHS  <- bkATAC_motif_enrichment$human_symbol

dup.TF <-  unique(bkATAC_motif_enrichment[, c("motif.hs", "motifHS")])
dup.TF <-  dup.TF$motifHS[duplicated(dup.TF$motifHS)]
oup.Motif = bkATAC_motif_enrichment[!(motifHS %in% dup.TF)]

for(i in dup.TF){
  tmp = bkATAC_motif_enrichment[motifHS == i]
  tmp = tmp[tmp[, .I[which.max(m.Log.P.Value)], by = c("cluster", "motifHS")]$V1]
  oup.Motif = rbindlist(list(oup.Motif, tmp))
}

oup.Motif$mLogP = oup.Motif$m.Log.P.Value
oup.Motif$target = oup.Motif$target_sequence_percentage * 100

oup.Motif$clustPlot = factor(oup.Motif$Cluster_Label, levels = rev(paste("C", 1:8, sep = " ")))

oup.Motif$motifHS = factor(oup.Motif$motifHS,
                           levels = unique(oup.Motif[order(-clustPlot,-mLogP)]$motifHS))

gg.Data = oup.Motif

gg.Data[target > 50]$target = 50

ggplot(gg.Data, aes(motifHS, clustPlot, color = mLogP, size = target)) +
  geom_point() +
  xlab("Transcription Factor") + ylab("ATAC Cluster") +
  scale_size_continuous(
  "Target sequence with motif (%)",
  range = c(0, 5.5),
  limits = c(0, 50),
  breaks = c(0, 25, 50),
  labels = c("0", "25", ">50")
  ) +
  scale_color_gradientn("-log (Enrichment P-value)",
  limits = c(50, 6000),
  breaks = c(0, 2000, 4000, 6000),
  colours = c("darkblue", "darkred")
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.text = element_text(size = 8), 
        legend.position = "bottom", legend.spacing.y = unit(0, "cm"),
        panel.grid.minor = element_blank())

```

## TFAP2C and GATA2 motif enrichment.

```{r tfap2c_motif_enrichment}
# Load object with motif enrichment data of GATA2 and TFAP2C (AP2-gamma).

base::load("data/bkATAC_gata2.ap2.gamma.motif.enrichment.Rdata")

tfap2c.motif.enrichment.plot <-
  ggplot(data = gata2.ap2.gamma.motif.enrichment %>% dplyr::filter(motif.hs == "TFAP2C"),
  aes(x = Cluster_Label,
  y = -log_P_value)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  xlab("ATAC-seq clusters") +
  ylab("-Log(Enrichment p-value)") +
  ggtitle("TFAP2C motif enrichment") +
  theme_classic()
  
gata2.motif.enrichment.plot <-
  ggplot(data = gata2.ap2.gamma.motif.enrichment %>% dplyr::filter(motif.hs == "GATA2"),
  aes(x = Cluster_Label,
  y = -log_P_value)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  xlab("ATAC-seq clusters") +
  ylab("-Log(Enrichment p-value)") +
  ggtitle("GATA2 motif enrichment") +
  theme_classic()

tfap2c.motif.enrichment.plot | gata2.motif.enrichment.plot
```

# Session information
### R 
```{r session_info}
sessionInfo()
```