---
title: "3-deazaadenosine (3DA) alleviates senescence to promote cellular fitness and therapy efficiency"

date: ""
abstract: ""
output:
   BiocStyle::html_document:
    theme: cosmo
    highlight:  "tango"
    toc: true
    toc_float:
      toc_collapsed: true
    toc_depth: 2
    number_sections: true
    code_folding: "show"
    df_print: kable

vignette: >
    %\VignetteIndexEntry{Vignette Title}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8} 
---

\pagebreak
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(cache=TRUE, eval = TRUE, results = "hide", message = FALSE, warning = FALSE, echo = TRUE, fig.wide = TRUE,  fig.align='center')
```

# H3K36me3 ChIP-seq data

Briefly, H3K36me3 ChIP-seq data were processed as follow : reads were cleaned and trimmed using fastq-mcf from the ea-utils suite v1.1.2 to remove adapters, low quality bases and reads, and discard reads shorter than 25 bp after filtering. Reads were then aligned to the human reference genome (hg19) with bowtie v1.1.1 using best matches parameters (bowtie -v 2 -m 1 --best --strata). Alignment files were further processed with samtools v1.2 and PicardTools v1.130 to flag PCR and optical duplicates and remove alignments located in Encode blacklisted regions. Fragment size was estimated *in silico* for each library using spp v1.10.1. Enriched regions were identified for each replicate independently with MACS v2.1.0 with non-IPed genomic DNA as a control (macs2 callpeak --nomodel --shiftsize --shift-control --gsize hs -p 1e-1). These relaxed peak lists were then processed through the irreproducible discovery rate (IDR) pipeline to generate an optimal and reproducible set of peaks for each histone modification type and each time point.

## Preprocessing of the data

### Import peaks

Using `r Biocpkg("DiffBind")`, we first generate a consensus set of peaks, merging all the peaks defined with the IDR for 3DA, 4OHT and DMSO. 

```{r}
library("DiffBind")
library("DESeq2")
library("biomaRt")
library("GenomicFeatures")
library("org.Hs.eg.db")

#--> Read experimental design
setwd("~/Documents/Post-Doc/B-Projects/JESUS_K36ME3/")
path_to_samples <- "./data/ChIP-seq/Plan_Expe_ChIP_DiffBind_"
sample_sheet <-  read.csv(paste(path_to_samples, "K36ME3", ".csv", sep = ""), sep = ";")

#--> Create the global DBA object
expe <- dba(sampleSheet = sample_sheet, skipLines = 0, bRemoveM=TRUE)
expe_consensus <- dba.peakset(expe, consensus = -DBA_REPLICATE)
expe_ol <- dba.overlap(expe_consensus,  expe_consensus$masks$Consensus)
```

### Explore commonalities and specificities

```{r, fig.cap = "Comparison of H3K36me3 peaks between DMSO, 4OHT and 3DA treated cells", fig.height=8, fig.width=8, fig.align='center', results = "show"}
#--> Get the data
DMSO <- dba.peakset(expe, peaks = expe$masks$DMSO, consensus=-DBA_REPLICATE, bRetrieve = TRUE)
DA3 <- dba.peakset(expe, peaks = expe$masks$`3DA`, consensus=-DBA_REPLICATE, bRetrieve = TRUE)
OHT4 <- dba.peakset(expe, peaks = expe$masks$`4OHT`, consensus=-DBA_REPLICATE, bRetrieve = TRUE)

#--> Load libraries
library(ChIPpeakAnno)
library(eulerr)
library(colorRamps)

#--> Format data
ol <- findOverlapsOfPeaks(DMSO, DA3, OHT4)
venn_data <- ol$venn_cnt[,"Counts"][-1]
names(venn_data) <- c("OHT4", "DA3", "DA3&OHT4", "DMSO", "DMSO&OHT4", "DMSO&DA3", "DMSO&DA3&OHT4")
venn_p <- euler(venn_data)

#--> Number of common / specific peaks
print(ol$venn_cnt)

#--> Draw Euler plot
plot(venn_p)
```

As depicted on the Euler above, while 19531 reproducible peaks are common to the 3 conditions, 2024 are specific to DMSO, 2370 to 3DA, and 5693 to 4OHT.

## TSS centric analysis 

Before starting this analysis, we need to retrieve the coordinates of TSS on a data-base such as UCSC or Ensembl. Here, we use the Ensembl database GRCh37 corresponding to the genome build *hg19*. We next consider the window +/- 10kb on each side of the TSS and split it in 1000 non-overlapping windows.

### Preparing the data 

```{r}
#--> Load the biomeRt and the biovizBase packages
library(biomaRt)
library(biovizBase)

#--> Connexion to BioMart hg19 (aka GRCh37)
ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL",
                   host="grch37.ensembl.org",
                   path="/biomart/martservice",
                   dataset="hsapiens_gene_ensembl")

#--> Retrive the coordinates of TSSs
tss <- getBM(mart = ensembl,
             attributes = c("chromosome_name", "transcript_start", "transcript_end", "hgnc_symbol", "strand"))

#--> Format the chromosome name column (to match the chromosome names in alignement files)
set.seed(1234)
tss <- tss[grep("^[1-9]{1,2}$", tss$chromosome_name),]

#--> Creating a GRange containing the information related to TSSs
# And add a +/- 1kb window
tss_gr <- GRanges(seqnames = paste0("chr", tss$chromosome_name),
                  ranges = IRanges(start = tss$transcript_start - 10000,
                                   end = tss$transcript_end + 10000,
                                   names = tss$hgnc_symbol),
                  strand = tss$strand)


#--> For each TSS +/- 5kb, tile the window in 100bp bins
tss_gr_plus  <- tss_gr[which(strand(tss_gr) == "+")]
tss_gr_plus_split <- tile(tss_gr_plus, n = 1000)
tss_gr_plus_split <- flatGrl(tss_gr_plus_split)
tss_gr_minus  <- tss_gr[which(strand(tss_gr) == "-")]
tss_gr_minus_split <- tile(tss_gr_minus, n = 1000)
tss_gr_minus_split <- flatGrl(tss_gr_minus_split)

```

Now that the object containing information related to TSSs is formated, we will first list the BigWig files we are interested in, and count the reads overlapping each tiles consituting each TSS window.

### Compute coverages across TSS

```{r}
library(edgeR)
library(ggplot2)
library(IRanges)

#--> List input averages and normalized BigWig files
bw <- list.files(path = "./data/ChIP-seq/", pattern = "*_ALL_DEDUP_NOBLACKLIST.bw", full = TRUE)[c(1:3)]

#--> Count signal in bins
counts_plus  <- matrix(NA, nrow = length(tss_gr_plus_split), ncol = 3)
counts_minus <- matrix(NA, nrow = length(tss_gr_minus_split), ncol = 3)

chrs <- c(paste("chr", seq(1, 22), sep = ""))

for(i in c(1:length(bw))) {
  
  print(i)
  
  coverage <- rtracklayer::import(bw[i], as = 'RleList')

    for(chr in chrs) {

      counts_plus[which(seqnames(tss_gr_plus_split) == chr) , i] <-
        mean(Views(coverage[[which(names(coverage) == chr)]], 
                   ranges(tss_gr_plus_split[seqnames(tss_gr_plus_split) == chr])))
        
        counts_minus[which(seqnames(tss_gr_minus_split) == chr) , i] <-
        mean(Views(coverage[[which(names(coverage) == chr)]],
                   ranges(tss_gr_minus_split[seqnames(tss_gr_minus_split) == chr])))
        
    }
}

