diffbind.Rmd

---
title: "ChIPSeq DiffBind"
author: "Harvard Chan Bioinformatics Core"
date: "`r Sys.Date()`"
output:
   html_document:
      code_folding: hide
      df_print: paged
      highlights: pygments
      number_sections: true
      self_contained: true
      theme: default
      toc: true
      toc_float:
         collapsed: true
         smooth_scroll: true
params:
  # Fill this file with the right paths to nfcore output
  # .qs file name for saving DiffBind Counts object
  # species = mouse or human
  params_file: params_diffbind-example.R
  project_file: ../information.R
  # chose the functions that matches the pipeline: chip or atac
  functions_file: ../libs/load_data_chip.R
  condition_of_interest: genotype
  numerator: cKO
  denominator: WT
  species: mouse
  counts_csv_fn: diffbind_counts.csv
  results_sig_anno_fn: diffbind_results_anno.csv
---
Template developed with materials from https://hbctraining.github.io/main/.

```{r, cache = FALSE, message = FALSE, warning=FALSE}
# This set up the working directory to this file so all files can be found
library(rstudioapi)
setwd(fs::path_dir(getSourceEditorContext()$path))

# NOTE: This code will check version, this is our recommendation, it may work
#.      other versions
stopifnot(R.version$major>= 4) # requires R4
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1") 
stopifnot(compareVersion(as.character(BiocManager::version()), "3.18")>=0)
```

This code is in this ![](https://img.shields.io/badge/status-alpha-yellow) revision.


```{r source_params, cache = FALSE, message = FALSE, warning=FALSE}
# 1. set up condition_of_interest parameter from parameter above or manually
#    this is used to color plots, it needs to be part of the metadata
# 2. Set input files in this file
source(params$params_file)
# 3. If you set up this file, project information will be printed below and
#.   it can be reused for other Rmd files.
source(params$project_file)
# 4. Load custom functions to load data from coldata/metrics/counts
source(params$functions_file)
```

# Overview

-   Project: `r project`
-   PI: `r PI`
-   Analyst: `r analyst`
-   Experiment: `r experiment`

# Methodology 

[DiffBind](https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf) 
is an R Bioconductor package which provides functions for processing 
DNA data enriched for genomic loci, including ChIPseq data enriched for sites 
where specific protein/DNA binding occurs or histone marks are enriched. 

DiffBind is mainly used for identifying sites that are differentially enriched 
between two or more sample groups. It works primarily with sets of peak calls 
('peaksets'), which are sets of genomic intervals representing candidate protein 
binding sites for each sample. It includes functions that support the processing 
of peaksets, including overlapping and merging peak sets across an entire dataset, 
counting sequencing reads in overlapping intervals in peak sets, and identifying 
statistically significantly differentially bound sites based on evidence of 
binding affinity (measured by differences in read densities). To this end it uses 
statistical routines developed  in an RNA-Seq context (primarily the Bioconductor packages [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) and [DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html)).

```{r load_libraries, cache = FALSE, message = FALSE, warning=FALSE}
library(tidyverse)
library(knitr)
# library(rtracklayer)
library(DESeq2)
library(DEGreport)
library(ggrepel)
# library(RColorBrewer)
library(DT)
library(pheatmap)
library(bcbioR)
library(janitor)
library(ChIPpeakAnno)
library(UpSetR)
library(DiffBind)
library(qs)
library(EnhancedVolcano)
library(ggprism)
library(ChIPseeker)
library(msigdbr)
library(fgsea)

if (params$species == 'mouse'){
  library(TxDb.Mmusculus.UCSC.mm10.knownGene)
  txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
  anno_db <- 'org.Mm.eg.db'
  library(org.Mm.eg.db)
} else if (params$species == "human"){
  library(TxDb.Hsapiens.UCSC.hg38.knownGene)
  txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  anno_db <- 'org.Hs.eg.db'
  library(org.Hs.eg.db)
}


colors=cb_friendly_cols(1:15)
ggplot2::theme_set(theme_prism(base_size = 14))
opts_chunk[["set"]](
    cache = FALSE,
    cache.lazy = FALSE,
    dev = c("png", "pdf"),
    error = TRUE,
    highlight = TRUE,
    message = FALSE,
    prompt = FALSE,
    tidy = FALSE,
    warning = FALSE,
    fig.height = 4)

# set seed for reproducibility
set.seed(1234567890L)

```


```{r sanitize-datatable}
sanitize_datatable = function(df, ...) {
 # remove dashes which cause wrapping
 DT::datatable(df, ..., rownames=gsub("-", "_", rownames(df)),
                   colnames=gsub("-", "_", colnames(df)))
}
```

```{r download example data, eval = params$params_file == 'params_diffbind-example.R'}
# if using example data, download it from github
new_samplesheet_fn <- paste('download', basename(diffbind_samplesheet_fn), sep = '_')
download.file(diffbind_samplesheet_fn, destfile = new_samplesheet_fn, mode = "wb")
diffbind_samplesheet_fn <- new_samplesheet_fn

new_counts_fn <- paste('download', basename(diffbind_counts_fn), sep = '_')
download.file(diffbind_counts_fn, destfile = new_counts_fn, mode = "wb")
diffbind_counts_fn <- new_counts_fn
```

# Samples and metadata
```{r load_coldata, message=F, warning=F}
coldata <- load_coldata(coldata_fn)
```

```{r make samplesheet, message = F, warning =F, eval = !file.exists(diffbind_samplesheet_fn)}
# DiffBind requires a very specific samplesheet in order to create the peak counts object, see https://www.rdocumentation.org/packages/DiffBind/versions/2.0.2/topics/dba for further details

# make_diffbind_samplesheet is a function provided by bcbioR to help assemble DiffBind's samplesheet 
# using the nf-core samplesheet and output. In the resulting DiffBind counts object, it
# encodes your condition of interest as "Condition" and the antibody as "Factor"

samplesheet <- make_diffbind_samplesheet(coldata, bam_dir, peaks_dir, params$condition_of_interest)
write_csv(samplesheet, diffbind_samplesheet_fn)

# if necessary, one additional covariate of interest can be encoded as "Tissue" 
```

```{r show_metadata}
samplesheet <- read_csv(diffbind_samplesheet_fn)
samplesheet %>% dplyr::select(SampleID, Replicate, Condition, Factor, ControlID) %>% sanitize_datatable()
```

# Calculate counts matrix

The first step is to read in a set of peaksets and associated metadata. 
This is done using the DiffBind sample sheet. Once the peaksets are read in, 
a merging function finds all overlapping peaks and derives a single set of unique 
genomic intervals covering all the supplied peaks (a consensus peakset for 
the experiment). A region is considered for the consensus set if it appears in 
more than two of the samples. This consensus set represents the overall set of 
candidate binding sites to be used in further analysis.

The next step is to take the alignment files and compute count information for 
each of the peaks/regions in the consensus set. In this step, for each of the 
consensus regions, DiffBind uses the number of aligned reads in the ChIP sample 
and the input sample to compute a normalized read count for each sample at every 
potential binding site. The peaks in the consensus peakset may be re-centered and 
trimmed based on calculating their summits (point of greatest read overlap) in 
order to provide more standardized peak intervals.

We then normalize the count matrix to adjust for varying 
library size, and we use the normalized counts for further analysis including PCA. 

```{r create diffbind counts object, eval = !file.exists(diffbind_counts_fn)}
diffbind_obj <- dba(sampleSheet = samplesheet, scoreCol = 5)

# This command may take several minutes. Recommend using multiple cores and lots of memory
diffbind_counts <- dba.count(diffbind_obj, bUseSummarizeOverlaps = TRUE, bParallel = T)

# save object when time-intensive command is finished, so that this cell only need run once
qsave(diffbind_counts, diffbind_counts_fn)
```

# PCA 

Principal Component Analysis (PCA) is a statistical technique used to simplify 
high-dimensional data by identifying patterns and reducing the number of variables. 
In the context of ChIPseq, PCA helps analyze large datasets containing information 
about thousands of binding locations across different samples (e.g., tissues, cells).

```{r PCA}
diffbind_counts <- qread(diffbind_counts_fn)

diffbind_norm <- dba.normalize(diffbind_counts)

norm_counts <- dba.peakset(diffbind_norm, bRetrieve=TRUE, DataType=DBA_DATA_FRAME) %>%
  mutate(peak = paste(CHR, START, END, sep = '_')) %>%
  dplyr::select(-CHR, -START, -END) 
rownames(norm_counts) <- norm_counts$peak
norm_counts <- norm_counts %>% dplyr::select(-peak) %>% as.matrix()
norm_counts_log <- log2(norm_counts + 1)
norm_counts_log_df <- norm_counts_log %>% as.data.frame() %>%
  rownames_to_column('peak')

write_csv(norm_counts_log_df, params$counts_csv_fn)

coldata_for_pca <- coldata[colnames(norm_counts), ]

stopifnot(all(colnames(norm_counts) == rownames(coldata_for_pca)))

degPCA(norm_counts_log, coldata_for_pca, condition = params$condition_of_interest) +
  scale_color_cb_friendly()
```


# Differentially Bound Peaks

A standardized differential analysis is performed using DiffBind and the DESeq2 package, 
including estimation of size factors and dispersions, fitting and testing the 
model, evaluating the supplied contrast, and shrinking the LFCs. A p-value and FDR 
is assigned to each candidate binding site indicating confidence that they are differentially bound.

We use [ChIPpeakAnno](https://bioconductor.org/packages/release/bioc/html/ChIPpeakAnno.html) 
to identify any gene features within 1000 bp of a differentially bound site. 


```{r DB analysis}
diffbind_norm <- dba.contrast(diffbind_norm, contrast = c('Condition', params$numerator, params$denominator))
results_obj <- dba.analyze(diffbind_norm, 
                           bBlacklist = F, # Use TRUE with your data
                           bGreylist = F)

results_report <- dba.report(results_obj, th = 1) 
results_report_sig <- dba.report(results_obj)

results <- results_report %>% as.data.frame()

```

```{r annotate DB peaks}

anno_data <- toGRanges(txdb, feature = 'gene')
results_anno_batch <- annotatePeakInBatch(results_report, 
                                              AnnotationData = anno_data,
                                              output = 'overlapping',
                                              maxgap = 1000)

results_anno_batch_df <- results_anno_batch %>% as.data.frame()

if(params$species == 'mouse'){
  entrez_to_symbol <- AnnotationDbi::select(org.Mm.eg.db, results_anno_batch_df$feature, 
                                            "ENTREZID", columns = 'SYMBOL') %>%
    filter(!is.na(ENTREZID)) %>% distinct()
} else if (params$species == 'human'){
  entrez_to_symbol <- AnnotationDbi::select(org.Hs.eg.db, results_anno_batch_df$feature, 
                                            "ENTREZID", columns = 'SYMBOL') %>%
    filter(!is.na(ENTREZID)) %>% distinct()
}

results_anno_batch_df <- results_anno_batch_df %>% 
  left_join(entrez_to_symbol %>% dplyr::select(feature = ENTREZID, gene_name = SYMBOL))

write_csv(results_anno_batch_df, params$results_sig_anno_fn)

```


## MA plot

This plot can help to:
- Identify Differential Binding: Sites that show a significant log-fold change (M value away from 0) indicate changes in binding between conditions.
- Assess Data Quality: The plot can help in identifying biases or systematic errors in the data. Ideally, most points should scatter around the M=0 line, indicating that there is no significant systematic difference between the conditions.
- Visualize data dispersion: The distribution of points along the A-axis gives a sense of the spread of binding levels and any patterns or anomalies in the dataset.

```{r MA plot}
results_for_ma <- results_anno_batch_df%>%
  mutate(peak = paste(seqnames, start, end, sep = '_')) %>%
  mutate(t = 0) %>% 
  dplyr::select(peak, AveExpr = Conc, logFC = Fold, P.Value = p.value, adj.P.Val = FDR, t)
degMA(as.DEGSet(results_for_ma, contrast = paste(params$numerator, params$denominator, sep = ' vs. ')))

```

## Table of differentially bound peaks

```{r DB table}

results_sig_anno_batch_df <- results_anno_batch_df %>% filter(FDR < 0.05)
results_sig_anno_batch_df %>% dplyr::select(names(results), feature, gene_name) %>%
  sanitize_datatable()

```


## Volcano plot

This volcano plot shows the binding sites that are significantly up- and down-regulated as a result of the analysis comparison. The points highlighted in purple are sites that have padj < 0.05 and a log2-fold change magnitude > 0.5. Points in blue have a padj > 0.05 and a log2-fold change magnitude > 0.5. Grey points are non-significant. The dashed lines correspond to the cutoff values of log2-fold change and padj that we have chosen.

```{r volcano, fig.height = 8}
results_mod <- results_sig_anno_batch_df %>% 
  mutate(Fold = replace(Fold, Fold < -5, -5)) %>% 
  mutate(Fold = replace(Fold, Fold > 5, 5)) %>%
  mutate(peak = paste(seqnames, start, end, sep = '_')) 
# show <- as.data.frame(results_mod[1:6, c("Fold", "FDR", "gene_name")])

show <- results_mod %>% filter(!is.na(gene_name)) %>% slice_min(n = 6, order_by = FDR)

results_mod <- results_mod %>% mutate(gene_name = ifelse(peak %in% show$peak  , gene_name, NA))
EnhancedVolcano(results_mod,
                lab= results_mod$gene_name,
                pCutoff = 0.05, 
                selectLab = c(show$gene_name),
                FCcutoff = 0.5,
                x = 'Fold',
                y = 'FDR', 
                title = paste(params$condition_of_interest, ':', params$numerator, 'vs', params$denominator),
                col=as.vector(colors[c("dark_grey", "light_blue",
                                       "purple", "purple")]),
                subtitle = "", drawConnectors = T,  max.overlaps = Inf) 

```

## Plot top peaks

We visualize the log2 normalized read counts at a few of the most differentially 
bound sites. 
```{r plot top peaks, fig.width = 8, fig.height = 6}
norm_counts_log_long <- norm_counts_log %>% as.data.frame() %>%
  rownames_to_column('peak') %>%
  pivot_longer(!peak, names_to = 'sample', values_to = 'norm_counts_log2') %>%
  left_join(coldata)

norm_counts_log_long_top <- norm_counts_log_long %>% filter(peak %in% show$peak)

ggplot(norm_counts_log_long_top, aes(x = .data[[params$condition_of_interest]], y = norm_counts_log2)) + 
  facet_wrap(~peak, scale = 'free_y') + geom_boxplot()
```

## Annotate DB peaks

We use the [ChIPseeker](https://www.bioconductor.org/packages/release/bioc/html/ChIPseeker.html) 
package to determine the genomic context of the differentially bound peaks and 
visualize these annotations. We consider the promoter region to be within 2000 bp in either direction of the TSS. 

```{r annotate, echo = F}

results_sig_anno <- annotatePeak(results_report_sig, 
                                 tssRegion = c(-2000, 2000), 
                                 TxDb = txdb, 
                                 annoDb = anno_db, 
                                 verbose = F)
results_sig_anno_df <- results_sig_anno %>% as.data.frame()

plotAnnoPie(results_sig_anno)

plotDistToTSS(results_sig_anno)

```

# Functional Enrichment

Over-Representation Analysis (ORA) is a statistical method used to determine whether a predefined set of genes (e.g., genes belonging to a specific biological pathway or function) is over-represented (or enriched) among a list of differentially bound genes (DEGs) from ChIP-seq. Adventages of ORA:

- Simplicity: Easy to perform and interpret.
- Biological Insight: Helps to identify pathways and processes that are significantly affected in the condition studied.
- Prior Knowledge Integration: Utilizes existing biological knowledge through predefined gene sets.

```{r get databases}
if(params$species == 'human'){
  all_in_life=get_databases()
} else if (params$species == 'mouse'){
  all_in_life = get_databases('Mus musculus')
}
```

```{r ora}

universe_mapping = results_anno_batch_df %>% 
    filter(!is.na(FDR), !is.na(feature)) %>% 
  dplyr::select(ENTREZID = feature, SYMBOL = gene_name) %>% distinct()

ora_input = results_anno_batch_df %>% 
  filter(!is.na(FDR), FDR < 0.01, abs(Fold) > 0.3, !is.na(feature)) %>% 
  dplyr::select(ENTREZID = feature, SYMBOL = gene_name) %>% distinct()
all = run_fora(ora_input, universe_mapping, all_in_life)

ora_input = results_anno_batch_df %>% 
  filter(!is.na(FDR), FDR < 0.01, Fold > 0.3, !is.na(feature)) %>% 
  dplyr::select(ENTREZID = feature, SYMBOL = gene_name) %>% distinct()
up = run_fora(ora_input, universe_mapping, all_in_life)

ora_input = results_anno_batch_df %>% 
  filter(!is.na(FDR), FDR < 0.01, Fold < -0.3, !is.na(feature)) %>% 
  dplyr::select(ENTREZID = feature, SYMBOL = gene_name) %>% distinct()
down = run_fora(ora_input, universe_mapping, all_in_life)
  
```


## Significant pathways using all DB genes

```{r all pathways}
all %>% sanitize_datatable()
```


## Significant pathways using increased DB genes

```{r up pathways}
up %>% sanitize_datatable()
```


## Significant pathways using decreased DB genes

```{r down pathways, results='asis'}
down %>% sanitize_datatable()
```

# R session

List and version of tools used for the report generation.

```{r}
sessionInfo()
```