#--> Norm factors
norm <- c(51217909, 50988684, 50578426)

counts_plus  <- cpm(counts_plus, lib.size = norm)
counts_minus <- cpm(counts_minus, lib.size = norm)


#--> Format outputs
count_3DA_R1  <- rbind(t(matrix(counts_plus[,1], nrow = 1000)),
                    t(matrix(counts_minus[,1], nrow = 1000))[,rev(1:1000)])
count_4OHT_R1  <- rbind(t(matrix(counts_plus[,2], nrow = 1000)),
                    t(matrix(counts_minus[,2], nrow = 1000))[,rev(1:1000)])    
count_DMSO_R1  <-  rbind(t(matrix(counts_plus[,3], nrow = 1000)),
                       t(matrix(counts_minus[,3], nrow = 1000))[,rev(1:1000)])  
```

### Build the metaprofiles 

```{r, cache.lazy = FALSE, fig.cap= "H3K36me3 metaprofiles around TSS", fig.height=7}

#--> Format outputs
count_all_1 <- list(count_3DA_R1, count_4OHT_R1, count_DMSO_R1)

count_all_2 <- lapply(count_all_1, function(x) colMeans(x, na.rm = TRUE))

count_all_3 <- data.frame(Pos = c(rep(rep(seq(from = -10000, to = 9999, by = 20),1),3)),
                          Condition = rep(c("DA3_1", "OHT4_1", "DMSO_1"), each = 1000),
                          Reads = unlist(count_all_2))

#--> Plot profiles
ggplot(count_all_3, aes(x = Pos, y = Reads)) +
  geom_line(aes(color = Condition)) +
  scale_x_continuous(name = "Distance to TSS (bp)") +
  scale_y_continuous(name = "Reads") +
  theme_bw()

```

## Generate the raw count matrix

To perform a differential analysis, we finally combined all those peaks in a consenus global set and computed the raw read coverage for each samples / conditions using `r Biocpkg("DiffBind")`. During the process, all H3K36me3 ChIP-seq peaks seperated by gaps smaller than 10kb were merged. 

```{r}
#--> Create consensus peakset
consensus <- dba.peakset(expe, consensus = DBA_FACTOR)
peaks <- dba.peakset(consensus, consensus$masks$ALL, bRetrieve=TRUE)

#--> Count binding in consensus peaks 
count <- dba.count(expe,
                   peaks = reduce(peaks, min.gapwidth = 10000),
                   score = DBA_SCORE_TMM_MINUS_EFFECTIVE,
                   bLog = FALSE,
                   bParallel = TRUE)

#--> Outputing the binding matrix containing filtered row counts
b_matrix <- dba.peakset(count, bRetrieve=TRUE, DataType=DBA_DATA_FRAME)
```

## Scaling and normalization

To consider differences in sequencing depth and immuno-precipitation efficiency, we applied :

* a quantile normalization using `r Biocpkg("EDASeq")`.
* a debatching strategy identifying surrogate variables using `r Biocpkg("RUVSeq")`.

```{r}
#--> Prepare the data
library(EDASeq)
library(RUVSeq)
library(RColorBrewer)
data_RUVs <- matrix(as.numeric(unlist(round(b_matrix[,-c(1:3)]))), nrow = nrow(b_matrix[,-c(1:3)]))

#-->  Quantile normalization
set <- newSeqExpressionSet(counts = data_RUVs)
colnames(set) <- sample_sheet$SampleID
qn_counts <- betweenLaneNormalization(set, which = "full")

#-->  Debatching
differences <- matrix(data = c(1,2,3,4,5, 6), ncol = 2, byrow = TRUE)
set_RUVs <- RUVs(x = qn_counts, cIdx = row.names(set), k = 1, scIdx = differences)

```

## Graphical QC

As depicted on the PCA plot below, when looking at the raw count matrix, the PC1 (71.5% of the variability) captures the variability introduce by the batch-of-origin of samples while the PC2 (17.2% of variability) captures the experimental variability. This suggest that the experiment is impaired by a non-neglectable batch effect (PC1 >> PC2) that we need to take into account.

Nevertheless, we can already notice that the biggest effect is driven by senescence induction (distance between DMSO and 4OHT on the PC2) while the effect of 3DA is rather mild (distance between 3DA and 4OHT on PC2). Moreover, the use of 3DA tend to push samples downward on the PC2 suggesting that the H3K36me3 profile of 3DA treated cells tend to be bend toward the one of DMSO growing cells.

```{r, fig.cap = "PCA on raw counts - H3K36me3", fig.height=10}
#--> PCA on raw data
colors <- brewer.pal(6, "Set1")
x <- as.factor(sample_sheet$Condition)
PCA_Raw <- plotPCA(set, col=colors[x], cex=1.2,
                   main = "PCA on raw counts")
```

Despite the quantile normalization, PC1 still captures a huge proportion of variability arising from the batch-of-origin. Further processing is required. 

```{r, fig.cap = "PCA on quantile-normalized pseucounts - H3K36me3", fig.height=10}
#--> PCA on quantile normalized data
PCA_UQ <- plotPCA(qn_counts, col=colors[x], cex=1.2,
                  main = "PCA on quantile-normalized pseucounts")
```

After non-linear debatching, we managed to get rid of the batch effect. Our initial conlusions are still valid in this sub-space : distance on PC1 between DMSO and 4OHT samples is much bigger than the distance between 3DA and 4OHT cells.

```{r, fig.cap = "PCA on quantile-normalized + RUVSeq pseudocounts -  H3K36me3", fig.height=10}
#--> PCA on quantile normalized and debatched data
PCA_RUVs <- plotPCA(set_RUVs, col=colors[x], cex=1.2,
                    main = "PCA on quantile-normalized + RUVSeq pseudocounts")
```

## Differential analysis

For the diffenrential analysis, we used `r Biocpkg("edgeR")` along with a linear model allowing to test for differences in H3K36me3 between DMSO and 4OHT (effect of senescence induction), and between 3DA and 4OHT (effect of the treatment on senescent cells).

The peaks were finally annotated using `r Biocpkg("bimoaRt")` and `r Biocpkg("ChIPpeakAnno")` considering the hg19 assembly. Here, a peaks was flag with the name of a gene when included in the gene body of the gene (either TSS, TTS, exons, introns).

```{r}
library(edgeR)

#--> Prepare the data
data_edger_chip <- data_RUVs
samples <- data.frame(sample_sheet$Condition)
design <- cbind(model.matrix(~ 0  + samples[,1]))
colnames(design) <- c("DA", "OHT", "DMSO")

#--> Build the DGEList object
data_edger_chip  <- DGEList(counts = data_edger_chip,
                            group = samples[,1],
                            genes = paste(b_matrix[,1],
                                          b_matrix[,2],
                                          b_matrix[,3], sep= "_"))

#--> Estimate dispersion and fit the model
data_edger_chip <- calcNormFactors(data_edger_chip, method = "none", lib.size = norm)
data_edger_chip <- estimateGLMRobustDisp(data_edger_chip, design)
fit_chip <- glmFit(data_edger_chip, design = design)

#--> DE peaks between 4OHT and DMSO
contrasts_4OHT_chip <- makeContrasts(Senescence =  OHT - DMSO, levels = design)
lrt_4OHT_chip <- glmLRT(fit_chip, contrast = contrasts_4OHT_chip)
top_4OHT_chip <- topTags(lrt_4OHT_chip, n = Inf)$table
colnames(top_4OHT_chip) <- paste("edgeR_ChIP_OHT",
                                 colnames(top_4OHT_chip), sep = ".")

#--> DE peaks between 4OHT and 3DA
contrasts_3DA_chip <- makeContrasts(Senescence =   DA - OHT, levels = design)
lrt_3DA_chip <- glmLRT(fit_chip, contrast = contrasts_3DA_chip)
top_3DA_chip <- topTags(lrt_3DA_chip, n = Inf)$table
colnames(top_3DA_chip) <- paste("edgeR_ChIP_3DA",
                               colnames(top_3DA_chip), sep = ".")

#--> Output
gr_edger_chip <- merge(top_4OHT_chip, top_3DA_chip, by = 1)
gr_edger_chip$chr <- strsplit2(gr_edger_chip[,1],"_")[,1]
gr_edger_chip$start <- strsplit2(gr_edger_chip[,1],"_")[,2]
gr_edger_chip$end <- strsplit2(gr_edger_chip[,1],"_")[,3]
gr_edger_chip <- GRanges(gr_edger_chip)

#--> Annotation
library(ChIPpeakAnno)
library(org.Hs.eg.db)

ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL",
                   host="grch37.ensembl.org",
                   path="/biomart/martservice",
                   dataset="hsapiens_gene_ensembl")

ensembl_gr <- getAnnotation(mart = ensembl)
gr_edger_chip <- annotatePeakInBatch(gr_edger_chip, AnnotationData = ensembl_gr,
                                     output = "overlapping")
gr_edger_chip <- addGeneIDs(gr_edger_chip, org.Hs.eg.db, IDs2Add=c("symbol", "entrez_id"), 
           feature_id_type="ensembl_gene_id", silence=TRUE, mart)

```

## Visual inspection

```{r, fig.cap="Differential analysis for H3K36me3 ChIP-seq : MA-plots", fig.height=15}
library(ggplot2)
library(ggrepel)
library(clusterProfiler)

#--> Import the senescence signatures
gmt_senescence <- qusage::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/SenescenceSignature.gmt")

#--> Build MA plot for the senescence induction
MA_OHT <- ggplot() +
         geom_point(data = data.frame(gr_edger_chip),
                    aes(y = edgeR_ChIP_OHT.logFC,
                        x = edgeR_ChIP_OHT.logCPM)) +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  edgeR_ChIP_OHT.logFC < -log2(1.2) & 
                                  edgeR_ChIP_OHT.FDR < .1),
                    aes(y = edgeR_ChIP_OHT.logFC,
                        x = edgeR_ChIP_OHT.logCPM),
                    color = "firebrick1") +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  edgeR_ChIP_OHT.logFC > log2(1.2) &
                                  edgeR_ChIP_OHT.FDR < .1),
                    aes(y = edgeR_ChIP_OHT.logFC,
                        x = edgeR_ChIP_OHT.logCPM),
                    color = "chartreuse3") +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$SASP)),
                    aes(y = edgeR_ChIP_OHT.logFC,
                        x = edgeR_ChIP_OHT.logCPM),
                    color = "steelblue") +
         geom_label_repel(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$SASP) &
                                  edgeR_ChIP_OHT.logFC > 1 &
                                  edgeR_ChIP_OHT.FDR < .1),
               aes(y = edgeR_ChIP_OHT.logFC,
                   x = edgeR_ChIP_OHT.logCPM,
                   label = symbol)) +
         scale_x_continuous(name = "logCPM") +
         scale_y_continuous(name = "logFC 4OHT over DMSO") +
         ggtitle("Senescence induction : Impact on H3K36me3") +
         theme_bw()

#--> Build MA plot for 3DA treatment
MA_DA <- ggplot() +
         geom_point(data = data.frame(gr_edger_chip),
                    aes(y = edgeR_ChIP_3DA.logFC,
                        x = edgeR_ChIP_3DA.logCPM)) +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  edgeR_ChIP_3DA.logFC < -log2(1.2) & 
                                  edgeR_ChIP_3DA.FDR < .1),
                    aes(y = edgeR_ChIP_3DA.logFC,
                        x = edgeR_ChIP_3DA.logCPM),
                    color = "firebrick1") +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  edgeR_ChIP_3DA.logFC > log2(1.2) &
                                  edgeR_ChIP_3DA.FDR < .1),
                    aes(y = edgeR_ChIP_3DA.logFC,
                        x = edgeR_ChIP_3DA.logCPM),
                    color = "chartreuse3") +
         geom_point(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$SASP)),
                    aes(y = edgeR_ChIP_3DA.logFC,
                        x = edgeR_ChIP_3DA.logCPM),
                    color = "steelblue") +
         geom_label_repel(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$SASP) &
                                  abs(edgeR_ChIP_3DA.logFC) > .2 &
                                  edgeR_ChIP_3DA.FDR < .1),
               aes(y = edgeR_ChIP_3DA.logFC,
                   x = edgeR_ChIP_3DA.logCPM,
                   label = symbol)) +
         scale_x_continuous(name = "logCPM") +
         scale_y_continuous(name = "logFC 3DA over 4OHT") +
         ggtitle("3DA : Impact on H3K36me3") +
         theme_bw()

#--> Build FC vs FC plot
FC_FC <- ggplot() +
    geom_point(data = data.frame(gr_edger_chip),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_ChIP_OHT.logFC),
               color = "grey") +
    geom_smooth(method = "gam", data = data.frame(gr_edger_chip),
               aes(x = edgeR_ChIP_3DA.logFC, y = edgeR_ChIP_OHT.logFC)) +
    geom_point(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$`OIS FC1.5`)),
                    aes(x = edgeR_ChIP_3DA.logFC,
                        y = edgeR_ChIP_OHT.logFC),
                    color = "steelblue") +
     geom_label_repel(data = subset(data.frame(gr_edger_chip),
                                  is.element(data.frame(gr_edger_chip)$symbol, 
                                             gmt_senescence$`OIS FC1.5`) &
                                  edgeR_ChIP_OHT.FDR < .1 &
                                  edgeR_ChIP_3DA.FDR < .1),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_ChIP_OHT.logFC,
                   label = symbol)) +
    xlab("log2 fold change in H3K36me3 (3DA over DMSO)") +
    ylab("log2 fold change in H3K36me3 (OHT over 4OHT)")

library(gridExtra)
grid.arrange(MA_OHT, MA_DA,  nrow = 2)
    
```

On the two MA-plot above, peaks gaining H3K36me3 during senescence or after 3DA treatment are colored in green, while the one losing H3K36me3 are colored in red. The overploted blue points depict SASP genes (as defined in the list provided by Jesus Gil).

As expected, we see that most of the SASP genes gain H3K36me3 after OIS (which know to be correlated to gene expression level). Nevertheless, these genes are rather evenly distributed among regions gaining / losing H3K36me3 upon 3DA treatment. 

```{r, results = "show", fig.cap="Differential analysis for H3K36me3 ChIP-seq : FC vs FC plot", fig.height=8}
FC_FC
cor.test(gr_edger_chip$edgeR_ChIP_3DA.logFC,
         gr_edger_chip$edgeR_ChIP_OHT.logFC)
```

Comparing H3K36me3 fold-changes during senescence induction and upon 3DA treatment reveales a mild (but significant) negative correlation (R2 = -0.35) : peaks gaining H3K36me3 after OIS tend to lose H3K36me3 upon 3DA treatment. But this rule doesn't apply well to SASP-related peaks. 

## Functionnal exploration

To depict in greater details the impact of senescence induction and post-senescence 3DA treatment on H3K36me3 landscape, we performed GSEA ranking peaks according to the fold-change in H3K36me3 (either after OIS or 3DA treatment). We considered two gene-set collections : the Reactome collection provided by MSigDB, and the Senescence collection we built based on the litterature. 

### GSEA for senescence induction - REACTOME

```{r, fig.cap="Top up/down pathways for senenscence induction - REACTOME", fig.height=8, fig.width=8, fig.align='center'}
library(fgsea)
library(qusage)

# Define the target 
target <- gr_edger_chip

# Import the GMT - REACTOME
gmt <- qusage::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/c2.cp.reactome.v6.0.symbols.gmt")

# Order peaks according to the logFC of OHT over DMSO
order_da <- sort(target$edgeR_ChIP_OHT.logFC)
names(order_da) <- target$symbol[order(target$edgeR_ChIP_OHT.logFC, decreasing = FALSE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Select Top NES (+/-)
topPathways<- rbind(fgsea_da[order(NES)[c(1:5)],],
                    fgsea_da[order(NES,decreasing = TRUE)[c(1:5)],])
topPathways <- topPathways[order(topPathways$NES, decreasing = TRUE), pathway]

# Plot
plotGseaTable(gmt[topPathways], order_da, fgsea_da)

plotEnrichment(gmt[["REACTOME_CELL_CYCLE"]], order_da)
```

The senescence induction tend to impact positively (gain in H3K36me3 in 4OHT samples) peaks connected to genes involved in cell-cell interaction, in Notch signaling, while impacting negatively (loss in H3K36me3 in 4OHT samples) peaks connected to genes involved in cell cycle. 

### GSEA for senescence induction - SENESCENCE

```{r, fig.cap="Top up/down pathways for senenscence induction - SENESCENCE", fig.height=8, fig.width=8, fig.align='center'}
# Define the target 
target <- gr_edger_chip

# Import the GMT - REACTOME
gmt <- gmt_senescence

# Order peaks according to the logFC of OHT over DMSO
order_da <- sort(target$edgeR_ChIP_OHT.logFC, decreasing = TRUE)
names(order_da) <- target$symbol[order(target$edgeR_ChIP_OHT.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Plot
plotGseaTable(gmt, order_da, fgsea_da)

plotEnrichment(gmt[["SASP"]], order_da)
```

The senescence induction tend to impact positively (gain in H3K36me3 in 4OHT samples) peaks connected to genes involved in the SASP as well as genes known to to be up-regulated during OIS (OIS FC2 and OIS FC1.5). 

### GSEA for 3DA treatment - REACTOME

```{r, fig.cap="Top up/down pathways with 3DA treatment - REACTOME", fig.height=8, fig.width=8, fig.align='center'}
library(fgsea)

# Define the target 
target <- gr_edger_chip

# Import the GMT - REACTOME
gmt <- qusage::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/c2.cp.reactome.v6.0.symbols.gmt")

# Order peaks according to the logFC of 3DA over OHT
order_da <- sort(target$edgeR_ChIP_3DA.logFC, decreasing = TRUE)
names(order_da) <- target$symbol[order(target$edgeR_ChIP_3DA.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Select Top NES (+/-)
topPathways<- rbind(fgsea_da[order(NES)[c(1:5)],],
                    fgsea_da[order(NES,decreasing = TRUE)[c(1:5)],])
topPathways <- topPathways[order(topPathways$NES, decreasing = TRUE), pathway]

# Plot
plotGseaTable(gmt[topPathways], order_da, fgsea_da)

plotEnrichment(gmt[["REACTOME_CELL_CYCLE"]], order_da)
```

The 3DA treatment on senescent cells tend to impact positively (gain in H3K36me3 in 3DA samples) peaks connected to genes involved in the cell cycle, while impacting negatively (loss in H3K36me3 in 3DA samples) peaks connected to genes involved in RAS and ERK pathways. 

### GSEA for 3DA treatment - SENESCENCE

```{r, fig.cap="Top up/down pathways with 3DA treatment - SENESCENCE", fig.height=8, fig.width=8, fig.align='center'}
# Define the target 
target <- gr_edger_chip

# Import the GMT - REACTOME
gmt <- gmt_senescence

# Order peaks according to the logFC of 3DA over OHT
order_da <- sort(target$edgeR_ChIP_3DA.logFC, decreasing = TRUE)
names(order_da) <- target$symbol[order(target$edgeR_ChIP_3DA.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Plot
plotGseaTable(gmt, order_da, fgsea_da)

plotEnrichment(gmt[["SASP"]], order_da)
```

On the other hand, the 3DA treatment on senescent cells tend to also impact positively (gain in H3K36me3 in 3DA samples) peaks connected to genes involved in SASP, while impacting negatively (loss in H3K36me3 in 3DA samples) peaks connected to genes known to to be up-regulated during OIS (OIS FC2 and OIS FC1.5).

# RNA-seq data

RNA-seq data was processed with STAR using standard parameters.

## Generate count matrix

Using `r Biocpkg("Rsubread")`, we first generate a gene-centric raw count matrix, considering a h19-formated annotation. 

```{r}
library(Rsubread)

#-->
pdata_rna <- read.table("./data/RNA-seq/JESUS_RNASEQ_PHENO.txt")

#--> List files
bam_rna <- list.files("./data/RNA-seq/", full.names = TRUE, pattern = "bam$")
count_rna <- featureCounts(files = bam_rna, annot.inbuilt = "hg19", nthreads = 6)
count_rna_2 <- count_rna$counts

colnames(count_rna_2) <- row.names(pdata_rna)
se_rna <- SummarizedExperiment(assays = count_rna_2, colData = pdata_rna)
```

## Normalization

To maintain homogeneity throughout our analytical workflow, we applied the same normalization scheme here as the one used for the H3K36me3 ChIP-seq data.
But overall, the batch effect is here less pronounced. 

```{r}

#-->  Quantile normalization
set <- newSeqExpressionSet(counts = count_rna_2)
colnames(set) <- row.names(pdata_rna)
qn_counts <- betweenLaneNormalization(set, which = "full")

#-->  Debatching
differences <- matrix(data = c(1:12), ncol = 4, byrow = TRUE)
set_RUVs <- RUVs(x = qn_counts, cIdx = row.names(set), k = 1, scIdx = differences)

```


```{r, fig.height=10, fig.cap =  "PCA on raw counts - RNA-seq"}
colors <- brewer.pal(6, "Set1")
x <- as.factor(pdata_rna$Treatment)

set <- newSeqExpressionSet(counts = count_rna_2)
PCA_Raw <- plotPCA(set, col=colors[x], cex=1.2,
                   main = "PCA on raw counts")
```


```{r, fig.height=10, "PCA on quantile-normalized pseucounts - RNA-seq"}
PCA_QN <- plotPCA(qn_counts, col=colors[x], cex=1.2,
                   main = "PCA on raw counts")
```


```{r, fig.height=10, fig.cap =  "PCA on quantile-normalized + RUVSeq pseudocounts - RNA-seq"}
PCA_QN <- plotPCA(set_RUVs , col=colors[x], cex=1.2,
                   main = "PCA on raw counts")
```


## DE analysis

For the diffenrential analysis, we used `r Biocpkg("edgeR")` along with a linear model allowing to test for differences in expression between DMSO and 4OHT (effect of senescence induction), and between 3DA and 4OHT (effect of the treatment on senescent cells).

```{r}
#--> Prepare the data
data_edger_rna <-  normCounts(set_RUVs)
design <- cbind(model.matrix(~ 0  + pdata_rna$Treatment))
colnames(design) <- c("DA", "OHT", "DMSO")

#--> Build the DGEList object
data_edger_rna  <- DGEList(counts = data_edger_rna, group = pdata_rna$Treatment,
                       genes = row.names(data_edger_rna))

#-->  Filter lowly epxressed genes
filter_g <- filterByExpr(data_edger_rna, group=pdata_rna$Treatment, min.count = 1)
data_edger_rna <- data_edger_rna[filter_g,, keep.lib.sizes=FALSE]

#--> Estimate dispersion and fit the model
data_edger_rna <- calcNormFactors(data_edger_rna, method = "none")
data_edger_rna <- estimateGLMRobustDisp(data_edger_rna, design)
fit_rna <- glmFit(data_edger_rna, design = design)

#--> DE peaks between 4OHT and DMSO
contrasts_4OHT_rna <- makeContrasts(Senescence =  OHT - DMSO, levels = design)
lrt_4OHT_rna <- glmLRT(fit_rna, contrast = contrasts_4OHT_rna)
top_4OHT_rna <- topTags(lrt_4OHT_rna, n = Inf)$table
colnames(top_4OHT_rna) <- paste("edgeR_RNA_OHT", colnames(top_4OHT_rna), sep = ".")

#--> DE peaks between 4OHT and 3DA
contrasts_3DA_rna <- makeContrasts(Senescence =   DA - OHT, levels = design)
lrt_3DA_rna <- glmLRT(fit_rna, contrast = contrasts_3DA_rna)
top_3DA_rna <- topTags(lrt_3DA_rna, n = Inf)$table
colnames(top_3DA_rna) <- paste("edgeR_RNA_3DA", colnames(top_3DA_rna), sep = ".")


#--> Output
gr_edger_rna <- merge(top_4OHT_rna, top_3DA_rna, by = 1)

#--> Annotation
symbols <- AnnotationDbi::select(org.Hs.eg.db, 
                  keys = gr_edger_rna $edgeR_RNA_OHT.genes,
                  columns = c("ENTREZID", "SYMBOL"),
                  keytype = "ENTREZID")
gr_edger_rna <- merge(gr_edger_rna, symbols,
                      by.x = "edgeR_RNA_OHT.genes", by.y = "ENTREZID")  
```


## Visual inspection

```{r, fig.cap= "Differential analysis for RNA-seq : MA-plots", fig.height=15}
#--> Build MA plot for the senescence induction
MA_OHT <- ggplot() +
         geom_point(data = data.frame(gr_edger_rna),
                    aes(y = edgeR_RNA_OHT.logFC,
                        x = edgeR_RNA_OHT.logCPM)) +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  edgeR_RNA_OHT.logFC < -log2(1.2) & 
                                  edgeR_RNA_OHT.FDR < .1),
                    aes(y = edgeR_RNA_OHT.logFC,
                        x = edgeR_RNA_OHT.logCPM),
                    color = "firebrick1") +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  edgeR_RNA_OHT.logFC > log2(1.2) &
                                  edgeR_RNA_OHT.FDR < .1),
                    aes(y = edgeR_RNA_OHT.logFC,
                        x = edgeR_RNA_OHT.logCPM),
                    color = "chartreuse3") +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP)),
                    aes(y = edgeR_RNA_OHT.logFC,
                        x = edgeR_RNA_OHT.logCPM),
                    color = "steelblue") +
         geom_label_repel(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP) &
                                  edgeR_RNA_OHT.logFC > 1 &
                                  edgeR_RNA_OHT.FDR < .1),
               aes(y = edgeR_RNA_OHT.logFC,
                   x = edgeR_RNA_OHT.logCPM,
                   label = SYMBOL)) +
         scale_x_continuous(name = "logCPM") +
         scale_y_continuous(name = "logFC 4OHT over DMSO") +
         ggtitle("Senescence induction : Impact on RNA") +
         theme_bw()

#--> Build MA plot for 3DA treatment
MA_DA <- ggplot() +
         geom_point(data = data.frame(gr_edger_rna),
                    aes(y = edgeR_RNA_3DA.logFC,
                        x = edgeR_RNA_3DA.logCPM)) +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  edgeR_RNA_3DA.logFC < -log2(1.2) & 
                                  edgeR_RNA_3DA.FDR < .1),
                    aes(y = edgeR_RNA_3DA.logFC,
                        x = edgeR_RNA_3DA.logCPM),
                    color = "firebrick1") +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  edgeR_RNA_3DA.logFC > log2(1.2) &
                                  edgeR_RNA_3DA.FDR < .1),
                    aes(y = edgeR_RNA_3DA.logFC,
                        x = edgeR_RNA_3DA.logCPM),
                    color = "chartreuse3") +
         geom_point(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP)),
                    aes(y = edgeR_RNA_3DA.logFC,
                        x = edgeR_RNA_3DA.logCPM),
                    color = "steelblue") +
         geom_label_repel(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP) &
                                  abs(edgeR_RNA_3DA.logFC) > .2 &
                                  edgeR_RNA_3DA.FDR < .1),
               aes(y = edgeR_RNA_3DA.logFC,
                   x = edgeR_RNA_3DA.logCPM,
                   label = SYMBOL)) +
         scale_x_continuous(name = "logCPM") +
         scale_y_continuous(name = "logFC 3DA over 4OHT") +
         ggtitle("3DA : Impact on RNA") +
         theme_bw()

#--> Build FC vs FC plot
FC_FC <- ggplot() +
    geom_point(data = data.frame(gr_edger_rna),
               aes(x = edgeR_RNA_3DA.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "grey") +
    geom_smooth(method = "gam", data = data.frame(gr_edger_rna),
               aes(x = edgeR_RNA_3DA.logFC, y = edgeR_RNA_OHT.logFC)) +
    geom_point(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP)),
                    aes(x = edgeR_RNA_3DA.logFC,
                        y = edgeR_RNA_OHT.logFC),
                    color = "steelblue") +
     geom_label_repel(data = subset(data.frame(gr_edger_rna),
                                  is.element(data.frame(gr_edger_rna)$SYMBOL, 
                                             gmt_senescence$SASP) &
                                  edgeR_RNA_OHT.FDR < .1 &
                                  edgeR_RNA_3DA.FDR < .1),
               aes(x = edgeR_RNA_3DA.logFC,
                   y = edgeR_RNA_OHT.logFC,
                   label = SYMBOL)) +
    xlab("log2 fold change in RNA (3DA over DMSO)") +
    ylab("log2 fold change in RNA (4OHT over 4OHT)")

library(gridExtra)
grid.arrange(MA_OHT, MA_DA, nrow = 2)

```

On the two MA-plot above, up-regulated genes during senescence or after 3DA treatment are colored in green, while the one down-regulated are colored in red. The overploted blue points depict SASP genes (as defined in the list provided by Jesus Gil).

As expected, we see that most of the SASP genes are up-regulated after OIS. Overall, these genes are mostly down-regulated genes upon 3DA treatment. The conclusion drown here for expression and H3K4me36 are also showing a "decoupling" between gene expression and H3K36me3 for SASP genes.

```{r, fig.cap = "Differential analysis for RNA-seq : FC vs FC plot", fig.height=8, results = "show"}
FC_FC
cor.test(gr_edger_rna$edgeR_RNA_3DA.logFC,
         gr_edger_rna$edgeR_RNA_OHT.logFC)
```

Comparing gene exoression fold-changes during senescence induction and upon 3DA treatment reveales a mild (but significant) negative correlation (R2 = -0.44) : genes upregulated after OIS tend to be down-regulated upon 3DA treatment. 

## Functionnal exploration

### GSEA for senenscence induction - REACTOME

```{r, fig.cap="Top up/down pathways for senenscence induction - REACTOME", fig.height=8, fig.width=8, fig.align='center'}
# Define the target 
target <- gr_edger_rna

# Import the GMT - REACTOME
gmt <-  qusage::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/c2.cp.reactome.v6.0.symbols.gmt")

# Order peaks according to the logFC of OHT over DMSO
order_da <- sort(target$edgeR_RNA_OHT.logFC, decreasing = TRUE)
names(order_da) <- target$SYMBOL[order(target$edgeR_RNA_OHT.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, BPPARAM = SerialParam(),  nperm = 1000)

# Select Top NES (+/-)
topPathways<- rbind(fgsea_da[order(NES)[c(1:5)],],
                    fgsea_da[order(NES,decreasing = TRUE)[c(1:5)],])
topPathways <- topPathways[order(topPathways$NES, decreasing = TRUE), pathway]


# Plot
plotGseaTable(gmt[topPathways], order_da, fgsea_da)

plotEnrichment(gmt[["REACTOME_CELL_CYCLE"]], order_da)
```

The senescence induction tend to impact positively (up-regulation in 4OHT samples) genes involved encoding for cytokines / chemokines, while impacting negatively (down-regulation in 4OHT samples) genes involved in cell cycle. 

### GSEA for senenscence induction - SENESCENCE

```{r, fig.cap="Top up/down pathways for senenscence induction - SENESCENCE", fig.height=8, fig.width=8, fig.align='center'}
# Define the target 
target <- gr_edger_rna

# Import the GMT - REACTOME
gmt <- gmt_senescence

# Order peaks according to the logFC of OHT over DMSO
order_da <- sort(target$edgeR_RNA_OHT.logFC, decreasing = TRUE)
names(order_da) <- target$SYMBOL[order(target$edgeR_RNA_OHT.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Plot
plotGseaTable(gmt, order_da, fgsea_da)

plotEnrichment(gmt[["SASP"]], order_da)
```

The senescence induction tend to impact positively (up-regulation in 4OHT samples) genes involved in the SASP. 

### GSEA for 3DA treatment - REACTOME

```{r, fig.cap="Top up/down pathways with 3DA treatment - REACTOME", fig.height=8, fig.width=8, fig.align='center'}
library(fgsea)

# Define the target 
target <- gr_edger_rna

# Import the GMT - REACTOME
gmt <- qusage::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/c2.cp.reactome.v6.0.symbols.gmt")

# Order peaks according to the logFC of 3DA over OHT
order_da <- sort(target$edgeR_RNA_3DA.logF, decreasing = TRUE)
names(order_da) <- target$SYMBOL[order(target$edgeR_RNA_3DA.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Select Top NES (+/-)
topPathways<- rbind(fgsea_da[order(NES)[c(1:10)],],
                    fgsea_da[order(NES,decreasing = TRUE)[c(1:10)],])
topPathways <- topPathways[order(topPathways$NES, decreasing = TRUE), pathway]

# Plot
plotGseaTable(gmt[topPathways], order_da, fgsea_da)


plotEnrichment(gmt[["REACTOME_CELL_CYCLE"]], order_da)

```

The 3DA treatment on senescent cells tend to impact positively  (up-regulation in 3DA samples)  genes involved in the cell cycle, while impacting negatively (down-regulation in 3DA samples) genes encoding for cytokines / chemokines. 

### GSEA for 3DA treatment - SENESCENCE

```{r, fig.cap="Top up/down pathways with 3DA treatment - SENESCENCE", fig.height=8, fig.width=8, fig.align='center'}
# Define the target 
target <- gr_edger_rna

# Import the GMT - REACTOME
gmt <- gmt_senescence

# Order peaks according to the logFC of 3DA over OHT
order_da <- sort(target$edgeR_RNA_3DA.logFC, decreasing = TRUE)
names(order_da) <- target$SYMBOL[order(target$edgeR_RNA_3DA.logFC, decreasing = TRUE)]

# Run GSEA
fgsea_da <- fgsea(gmt, order_da, nperm = 1000, BPPARAM = SerialParam(), minSize = 20)

# Plot
plotGseaTable(gmt, order_da, fgsea_da)

plotEnrichment(gmt[["SASP"]], order_da)
```

The 3DA treatment on senescent cells tend to impact nagatively  (down-regulation in 3DA samples)  genes involved in the SASP. 

# Integration

To explore the possible coupling between H3K36me3 signal in gene body and gene expression, we integrated the two data-sets and explore correlation between :

* Fold-change in H3K36me3 and gene expression after OIS.
* Fold-change in H3K36me3 and gene expression after 3DA treatment on senescent cells.

## Correlations between ChIP-seq and RNA-seq during OIS

```{r, fig.cap="FC vs FC plot for RAS OIS induction", fig.height=8, fig.width=8, fig.align='center', fig.height=15, results = "show"}

library(ggplot2)
library(ggrepel)

chip_rna <- merge(gr_edger_chip, gr_edger_rna, by.x = "entrez_id", by.y = "edgeR_RNA_OHT.genes")

fc_fc_ois <-ggplot() +
    geom_point(data = chip_rna,
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "grey") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_OHT.logFC > log2(1.5) &
                             edgeR_RNA_OHT.logFC >log2(1.5)),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "chartreuse3") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_OHT.logFC < - log2(1.5) &
                             edgeR_RNA_OHT.logFC <  - log2(1.5)),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "firebrick") +
    geom_label_repel(data = subset(chip_rna, edgeR_ChIP_OHT.logFC < - 2 &
                                   edgeR_RNA_OHT.logFC <  - 4),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC,
                   label = symbol), cex = 3) +
    geom_label_repel(data = subset(chip_rna, edgeR_ChIP_OHT.logFC >  2 &
                                   edgeR_RNA_OHT.logFC >   4),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC,
                   label = symbol), cex = 3) +
    geom_hline(yintercept = c(log2(1.5), -log2(1.5))) +
    geom_vline(xintercept = c(log2(1.5), -log2(1.5))) +
    geom_smooth(method = "lm", data = chip_rna,
                aes(x = edgeR_ChIP_OHT.logFC,
                    y = edgeR_RNA_OHT.logFC)) +
    ggtitle("Senescence induction : correlation between FC in expression and H3K36me3") +
    xlab("log2 fold change in H3K36me3 (OHT vs DMSO)") +
    ylab("log2 fold change in expression (OHT vs DMSO)")

fc_fc_ois_sasp <- ggplot() +
    geom_point(data = chip_rna,
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "grey") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_OHT.logFC > log2(1.2) &
                             edgeR_RNA_OHT.logFC >log2(1.2)),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "chartreuse3") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_OHT.logFC < - log2(1.2) &
                             edgeR_RNA_OHT.logFC <  - log2(1.2)),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC),
               color = "firebrick") +
    geom_point(data = subset(chip_rna,
                             is.element(chip_rna$SYMBOL, 
                                        gmt_senescence$SASP)),
                aes(x = edgeR_ChIP_OHT.logFC,
                    y = edgeR_RNA_OHT.logFC),
                color = "steelblue") +
    geom_label_repel(data = subset(chip_rna,
                             is.element(chip_rna$SYMBOL, 
                                        gmt_senescence$SASP) &
                                  edgeR_ChIP_OHT.logFC > 1 &
                                  edgeR_RNA_OHT.logFC > 1),
               aes(x = edgeR_ChIP_OHT.logFC,
                   y = edgeR_RNA_OHT.logFC,
                   label = SYMBOL)) + 
    geom_hline(yintercept = c(log2(1.2), -log2(1.2))) +
    geom_vline(xintercept = c(log2(1.2), -log2(1.2))) +
    geom_smooth(method = "lm", data = chip_rna,
                aes(x = edgeR_ChIP_OHT.logFC,
                    y = edgeR_RNA_OHT.logFC)) +
    ggtitle("Senescence induction : correlation between FC in expression and H3K36me3") +
    xlab("log2 fold change in H3K36me3 (OHT vs DMSO)") +
    ylab("log2 fold change in expression (OHT vs DMSO)")

grid.arrange(fc_fc_ois, fc_fc_ois_sasp, nrow = 2)

cor(chip_rna$edgeR_ChIP_OHT.logFC, chip_rna$edgeR_RNA_OHT.logFC)
```

The two plots above show the correlation between FC in H3K36me3 and gene expression during OIS. Overall, during RAS induction, there is a coordinate change in H3K36me3 and gene expression (R2=0.51) : genes gaining H3K36me3 after OIS tend to get up-regulated.

The graph at the shows top up / down coordinated genes. The graph at the bottom is the same, with SASP genes highlited. 

## Correlations between ChIP-seq and RNA-seq after 3DA

```{r, fig.cap="FC vs FC plot for 3DA treatment on senescent cells", fig.width=8, fig.align='center', fig.height=15, results = "show"}

fc_fc_da <- ggplot() +
    geom_point(data = chip_rna,
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "grey") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_3DA.logFC > log2(1.2) &
                             edgeR_RNA_3DA.logFC >log2(1.2)),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "chartreuse3") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_3DA.logFC < - log2(1.2) &
                             edgeR_RNA_3DA.logFC <  - log2(1.2)),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "firebrick") +
    geom_label_repel(data = subset(chip_rna, edgeR_ChIP_3DA.logFC < - 1 &
                                   edgeR_RNA_3DA.logFC <  - 1)[1:10,],
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC,
                   label = symbol), cex = 3) +
    geom_label_repel(data = subset(chip_rna, edgeR_ChIP_3DA.logFC >  1 &
                                   edgeR_RNA_3DA.logFC >   1)[1:10,],
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC,
                   label = symbol), cex = 3) +
    geom_hline(yintercept = c(log2(1.2), -log2(1.2))) +
    geom_vline(xintercept = c(log2(1.2), -log2(1.2))) +
    geom_smooth(method = "lm", data = chip_rna,
                aes(x = edgeR_ChIP_3DA.logFC,
                    y = edgeR_RNA_3DA.logFC)) +
    ggtitle("3DA treatment in senescent cells : correlation between FC in expression and H3K36me3") +
    xlab("log2 fold change in H3K36me3 (OHT vs DA)") +
    ylab("log2 fold change in expression (OHT vs DA)") 

fc_fc_da_sasp <-ggplot() +
    geom_point(data = chip_rna,
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "grey") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_3DA.logFC > log2(1.2) &
                             edgeR_RNA_3DA.logFC >log2(1.2)),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "chartreuse3") +
    geom_point(data = subset(chip_rna, edgeR_ChIP_3DA.logFC < - log2(1.2) &
                             edgeR_RNA_3DA.logFC <  - log2(1.2)),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC),
               color = "firebrick") +
    geom_point(data = subset(chip_rna,
                             is.element(chip_rna$SYMBOL, 
                                        gmt_senescence$SASP)),
                aes(x = edgeR_ChIP_3DA.logFC,
                    y = edgeR_RNA_3DA.logFC),
                color = "steelblue") +
    geom_label_repel(data = subset(chip_rna,
                             is.element(chip_rna$SYMBOL, 
                                        gmt_senescence$SASP) &
                                  edgeR_ChIP_3DA.logFC > 1 &
                                  edgeR_RNA_3DA.logFC < .1),
               aes(x = edgeR_ChIP_3DA.logFC,
                   y = edgeR_RNA_3DA.logFC,
                   label = SYMBOL)) +
    geom_hline(yintercept = c(log2(1.2), -log2(1.2))) +
    geom_vline(xintercept = c(log2(1.2), -log2(1.2))) +
    geom_smooth(method = "lm", data = chip_rna,
                aes(x = edgeR_ChIP_3DA.logFC,
                    y = edgeR_RNA_3DA.logFC)) +
    ggtitle("3DA treatment in senescent cells : correlation between FC in expression and H3K36me3") +
    xlab("log2 fold change in H3K36me3 (OHT vs DA)") +
    ylab("log2 fold change in expression (OHT vs DA)")

grid.arrange(fc_fc_da, fc_fc_da_sasp, nrow = 2)
cor.test(chip_rna$edgeR_ChIP_3DA.logFC, chip_rna$edgeR_RNA_3DA.logFC)

```

The two plots above show the correlation between FC in H3K36me3 and gene expression after 3DA treatment on senescent cells. Here, the coupling between H3K36me3 and gene expression is lost (R2=0.19) for almost all genes, except cell cycle-related genes, for which an up-regulation at the RNA level is supported by an increase in H3K36me3.

The graph at the shows top up / down coordinated genes. The graph at the bottom is the same, with SASP genes highlited.

## List of up-regulated genes during 3DA with a concordant H3K36me3 profile

Considering a threshold at log2(1.2) for both gene expression and H3K36me3 fold-change after 3DA treatment, here is the list of up-regulated genes with a concordant H3K36me3 profile.

```{r, results = "show"}
up_list <- noquote(unique(subset(chip_rna, edgeR_ChIP_3DA.logFC > log2(1.2) &
                              edgeR_RNA_3DA.logFC > log2(1.2))$SYMBOL))
up_list
```


## List of down-regulated genes with a concordant H3K36me3 profiles

```{r, results = "show"}
down_list <- noquote(unique(subset(chip_rna, edgeR_ChIP_3DA.logFC < - log2(1.2) &
                              edgeR_RNA_3DA.logFC < - log2(1.2))$SYMBOL))
down_list
```

Considering a threshold at -log2(1.2) for both gene expression and H3K36me3 fold-change after 3DA treatment, here is the list of down-regulated genes with a concordant H3K36me3 profile. Note the presence of IL1B, IL1R, TNFSF15 and NFKB in this list.


# Examples

## Prepare the Gviz session

```{r}
#--> Loading libraries
library(Gviz)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)

#--> Define the GenomeAxis layout
g_axe <- GenomeAxisTrack(col="black")

#--> Define the GeneRegion characteristics
g_reg <- BiomartGeneRegionTrack(genome = "hg19")

#--> List histone ChIP-seq normalized signal tracks
samples_chip <- bw[c(3,2,1)]
names <- c("DMSO", "4OHT", "3DA")

#--> Create DataTrack for each samples
colors  <- c("#FD8800", "#8FBD00", "#BA008C")
type <- "histogram"

s1 <- DataTrack(samples_chip[1], name = names[1], ylim = c(0,15), type = type, col.histogram= colors[1])
s2 <- DataTrack(samples_chip[2], name = names[2], ylim = c(0,15), type = type, col.histogram = colors[2])
s3 <- DataTrack(samples_chip[3], name = names[3], ylim = c(0,15), type = type, col.histogram = colors[3])

samples <- list(s1, s2, s3)
```

```{r,fig.cap="Normalized profiles of H3K4me1, H3K27ac, H3K4me3 and H3K27me3 ChIP-seq signal in the IL12RB2 locus"}
#--> Define the region of interest
#region <- GRanges("chr4", IRanges(103420486,103540459)) #NKB1
#region <- GRanges("chr2", IRanges(136595196,136636047)) #MCM6
region <- GRanges("chr1", IRanges(67771047,67864583)) #IL12RB2


#--> Setup the highlight
hl <- HighlightTrack(trackList = samples, 
                     start = as.numeric(start(region)),
                     end = as.numeric(end(region)),
                     chromosome = as.character(seqnames(region)))

#--> Setup the chromosome ideogram
ideo <- IdeogramTrack(genome = "hg19",
                      chromosome = as.character(seqnames(region )))

#--> Plot
plotTracks(c(hl, g_reg, ideo),
           sizes=c(rep(1,3),1, 1),
           from = as.numeric(start(region)) - 1000,
           to = as.numeric(end(region)) + 1000,
           chromosome = as.character(seqnames(region)), shape = "arrow",
           transcriptAnnotation = "symbol", collapseTranscripts = "meta")


```


# Clustering


## Perform Clustering

```{r, fig.cap="K-means clustering of H3K36me3 data", fig.width=6, fig.align='center', fig.height=10}
library(pheatmap)
library(viridis)

de_all_chip_1 <- gr_edger_chip[
                    which((abs(gr_edger_chip$edgeR_ChIP_3DA.logFC) > log2(2) & gr_edger_chip$edgeR_ChIP_3DA.FDR < .1) |
                      abs(gr_edger_chip$edgeR_ChIP_OHT.logFC) > log2(2) & gr_edger_chip$edgeR_ChIP_OHT.FDR < .1),]
    

data_full <- cpm(data_edger_chip)
row.names(data_full) <- c(data_edger_chip$genes$genes)
colnames(data_full) <- data_edger_chip$samples$group

de_all_chip_2  <- data_full[match(de_all_chip_1$edgeR_ChIP_OHT.genes, row.names(data_full), nomatch = 0),]
#de_all_chip <- do.call("cbind", by(data = t(de_all_chip), INDICES = colnames(de_all_chip), colMeans))
de_all_chip_2 <- t(scale(t(de_all_chip_2)))

set.seed(123)
k <- kmeans(de_all_chip_2, centers = 5)
row.names(de_all_chip_2[order(k$cluster),])

anno_r <- data.frame(Cluster = as.factor(k$cluster), row.names = make.names(names(k$cluster), unique = TRUE))

pheatmap(unique(de_all_chip_2[order(k$cluster),]),
         cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = FALSE, color = viridis(100),
         annotation_row = anno_r)

```

## Functional over-representation

```{r, fig.cap= "Functional over-representation analysis (MSigDB hallmark)", fig.height= 8}
#---> Get clusters
probes <- data.frame(k$cluster, peak = names(k$cluster))
probes <- merge(probes, gr_edger_chip,
                by.x = "peak", by.y = "edgeR_ChIP_OHT.genes")

#---> Get the corresponding Entrez_ID
fun_enrich_WGCNA <- split(probes$entrez_id, probes$k.cluster)

#--->  Over-representation analysis (MSigDB Hallmark)
library(ReactomePA)
library(clusterProfiler)
            
# Load the gmt file
GMT <- clusterProfiler::read.gmt("~/Documents/Post-Doc/A-Data/OUTILS/h.all.v6.0.entrez.gmt")

# Run the over-representation analysis 
fun_enrich_res <- data.frame(compareCluster(fun_enrich_WGCNA, fun='enricher',
                                            TERM2GENE=GMT, pvalueCutoff=1, qvalueCutoff = 1,
                                            minGSSize = 20))
# Filter the data
fun_enrich_res_2 <- data.frame()
for (i in unique(fun_enrich_res$Cluster)) {
  x <- fun_enrich_res[which(fun_enrich_res$Cluster == i),]
  x[which(x$qvalue > .2),]$qvalue <- NA
  x <- x[order(x$qvalue),]
  fun_enrich_res_2 <- rbind(fun_enrich_res_2, x)
}

fun_enrich_res_2 <- na.omit(fun_enrich_res_2)

# Clustering
library(reshape)
fun_enrich_res_tmp <- fun_enrich_res_2[,c("ID", "qvalue", "Cluster")]
fun_enrich_res_tmp <- cast(fun_enrich_res_tmp, Cluster~ID, value = "qvalue")
fun_enrich_res_tmp <- data.matrix(fun_enrich_res_tmp)
fun_enrich_res_tmp[is.na(fun_enrich_res_tmp)] <- 0
row.names(fun_enrich_res_tmp) <- fun_enrich_res_tmp[,1]
fun_enrich_res_tmp <- fun_enrich_res_tmp[,-1]
order_fun_enrich <- hclust(as.dist(1-cor(fun_enrich_res_tmp)), method = "ward.D2")$order

fun_enrich_res_2$Description <- factor(fun_enrich_res_2$Description,
                                       levels =colnames(fun_enrich_res_tmp)[order_fun_enrich] )
fun_enrich_res_2$Percent <- unlist(lapply(strsplit(fun_enrich_res_2$GeneRatio, "/"), function(x){as.numeric(x[1])/as.numeric(x[2])}*100))

# Plot
ggplot(fun_enrich_res_2, aes(y = ID, x = Cluster, fill = -log10(pvalue), size = Percent)) +
  geom_point(shape = 21, stroke = .5) +
  scale_fill_gradient(low = "firebrick1", high = "chartreuse3") +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(size = 8))

```

## Integration

```{r, fig.cap="Projection of RNA-seq data on K-means clustering obtained for H3K36me3 data", fig.width=6, fig.align='center', fig.height=10}
data_all <- unique(merge(probes,  normCounts(set_RUVs), by.x = "entrez_id", by.y = "row.names"))
data_all <- data_all[order(data_all$k.cluster),]
data_all[,c(29:40)] <- t(scale(t(log10(data_all[,c(29:40)]+1))))

anno_r <- data.frame(as.factor(data_all$k.cluster), row.names = row.names(data_all))

pheatmap(unique(data_all[,c(29:40)]), cluster_rows = FALSE, cluster_cols = FALSE, scale = "none", show_rownames = FALSE, color = viridis(100), annotation_row = anno_r)
```

# Session

```{r,results = "show"}
sessionInfo()
```