1.1 Processing of RNA-seq data

1.2 Processing of ChIP-seq data

1.3 Start a project and define the paths to the data

Before starting any project, it is a good idea to setup the path to the project as a variable. For example, we could use:

projPath <- "/shared/projects/2514_ebiii_n2/chipseq"
## If you are not using RStudio, you can change the working directory with:
setwd(projPath) ##NOT RUN

If you use an Rmd file and RStudio, when you open the file, your working directory should already be set as the folder where this Rmd file is located. You can check it with:

getwd()

Here, we will use different files for our analyses. We can create a list containing the absolute paths or names to the different files and folders that we will use during the practice session.

params <- list(
  workdir = "/shared/projects/2514_ebiii_n2/chipseq",
  mapping.chipseq = "stats_mappingChIPseq.tsv",
  peakcalling.chipseq = "stats_peakCalling.tsv",
  brg1 = "BRG1siCTRL_CHIP-seq_peaks.narrowPeak",
  mitf = "MITF_CHIP-seq_peaks.narrowPeak",
  sox10 = "SOX10_CHIP-seq_peaks.narrowPeak",
  brg1.bw = "BRG1siCTRL_CHIP-seq.bw",
  mitf.bw = "MITF_CHIP-seq.bw",
  sox10.bw = "SOX10_CHIP-seq.bw",
  rnaseq = "RNAseq_diff_norm.rds",
  annot = "annot.rds"
  )

2.1 Importing peak data

Peak files are in narrowPeak format which is of the form:

1. chrom - Name of the chromosome (or contig, scaffold, etc.).
   
2. chromStart - The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0 (i.e. 0-based).
   
3. chromEnd - The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as chromStart=0, chromEnd=100, and span the bases numbered 0-99.
   
4. name - Name given to a region (preferably unique). Use “.” if noname is assigned.
   
5. score - Indicates how dark the peak will be displayed in the browser (0-1000). If all scores were “‘0”’ when the data were submitted to the DCC, the DCC assigned scores 1-1000 based on signal value. Ideally the average signalValue per base spread is between 100-1000.
   
6. strand - +/- to denote strand or orientation (whenever applicable). Use “.” if no orientation is assigned.
   
7. signalValue - Measurement of overall (usually, average) enrichment for the region.
   
8. pValue - Measurement of statistical significance (-log10). Use -1 if no pValue is assigned.
   
9. qValue - Measurement of statistical significance using false discovery rate (-log10). Use -1 if no qValue is assigned.
   
10. peak - Point-source called for this peak; 0-based offset from chromStart. Use -1 if no point-source called.

There are different ways to import such data. We could use a simple read.table function but we know that it would not convert the 0-based coordinates to the 1-based coordinates used in R/Bioconductor.

The ChIPseeker package provides a readPeakFile function that we could use like this:

##NOT RUN
library(ChIPseeker)
peaks <- list()
peaks[["BRG1"]] <- readPeakFile(file.path(params$workdir,"data", "ChIPseq", params$brg1), 
                                as="GRanges")

However, the function is essentially a wrapper to the read.delim function and does not assign useful column names to the imported data. So instead we will use the readNarrowPeak function from the genomation package to import the narrowPeak files obtained for the 3 factors (BRG1, MITF and SOX10):

library(genomation)
peaks <- list()
peaks[["BRG1"]] <- genomation::readNarrowPeak(file.path(params$workdir,
                                                        "data", "ChIPseq", 
                                                        params$brg1))
peaks[["MITF"]] <- readNarrowPeak(file.path(params$workdir,
                                            "data", "ChIPseq", 
                                            params$mitf))
peaks[["SOX10"]] <- readNarrowPeak(file.path(params$workdir,
                                             "data", "ChIPseq", 
                                             params$sox10))

The package also provides a readGeneric function that can adapt to different peak file formats. For example, the following command is equivalent to what we did above for BRG1:

## NOT RUN
peaks[["BRG1"]] <- readGeneric(file.path(params$workdir,"data", "ChIPseq", params$mitf),
                               meta.cols = list(name = 4, score = 5, signalValue = 7,
                                                pValue = 8, qValue = 9, peak = 10),
                               zero.based = TRUE)

Let’s take a look at how the data was imported:

peaks[["BRG1"]]

## GRanges object with 72024 ranges and 6 metadata columns:
##           seqnames            ranges strand |                   name     score
##              <Rle>         <IRanges>  <Rle> |            <character> <integer>
##       [1]     chr1     980402-981178      * | BRG1siCTRL_CHIP-seq_..        77
##       [2]     chr1     983476-984526      * | BRG1siCTRL_CHIP-seq_..       108
##       [3]     chr1   1000729-1001179      * | BRG1siCTRL_CHIP-seq_..        90
##       [4]     chr1   1001762-1002054      * | BRG1siCTRL_CHIP-seq_..        40
##       [5]     chr1   1020914-1021207      * | BRG1siCTRL_CHIP-seq_..       114
##       ...      ...               ...    ... .                    ...       ...
##   [72020]     chrY 18936757-18937116      * | BRG1siCTRL_CHIP-seq_..        37
##   [72021]     chrY 19577774-19578316      * | BRG1siCTRL_CHIP-seq_..        73
##   [72022]     chrY 19891214-19892094      * | BRG1siCTRL_CHIP-seq_..       114
##   [72023]     chrY 19892553-19893192      * | BRG1siCTRL_CHIP-seq_..        55
##   [72024]     chrY 21837647-21838039      * | BRG1siCTRL_CHIP-seq_..        55
##           signalValue    pvalue    qvalue      peak
##             <numeric> <numeric> <numeric> <integer>
##       [1]     5.81321  10.22971   7.71662       166
##       [2]     5.97792  13.59661  10.83201       384
##       [3]     6.65677  11.62894   9.01201       199
##       [4]     3.78100   6.17367   4.02109       209
##       [5]     8.11634  14.29797  11.46427       158
##       ...         ...       ...       ...       ...
##   [72020]     4.53175   5.93549   3.79936       167
##   [72021]     6.32404   9.87728   7.38650       165
##   [72022]     7.86030  14.29797  11.46427       711
##   [72023]     5.29987   7.83781   5.50942       390
##   [72024]     5.55592   7.83781   5.50942       152
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

Peaks are stored as GenomicRanges (or GRanges) objects; this is an R format which looks like the bed format, but is optimized in terms of memory requirements and speed of execution.

NOTE:: when using GRanges it is always a good idea to set the seqinfo metadata (info about chromosome names and length essentially). We will do that later below.

We can start by computing some basic statistics on the peak sets.

2.2 How many peaks were called?

Compute the number of peaks per dataset. We use here a function from the apply family which help to apply recursively a given function on elements of the object. We will often use these functions along the course.

The sapply() function typically takes a list or vector as input and applies the same function (here length) to each element:

sapply(peaks,length)

##  BRG1  MITF SOX10 
## 72024  9702  4538

Note that in this specific case, base R provides a lengths function (lengths with an s at the end) that does the same, i.e. gives the length of each element of a list object:

lengths(peaks)

Make a simple barplot showing the number of BRG1 peaks.

barplot(lengths(peaks), ylab = "Number of peaks")

# or:
# barplot(lengths(peaks), ylab = "Number of peaks", log="y")

Let’s create a barplot with ggplot2 out of this data.
Step by step:
Do not forget to load the ggplot2 library

# Load ggplot2 library
library(ggplot2)

• Create a data.frame with two columns: IP contains names of the chipped TFs and NbPeaks contains the number of peaks.

# create a table with the data to display
peak.lengths <- data.frame(IP = names(peaks),
                           NbPeaks = lengths(peaks))

• use ggplot2 with the appropriate geometric object geom_*

# make the barplot
ggplot(peak.lengths, aes(x = IP, y = NbPeaks)) +
         geom_col()

• Improve the barplot by changing axis labels and the theme

#make the barplot a bit nicer:
ggplot(peak.lengths, aes(x = IP, y = NbPeaks)) +
  geom_col() +
  labs(x = "ChIP", y = "Number of peaks") +
  theme_bw()

We can customize the plots by adding colors.

There are a number of color palettes that you can use in R and different ways to specify colors. For example, the Rcolorbrewer package is widely used and offers several color palettes:

# See:
library(RColorBrewer)
par(mar = c(3,4,2,2))
display.brewer.all()

The ggsci package offers color palettes based on scientific journals.

If you are searching for a specific color, take a look at this site.

Now lets add colors to the barplot.
Step by step:
- Add the new information fill=IP to let ggplot know that colors change based on chipped protein

ggplot(peak.lengths, aes(x = IP, y = NbPeaks, fill = IP)) +
         geom_col() +
         labs(x = "ChIP", y = "Number of peaks") +
         theme_bw()

We want to use colors from RColorBrewer library with the “Set1” color palette. To set your own colors to fill the plot, several functions are available (scale_fill_*), here we use scale_fill_brewer.

ggplot(peak.lengths, aes(x = IP, y = NbPeaks, fill = IP)) +
         geom_col()+
         scale_fill_brewer(palette="Set1") +
         labs(x = "ChIP", y = "Number of peaks") +
         theme_bw()

And now use one of the ggsci palette.

## Load ggsci library
library(ggsci)

## Make the plot with the NPG palette
ggplot(peak.lengths, aes(x = IP, y = NbPeaks, fill = IP)) +
         geom_col()+
         scale_fill_npg() +
         labs(x = "ChIP", y = "Number of peaks") +
         theme_bw()

2.3 How large are these peaks?

The peaks were loaded as a GRanges objects.
- To manipulate them we load the GenomicRanges library. Then we use the appropriate function to retrieve the width of the peaks from the GRanges and use a function to summarize the peak length (display the statistics like minimum, maximum, quartiles, etc.)

## we use the function width() from GenomicRanges
library(GenomicRanges)
summary(width(peaks$BRG1))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   254.0   308.0   427.0   578.7   671.0  9186.0

# you can also try with the skim function from the skimr package:
# skimr::skim(width(peaks$BRG1))

Create a simple boxplot of the peak sizes.

# use a log scale on the y axis to facilitate the comparisons
peak.width <- lapply(peaks, width)
boxplot(peak.width, ylab = "peak width", log = "y")

Now, create a nice looking boxplot with ggplot2. ggplot takes a data frame as input. We can either create a date frame, this is what we have done when we created the barplot. Here, we are going to use the package reshape2 which, among all its features, can create a data frame from other types of data, in particular list objects, via the melt function:

# Load the package
library(reshape2)
peak.width.table <- melt(peak.width)
head(peak.width.table)

##   value   L1
## 1   777 BRG1
## 2  1051 BRG1
## 3   451 BRG1
## 4   293 BRG1
## 5   294 BRG1
## 6   351 BRG1

Create a ggplot object with correct aesthetics to display a boxplot according the chipped TF (x-axis) and their width (y-axis). Then use appropriate geometric object geom_*.

# create boxplot
ggplot(peak.width.table, aes(x = L1, y = value)) +
         geom_boxplot()

By default, the background color is grey and often does not allow to highlight properly the graphs. Many theme are available defining different graphics parameters, they are often added with a function theme_*(), the function theme() allows to tune your plot parameter by parameter. We propose here to use theme_classic() that changes a grey background to white background. Some peaks are very long and squeeze the distribution, one way to make the distribution easier to visualize is to transform the y axis to log scale. Finally, we want to remove the x label, change the y label to ‘Peak widths’ and the legend relative to the fill color from “L1” to “TF”.
Using the ggplot2 documentation (e.g. tidyverse or google it, try to enhance the boxplot.

# - theme_classic() change grey background to white background
# - fill=L1 and scale_fill_brewer(palette="Set1") colors boxplots
# based on chipped protein and with colors from RColorBrewer Set1 palette
# - labs changes x and y axis labels and legend title
# - scale_y_log10() set y axis to a log scale so that we can have a nice
# view of the data in small values
ggplot(peak.width.table, aes(x = L1, y = value, fill = L1)) +
         geom_boxplot()+
         theme_classic()+
         scale_fill_brewer(palette = "Set1")+
         labs(x = "", y = "Peak sizes", fill = "TF")+
         scale_y_log10()

2.4 Peak filtering

To make sure we keep only high quality data. We are going to select those peaks having a \(qvalue >= 4\). In the GRanges metadata columns (mcols) one of them is the qvalue. So, we are going to set a threshold on this column. How many peaks are selected ?

## Select BRG1 peaks with qvalue>=4 and count them
length( peaks$BRG1[peaks$BRG1$qvalue >= 4] )

## [1] 60184

Now we want to apply this function to all IPs.

## Select BRG1 peaks with qvalue>=4 and count them
peakfilt <- lapply(peaks,
                   function(x) {
                     x[x$qvalue >= 4]
                   })
lengths(peakfilt)

##  BRG1  MITF SOX10 
## 60184  9279  4214

3.1 Where are the peaks located across the whole genome?

Sometime, peaks may occur more on some chromosomes than others or there can be regions where peaks are absent (example: repetitive regions that have been filtered out). We can display the genomic distribution of peaks along the chromosomes, using the covplot function from ChIPSeeker. Height of peaks is drawn based on the peak scores.
Let’s use this function on SOX10 for example:

3.2 Where are the peaks located relative to genomic features?

Here, the aim is to identify where the peaks tend to be located relative to typical genome annotations such as promoters, exons, introns, intergenic regions, etc.
We will assign peaks to genomic features (genes, introns, exons, promoters, distal regions, etc.) by checking if they overlap with these features.
We will illustrate this task using functions from the ChIPseeker package.
But first we need to retrieve the annotations

## org.Hs.eg.db is an R object that contains mappings between Entrez Gene identifiers and GenBank accession numbers.
library(org.Hs.eg.db)
## Load transcript-oriented annotations
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

peaks <- lapply(peaks, function(x) {seqinfo(x) <- seqinfo(txdb); return(x)})
peakfilt <- lapply(peakfilt, function(x) {seqinfo(x) <- seqinfo(txdb); return(x)})

## Annotate peaks for all datasets and store it in a list
## Here TSS regions are regions -1000Kb/+100b arount TSS positions
## Peak annotations are stored in a list
peakAnno <- list()
peakAnno[["BRG1"]] = annotatePeak(peakfilt$BRG1, 
                                  tssRegion = c(-1000, 100), 
                                  TxDb = txdb, 
                                  annoDb = "org.Hs.eg.db")

## >> preparing features information...      2025-05-29 11:10:43 AM 
## >> identifying nearest features...        2025-05-29 11:10:44 AM 
## >> calculating distance from peak to TSS...   2025-05-29 11:10:46 AM 
## >> assigning genomic annotation...        2025-05-29 11:10:46 AM 
## >> adding gene annotation...          2025-05-29 11:11:12 AM

## >> assigning chromosome lengths           2025-05-29 11:11:12 AM 
## >> done...                    2025-05-29 11:11:12 AM

class(peakAnno$BRG1)

## [1] "csAnno"
## attr(,"package")
## [1] "ChIPseeker"

## Visualize and export annotation as a data table
# as.data.frame(peakAnno$BRG1)
head(as.data.frame(peakAnno$BRG1))

##   seqnames   start     end width strand                       name score
## 1     chr1  980402  981178   777      * BRG1siCTRL_CHIP-seq_peak_1    77
## 2     chr1  983476  984526  1051      * BRG1siCTRL_CHIP-seq_peak_2   108
## 3     chr1 1000729 1001179   451      * BRG1siCTRL_CHIP-seq_peak_3    90
## 4     chr1 1001762 1002054   293      * BRG1siCTRL_CHIP-seq_peak_4    40
## 5     chr1 1020914 1021207   294      * BRG1siCTRL_CHIP-seq_peak_5   114
## 6     chr1 1059581 1059989   409      * BRG1siCTRL_CHIP-seq_peak_8    68
##   signalValue   pvalue   qvalue peak
## 1     5.81321 10.22971  7.71662  166
## 2     5.97792 13.59661 10.83201  384
## 3     6.65677 11.62894  9.01201  199
## 4     3.78100  6.17367  4.02109  209
## 5     8.11634 14.29797 11.46427  158
## 6     5.43367  9.31879  6.86163  239
##                                          annotation geneChr geneStart geneEnd
## 1                                            5' UTR       1    975198  982093
## 2                                 Distal Intergenic       1    976121  982117
## 3                                          Promoter       1   1001138 1014540
## 4    Intron (ENST00000624697.4/9636, intron 1 of 2)       1   1001145 1014435
## 5 Intron (ENST00000379370.7/375790, intron 1 of 35)       1   1020123 1056118
## 6                                          Promoter       1   1059687 1066449
##   geneLength geneStrand    geneId      transcriptId distanceToTSS
## 1       6896          2     84808 ENST00000433179.4           915
## 2       5997          2     84808 ENST00000694917.1         -1359
## 3      13403          1      9636 ENST00000624697.4             0
## 4      13291          1      9636 ENST00000624652.1           617
## 5      35996          1    375790 ENST00000620552.4           791
## 6       6763          1 100288175 ENST00000688131.1             0
##           ENSEMBL       SYMBOL                                     GENENAME
## 1 ENSG00000187642        PERM1 PPARGC1 and ESRR induced regulator, muscle 1
## 2 ENSG00000187642        PERM1 PPARGC1 and ESRR induced regulator, muscle 1
## 3 ENSG00000187608        ISG15                ISG15 ubiquitin like modifier
## 4 ENSG00000187608        ISG15                ISG15 ubiquitin like modifier
## 5 ENSG00000188157         AGRN                                        agrin
## 6 ENSG00000291156 LOC100288175                 uncharacterized LOC100288175

## NOT RUN
## download GTF file
download.file("https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/93934/101/GCF_001577835.2_Coturnix_japonica_2.1/GCF_001577835.2_Coturnix_japonica_2.1_genomic.gtf.gz", "Coturnix_japonica_2.1.annotation.gtf.gz")
## Build TxDb object
library(GenomicFeatures)
txdb = makeTxDbFromGFF("Coturnix_japonica_2.1.annotation.gtf.gz")
## To save the txdb database
library(AnnotationDbi)
saveDb(txdb, 'txdb.Coturnix_japonica_2.1.sqlite')
## load it when needed
library(AnnotationDbi)
txdb = loadDb(file = 'txdb.Coturnix_japonica_2.1.sqlite')

## distribution of genomic features for BRG1 peaks
# as a pie chart - which is the most widely used representation in publication
plotAnnoPie(peakAnno$BRG1)

## Use the annotatePeak function
peakAnno[["SOX10"]] <- annotatePeak(peakfilt$SOX10, 
                                    tssRegion = c(-1000, 100), 
                                    TxDb = txdb, 
                                    annoDb = "org.Hs.eg.db")

## >> preparing features information...      2025-05-29 11:11:16 AM 
## >> identifying nearest features...        2025-05-29 11:11:16 AM 
## >> calculating distance from peak to TSS...   2025-05-29 11:11:17 AM 
## >> assigning genomic annotation...        2025-05-29 11:11:17 AM 
## >> adding gene annotation...          2025-05-29 11:11:20 AM

## >> assigning chromosome lengths           2025-05-29 11:11:20 AM 
## >> done...                    2025-05-29 11:11:20 AM

peakAnno[["MITF"]] <- annotatePeak(peakfilt$MITF, 
                                   tssRegion = c(-1000, 100), 
                                   TxDb = txdb, 
                                   annoDb = "org.Hs.eg.db")

## >> preparing features information...      2025-05-29 11:11:20 AM 
## >> identifying nearest features...        2025-05-29 11:11:20 AM 
## >> calculating distance from peak to TSS...   2025-05-29 11:11:20 AM 
## >> assigning genomic annotation...        2025-05-29 11:11:20 AM 
## >> adding gene annotation...          2025-05-29 11:11:23 AM

## >> assigning chromosome lengths           2025-05-29 11:11:24 AM 
## >> done...                    2025-05-29 11:11:24 AM

## We could also have generated annotations for all factors at the same time using:
# peakAnno <- lapply(peakfilt, annotatePeak,  tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")

## Use the plotAnnoBar function
plotAnnoBar(peakAnno)

upsetplot(peakAnno$BRG1)

3.3 Do the peaks of the different factors overlap?

library(ChIPpeakAnno)
ovl <- findOverlapsOfPeaks(peakfilt, connectedPeaks = "keepAll")

# Estimate the average size of the peaks ...
averagePeakWidth <- mean(width(unlist(GRangesList(ovl$peaklist))))
# ... to count how many potential sites could have been bound in coding regions.
tot <- ceiling(3.3e+9 * 0.1 / averagePeakWidth)

makeVennDiagram(ovl, totalTest = tot, connectedPeaks = "keepAll",
                fill = brewer.pal(3, "Set1"), # circle fill color
                col = brewer.pal(3, "Set1"), #circle border color
                cat.col = brewer.pal(3, "Set1"))

# Prepare a pool of random peaks following the characteristics of our peak sets
pool <- preparePool(txdb, peaks$SOX10,
                    bindingType = "TSS",
                    featureType = "transcript",
                    seqn = paste0("chr",c(1:22, "X", "Y")))

# Create the permPool object needed for peakPermTest
pool <- new("permPool", grs=pool$grs[1], N=length(peaks$SOX10))
SOX10.MITF <- peakPermTest(peaks$SOX10, peaks$MITF, 
                           pool = pool, 
                           seed = 123, 
                           force.parallel = FALSE)
SOX10.MITF

## $cntOverlaps
## P-value: 0.0099009900990099
## Z-score: 149.1501
## Number of iterations: 100
## Alternative: greater
## Evaluation of the original region set: 535
## Evaluation function: cntOverlaps
## Randomization function: randPeaks
## 
## attr(,"class")
## [1] "permTestResultsList"

plot(SOX10.MITF)

coocs <- as(ovl$peaklist, "GRangesList")
# apply annotatePeak to each set of peaks in the list coocs
peakAnnoList <- lapply(coocs, annotatePeak, 
                       TxDb = txdb, annoDb = "org.Hs.eg.db",
                       tssRegion = c(-1000, 100),
                       verbose = FALSE)

# Plot peak distribution relative to gene features
plotAnnoBar(peakAnnoList)

# Plot peak distribution relatively to their distance to the TSS
plotDistToTSS(peakAnnoList)

4.1 Distribution of BRG1 peaks

We want to know if BRG1 is binding large and/or narrow regions, unique and/or tandem etc. So we will study the distribution of the BRG1 peaks that have been called in a region spanning +/-2Kb around the peak centers.

First of all, for calculating the profile of ChIP peaks around the BRG1 peak center, we need to define BRG1 peak centers and extend them 2Kb on each side. We will then select the 10000 best peaks using the q-value to get a lighter matrix. Then, we need to “align” the peaks that are mapped to these regions, and thus generate a matrix with the different regions in row and the position around the peak centers in columns. The resulting matrix will be 2000*2+1 columns and 10000 rows.

Step by step:
- Compute the centers of BRG1 peaks using the GenomicRanges `resize`` function.

# Compute the centers of BRG1 peaks
peak_centers_BRG1 <- resize(peakfilt$BRG1, width = 1, fix = "center")

• Now select the BRG1 peak centers that correspond to the 10,000 most significant BRG1 peaks (based on qvalue)

# For computation and memory efficiency reasons,
# we subset the top 10K peaks according to the qvalue
top10k_BRG1 <- peak_centers_BRG1[order(peak_centers_BRG1$qvalue, 
                                       decreasing = TRUE)][1:10000]

• Create a GRanges object with the regions corresponding to +/-2kb around these peak centers.

# Generate GRanges of peak centers +/-2kb for the 10K top peaks
top10k_BRG1_extend2kb <- top10k_BRG1 + 2000
# Note that we could have obtained this GRanges for all peaks directly using:
# BRG1_extend2kb <- resize(peakfilt$BRG1, width = 4001, fix = "center"))

BRG1_extend_2kb <- makeBioRegionFromGranges(top10k_BRG1,
                                            by = "peaks",
                                            type = "start_site",
                                            upstream = 2000,
                                            downstream = 2000)

## compute the coverage of peaks in these regions
BRG1_tagmat <- getTagMatrix(peakfilt$BRG1, 
                            windows = BRG1_extend_2kb)

## >> preparing start_site regions by peaks... 2025-05-29 10:56:26 AM
## >> preparing tag matrix...  2025-05-29 10:56:26 AM

## plot the heatmap
tagHeatmap(BRG1_tagmat)

## plot the heatmap of BRG1 peaks around the TSS
peakHeatmap(peakfilt$BRG1, TxDb = txdb, 
            upstream = 2000, downstream = 2000)

## >> preparing promoter regions...  2025-05-29 11:12:09 AM 
## >> preparing tag matrix...        2025-05-29 11:12:09 AM 
## >> preparing start_site regions by gene... 2025-05-29 11:12:10 AM
## >> preparing tag matrix...  2025-05-29 11:12:10 AM 
## >> generating figure...       2025-05-29 11:12:17 AM 
## >> done...            2025-05-29 11:12:20 AM

## compute the coverage of peaks in these regions
plotAvgProf(BRG1_tagmat,
            xlim = c(-2000, 2000),
            xlab = "Distance to peak center",
            ylab = "Peak Count Frequency")

## >> plotting figure...             2025-05-29 11:14:31 AM

## Prepare the matrix of BRG1 peaks across BRG1 peak centers +/-2kb
sm_BRG1 <- ScoreMatrix(target = peakfilt$BRG1, 
                       windows = top10k_BRG1_extend2kb)

## Plot the heatmap
heatMatrix(sm_BRG1, xcoords = c(-2000, 2000))

## Plot the heatmap using qvalue as weight
heatMatrix(ScoreMatrix(target = peakfilt$BRG1, 
                       windows = top10k_BRG1_extend2kb, 
                       weight.col = "qvalue"), 
           xcoords = c(-2000, 2000))

## Plot the heatmap with reordered rows
heatMatrix(sm_BRG1, xcoords = c(-2000, 2000), order = TRUE)

## we also add a confidence interval around the mean
heatMatrix(ScoreMatrix(target = peakfilt$BRG1, 
                       windows = promoters(genes(txdb),
                                           upstream = 2000,
                                           downstream = 2001)), 
           xcoords = c(-2000, 2000),
           order = TRUE)

## Average profile with genomation
## we also add a confidence interval around the mean
plotMeta(sm_BRG1, xcoords = c(-2000, 2000), dispersion = "se")

## Obtain regions of interest (SOX10 peak centers +/- 2kb)
SOX10_extend2kb <- resize(peakfilt$SOX10, width = 4001, fix = "center")

## Get the coverage of peaks for all factors around SOX10 peak centers
sm_sox10peaks <-  ScoreMatrixList(target = peakfilt[c("SOX10", "BRG1", "MITF")], 
                                  windows = SOX10_extend2kb)

## Get the coverage of peaks for all factors around SOX10 peak centers
## We order the rows by the average signal across all matrices and use a common color scale for all heatmaps
multiHeatMatrix(sm_sox10peaks, 
                xcoords = c(-2000, 2000),
                 order = TRUE,
                common.scale = TRUE)

4.2 Read enrichment in regions of interest

So far, we’re only looking at heatmaps of peak presence/absence. Even if we sort the regions by their intensity or by qvalue, this gives little quantitative information on the actual links between the factors.
To dig into this, we will look at the coverage of reads, instead of peaks, in specific regions.

## Obtain regions of interest (BRG1 peak centers +/- 5kb)
extended.5K.BRG1 <- top10k_BRG1 + 5000

## Import BRG1 bigwig file as an RleList
cvg.BRG1 <- import(file.path(params$workdir, "data", "ChIPseq", params$brg1.bw),
                      format = "BigWig",
                      which = extended.5K.BRG1,
                      as = "RleList")
cvg.BRG1

## RleList of length 25
## $chr1
## numeric-Rle of length 248956422 with 1776460 runs
##   Lengths: 975789     11      5     35      2 ...   1288     73     36  22594
##   Values :   0.00   0.32   0.81   1.13   0.97 ...   0.00   0.17   0.33   0.00
## 
## $chr10
## numeric-Rle of length 133797422 with 693950 runs
##   Lengths: 270215     30     64     37      6 ...     90     58    154 272220
##   Values :   0.00   0.65   0.97   0.81   0.65 ...   0.65   0.00   0.33   0.00
## 
## $chr11
## numeric-Rle of length 135086622 with 740549 runs
##   Lengths: 202802     21    117    221     41 ...      9     70    154  74826
##   Values :   0.00   0.49   0.16   0.00   0.16 ...   0.16   0.00   0.16   0.00
## 
## $chr12
## numeric-Rle of length 133275309 with 769359 runs
##   Lengths: 159158    105    616     25    129 ...     24    372    154  89209
##   Values :   0.00   0.16   0.00   0.16   0.49 ...   0.16   0.00   0.16   0.00
## 
## $chr13
## numeric-Rle of length 114364328 with 467206 runs
##   Lengths: 18172785       58        1       21 ...       44      110    40729
##   Values :     0.00     0.16     0.65     0.97 ...     0.49     0.16     0.00
## 
## ...
## <20 more elements>

## Obtain the matrix with 200 bins
smr_BRG1 <- ScoreMatrixBin(target = cvg.BRG1, 
                           windows = extended.5K.BRG1, 
                           bin.num = 200)

## Plot the heatmap, ordering the regions by their average intensity
heatMatrix(smr_BRG1, 
           xcoords = c(-5000, 5000),
           order = TRUE) 

## Plot the heatmap, ordering the regions by their average intensity
heatMatrix(smr_BRG1, 
           xcoords = c(-5000, 5000),
           winsorize = c(0, 99.9),
           order = TRUE) 

## Plot the heatmap for the scaled matrix
heatMatrix(scaleScoreMatrix(smr_BRG1), 
           xcoords = c(-5000, 5000), 
           order = TRUE) 

## Make 5 groups of peaks based on the pattern of read coverage around the peaks
set.seed(123) # for reproducibility since we're using k-means clustering
heatMatrix(smr_BRG1, 
           xcoords = c(-5000, 5000), 
           winsorize = c(0, 99.9),
           clustfun = function(x) kmeans(x, centers=5)$cluster) 

## Obtain regions of interest (SOX10 peak centers +/- 2kb)
SOX10_extend3kb_sorted <- resize(peakfilt$SOX10[order(peakfilt$SOX10$qvalue, 
                                                      decreasing = TRUE)],
                                 width = 6001,
                                 fix = "center")

## Import read coverage for the 3 factors
cvg <- list()
cvg[["SOX10"]] <- readBigWig(file.path(params$workdir, "data", "ChIPseq", params$sox10.bw),
                             windows = SOX10_extend3kb_sorted)
cvg[["BRG1"]] <- readBigWig(file.path(params$workdir, "data", "ChIPseq", params$brg1.bw),
                            windows = SOX10_extend3kb_sorted)
cvg[["MITF"]] <- readBigWig(file.path(params$workdir, "data", "ChIPseq", params$mitf.bw),
                            windows = SOX10_extend3kb_sorted)

## Obtain the score matrices
smr_at_SOX10peaks <- ScoreMatrixList(target = cvg, 
                                     windows = SOX10_extend3kb_sorted, 
                                     bin.num = 300)

## Plot the heatmaps side-by-side
multiHeatMatrix(smr_at_SOX10peaks, 
                xcoords = c(-3000, 3000), 
                winsorize = c(0, 99.9))

4.2.3 Customized heatmaps with the EnrichedHeatmap package

So far, we have used functions from ChIPseeker and genomation to plot the heatmaps. While these functions offer some level of customization, you may need further customization in certain cases. The EnrichedHeatmap package offers an almost complete level of customization to plot heatmaps and heatmap annotations from functional genomic experiments.
The package is built on the ComplexHeatmap package, which has an extensive documentation and examples.

As an example, let’s plot the enrichment of SOX10 and BRG1 reads around SOX10 peaks.

Step by step:
- First convert the ScoreMatrix objects to normalizedMatrix objects that are used by EnrichedHeatmap. You can create a normalizedMatrix using the normalizeToMatrix function but since we already computed the matrices, we will instead convert them using as.normalizedMatrix

Show code: Convert a ScoreMatrix to a normalizedMatrix

## Convert signal matrices to normalizedMatrix objects
normat <- list()

## For BRG1 signal:
normat[["BRG1"]] <- as.normalizedMatrix(
  as(smr_at_SOX10peaks[["BRG1"]], "matrix"),
  k_upstream = 150,
  k_downstream = 150,
  k_target = 0,
  extend = c(3000, 3000),
  signal_name = "BRG1",
  target_name = "SOX10 peak centers"
)

## For SOX10 signal:
normat[["SOX10"]] <- as.normalizedMatrix(
  as(smr_at_SOX10peaks[["SOX10"]], "matrix"),
  k_upstream = 150,
  k_downstream = 150,
  k_target = 0,
  extend = c(3000, 3000),
  signal_name = "SOX10",
  target_name = "SOX10 peak centers"
)

• Plot BRG1 heatmap using default parameters

Show code: Plot BRG1 EnrichedHeatmap

## BRG1 EnrichedHeatmap
EnrichedHeatmap(normat$BRG1, name = "BRG1")

We notice that the heatmap is sorted by increasing average signal across the regions.

• Now use the syntax of EnrichedHeatmap to plot both SOX10 and BRG1 heatmaps

Show code: Plot SOX10 + BRG1 EnrichedHeatmaps

## SOX10 + BRG1 EnrichedHeatmap
EnrichedHeatmap(normat$SOX10, name = "SOX10") +
  EnrichedHeatmap(normat$BRG1, name = "BRG1")

The first heatmap (SOX10) is ordered by increasing average signal and the second heatmap (BRG1) follows the same order as the first one so that we can compare the signal of the 2 factors in the same regions.

• Now let’s cluster the BRG1 signal into 5 groups using kmeans

Show code: cluster BRG1 signal into 8 groups

## make 5 clusters of row using kmeans
set.seed(1234) #for reproducibility
km5 <- kmeans(as.matrix(normat$BRG1),
              centers=5)
# keep the asignment of the rows to clusters
part5clust = km5$cluster

• Take a look at the quantiles of the 2 signal matrices to choose the color scales

Show code: Based on the quantiles of the signal matrices create color scales for the heatmaps. The color scales are created with the colorRamp2 function from the circlize package which linearly interpolates colors according to break values.

## Look at the quantiles of the BRG1 matrix
quantile(as(smr_at_SOX10peaks[["BRG1"]], "matrix"), seq(0.95,1,0.01))

##      95%      96%      97%      98%      99%     100% 
##   8.8985   9.7420  10.8640  12.4920  15.5355 105.4010

## and for SOX10
quantile(as(smr_at_SOX10peaks[["SOX10"]], "matrix"), seq(0.95,1,0.01))

##        95%        96%        97%        98%        99%       100% 
##   3.436025   4.040000   4.770000   5.814500   7.652500 162.255500

## We pick 15 as the max color value (~99th quantile of the matrix) for BRG1
col_brg1 <- colorRamp2(breaks = c(0,10,15), 
                       colors = c("white", "blue","black"))
## and 8 for SOX10
col_sox10 <- colorRamp2(breaks = c(0,4,8), 
                        colors = c("white", "orange","red"))

• Plot the heatmaps of the partitioned SOX10 and BRG1 signals at SOX10 binding sites using EnrichedHeatmap

Show code: EnrichedHeatmaps of SOX10 and BRG1

## Create a vector of the colors that will be used for the clusters
clustcolors <- brewer.pal(5,"Set1")
names(clustcolors) <- 1:5

## Create a legend object for the partition
lgd <- Legend(at = paste0("cluster", 1:5),
              title = "Clusters",
              type = "lines", legend_gp = gpar(col = clustcolors))


## In the plot we want:
## 1) the partition in cluster 1 to 8
## 2) the heatmap of BRG1 with the average profiles for the different clusters
## 3) the heatmap of SOX10 with the average profiles for the different clusters
## Define the heatmap list with these 3 objects:
ht_list <- Heatmap(part5clust, 
                   col = clustcolors, 
                   cluster_rows = FALSE,
                   name = " ", 
                   width = unit(3, "mm"),
                   show_heatmap_legend = FALSE) +
  EnrichedHeatmap(normat$BRG1, 
                  name = "BRG1", 
                  col=col_brg1,
                  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = clustcolors))),
                  column_title = "BRG1",
                  axis_name = c("-3kb", "SOX10\npeak center", "+3kb")) +
  EnrichedHeatmap(normat$SOX10, 
                  name = "SOX10", 
                  col=col_sox10,
                  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = clustcolors))),
                  column_title = "SOX10",
                  axis_name = c("-3kb", "SOX10\npeak center", "+3kb"))

## Draw the heatmap while spliting the rows based on the clusters:
draw(ht_list, 
     split = part5clust,
     annotation_legend_list = list(lgd),
     ht_gap = unit(c(2, 8, 8), "mm"))

6.1 GO enrichment analysis of genes in cluster 2

We want to know if the genes in cluster 2, which are characterized by a diffuse binding pattern of BRG1 at their TSS and a tendency to be upregulated upon BRG1 knockdown, are enriched for specific Gene Ontology (GO) Biological Process annotations.

The clusterprofiler package provides very useful functions to perform either classical enrichment analysis of functional categories in gene lists or gene set enrichment analyses.
Here we will use the classical approach, based on the hypergeometric distribution, called “enrichment analysis” or “GO over-representation analysis” in the terminology of clusterprofiler documentation.

Représentation simplifiée graphique du test hypergéometrique

We can easily extract the IDs of the genes in cluster2 using for example:

• names(tss)[tss$clusters == "cluster2"]

The main question that we need to answer is what is the gene universe? In general, you should think as the “gene universe” as the set of genes that could potentially have been selected in your gene list. Since we did different filters, e.g. by selecting genes present in the RNAseq dataset or that had a match between the ENSEMBLID and the ENTREZID, our gene universe does NOT encompass all the genes in the genome. In our case, a good definition of the gene universe would be all the genes in the tss object, since we performed kmeans on these genes to identify genes in cluster 2.

So now, let’s use this info to compute the GO over-representation statistics:

# load the library
library(clusterProfiler)
# Run the enrichment analysis
ego <- enrichGO(gene          = names(tss)[tss$clusters == "cluster2"],
                universe      = names(tss),
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                readable      = TRUE)

And take a look at the results:

head(ego)

##                    ID                                     Description GeneRatio
## GO:0042438 GO:0042438                    melanin biosynthetic process   10/1041
## GO:0046189 GO:0046189 phenol-containing compound biosynthetic process   12/1041
## GO:0006582 GO:0006582                       melanin metabolic process   10/1041
## GO:0044550 GO:0044550       secondary metabolite biosynthetic process   10/1041
## GO:0002385 GO:0002385                         mucosal immune response    8/1041
## GO:0032922 GO:0032922         circadian regulation of gene expression   17/1041
##             BgRatio RichFactor FoldEnrichment   zScore       pvalue    p.adjust
## GO:0042438 21/11891  0.4761905       5.439367 6.306766 3.642124e-06 0.005906749
## GO:0046189 31/11891  0.3870968       4.421679 5.908526 5.578739e-06 0.005906749
## GO:0006582 22/11891  0.4545455       5.192123 6.095926 6.154427e-06 0.005906749
## GO:0044550 22/11891  0.4545455       5.192123 6.095926 6.154427e-06 0.005906749
## GO:0002385 14/11891  0.5714286       6.527240 6.409445 6.247561e-06 0.005906749
## GO:0032922 59/11891  0.2881356       3.291278 5.464823 7.365024e-06 0.005906749
##                 qvalue
## GO:0042438 0.005694326
## GO:0046189 0.005694326
## GO:0006582 0.005694326
## GO:0044550 0.005694326
## GO:0002385 0.005694326
## GO:0032922 0.005694326
##                                                                                                            geneID
## GO:0042438                                                PMEL/TYRP1/DCT/TYR/MFSD12/RAB38/ZEB2/GIPC1/SLC7A11/MC1R
## GO:0046189                                     PMEL/TYRP1/DCT/TYR/MFSD12/RAB38/ZEB2/SNCA/GIPC1/SLC7A11/NR4A2/MC1R
## GO:0006582                                                PMEL/TYRP1/DCT/TYR/MFSD12/RAB38/ZEB2/GIPC1/SLC7A11/MC1R
## GO:0044550                                                PMEL/TYRP1/DCT/TYR/MFSD12/RAB38/ZEB2/GIPC1/SLC7A11/MC1R
## GO:0002385                                                    OTUD7B/RAB17/H2BC21/H2BC8/H2BC12/H2BC11/H2BC4/H2BC7
## GO:0032922 ATF4/PPP1CB/PPP1CC/HDAC2/BHLHE41/CSNK1E/PER1/TOP1/GFPT1/BHLHE40/RAI1/MYCBP2/NPAS2/NCOA2/ID2/NOCT/NR1D1
##            Count
## GO:0042438    10
## GO:0046189    12
## GO:0006582    10
## GO:0044550    10
## GO:0002385     8
## GO:0032922    17

6.2 KEGG enrichment analysis in DE genes from the different clusters

Now let’s perform the same type of analysis but focusing on the differentially expressed (DE) genes present in each of the 4 clusters. Here, we’re asking if DE genes that have different BRG1 binding patterns at their TSS are enriched in specific KEGG pathways.
The clusterProfiler package provides the function compareCluster that allows to analyze and compare enrichment in different group of genes. This function recognizes ENTREZ identifiers.

# Select genes whose gene expression is significantly changing
# between conditions
geneList <- tss[tss$signif != "stable",]

# The function split can split a vector to a list following categories
# in an other vector (take a look at the geneList object after running split).
geneList <- split(names(geneList), geneList$clusters)

# Now perform a KEGG pathway over-representation analysis 
compKEGG <- compareCluster(geneCluster   = geneList,
                           fun           = "enrichKEGG",
                           pvalueCutoff  = 0.05,
                           pAdjustMethod = "BH")
# Plot the results
dotplot(compKEGG, 
        showCategory = 15, 
        title = "KEGG Pathway Enrichment Analysis")

Let’s focus on cluster 4 and visualize which genes of the the KEGG pathway “Apoptosis” are both differentially expressed and bound by BRG1 and MITF.
We obtain the identifier of the pathway by looking at the compKEGG object:

• compKEGG@compareClusterResult

So the ID of the KEGG pathway “Apoptosis” is hsa04210.
We use the pathview package which maps differential gene expression values to KEGG pathway maps and create a PNG file in your current directory that you get with getwd().

# load the library
library(pathview)

# Get the names of the genes in cluster 4
EntrezID.cl4 <- names(tss)[tss$clusters == "cluster4"]

# match these EntrezIDs to the RNAseq dataset
mmid <- match(EntrezID.cl4, RNAseq$ENTREZID)

# Get the log fold-change for these genes from the RNAseq table
logFC.cl4 <- RNAseq$logFC[mmid]
names(logFC.cl4) <- RNAseq$ENTREZID[mmid]

# Create an image of the pathway with colors for differential gene expression
pathview(logFC.cl4, pathway.id = "hsa04210", species = "hsa")

Pathview

7.1 De novo motif discovery

Here, we will search for motifs enriched under the SOX10 peaks. Since we do not start from a pre-defined list of motifs, e.g. available in a database like JASPAR or CIS-BP, this is called a de novo search. We start by sorting SOX10 peaks by qvalue and keep the sequences corresponding to the peak centers +/- 150bp:

sox10peaks_300bp <- resize(peakfilt[["SOX10"]][order(peakfilt[["SOX10"]]$qvalue, 
                                                     decreasing = TRUE)],
                           width = 300,
                           fix = "center")

Next, we extract the underlying sequences and obtain a DNAStringSet object:

sox10_seq <- getSeq(Hsapiens, sox10peaks_300bp)

Most of the time, you will want to export these sequences as a fasta file and run an external motif discovery program such as:

• MEME or
• RSAT

Note that the TFBStools package has a runMEME function that is a wrapper to MEME.
To export your the sequences you can use:

## Define the path to the fasta file:
myfastafile = file.path("data", "sox10peaks_300bp.fa")
## NOT RUN (already exported): export sequences in fasta format
writeXStringSet(sox10_seq, myfastafile)

Here, we will use the BCRank package which can also perform a de novo search for motifs within R. The method takes a list of DNA sequences sorted by some measure of the probability that a relevant motif is present. Here, the qvalue plays the role of this measure, which is why we sorted the peaks by qvalue. BCrank takes a long time to run, so we will not run the commands below.

##NOT RUN (too long...)
set.seed(123) #for reproducibility
BCRANKout <- bcrank(myfastafile,
                    restarts = 25, 
                    silent = TRUE,
                    use.P1 = TRUE,
                    use.P2 = TRUE)
saveRDS(BCRANKout, file.path("data", "BCRANKout.rds"))

Instead, we ran the commands for you and saved the object as an rds file. So we can just import this file and take a look at it:

BCRANKout <- readRDS(file.path("data", "BCRANKout.rds"))
BCRANKout

## 
## An object of class "BCRANKresult"
## 
## Top 25 DNA motifs predicted by BCRANK:
## 
##       Consensus     Score
## 1       GCWTTGT 122.34709
## 2  VBRVCCTGCTTT  84.48699
## 3   GSRKCANTGTT  75.15736
## 4     VCTGTGAGB  71.66138
## 5   ATKGCACCWDD  69.77417
## 6    GVVTTGRGAB  69.57777
## 7     TYCAMAKCC  65.85154
## 8     TTSCCVAAG  65.51148
## 9   ASWTTTTGCHC  65.39492
## 10   SNCAASHCGG  64.43600
## 11   MCGGVYCAGC  64.12572
## 12  ACAAHKCAVGY  62.69532
## 13   TGGGGCSDDG  62.05447
## 14 TTTTGDAGSTCT  61.36914
## 15     TCGAKASC  59.74472
## 16  WTBCAATRATA  59.70302
## 17     CRKGTCDG  59.32284
## 18  AHAATDBGRCM  58.34722
## 19  SCYTCACBNGY  58.12427
## 20  CYAKTGNNVAC  56.36697
## 21   DTGTVHATTY  56.26696
## 22 HHTARTCTGTCY  55.40342
## 23   BCMVCGNGCA  55.34975
## 24  KCWSWTCCBTT  51.94605
## 25   CBCNVCGAGD  51.21762
## 
## 
## Use methods toptable(object) and fname(object) to access object slots.

We can extract the top motif found by BCRANK with:

topMotif <- toptable(BCRANKout, 1)

and extract the Position Weight Matrix (PWM) for this motif and plot the sequence logo with:

pwmnorm_topmotif <- BCRANK::pwm(topMotif, normalize = TRUE)  #with normalize=TRUE the columns will sum to 1 (PFM)
library(seqLogo)
seqLogo(reverseComplement(pwmnorm_topmotif))

You can go to the JASPAR website and search “SOX10”. You will see that the first motif found by BCRANK is very similar to the motifs that have been found for SOX10. We will actually retrieve and take a look at one of these motifs below.

7.2 Searching for known motifs

Another approach consists in searching for known motifs in the ChIP-seq peaks. Since the SOX10 has been extensively studied, there are already SOX10 binding motifs available in the databases.

library(MotifDb)

Search for SOX10 motifs in JASPAR databases

sox10_motifs = MotifDb::query(MotifDb, c("hsapiens", "jaspar", "sox10"))
## We keep the motif from JASPAR 2022
Sox10JASP22 <- sox10_motifs[["Hsapiens-jaspar2022-SOX10-MA0442.2"]]
seqLogo(Sox10JASP22)

MotifDb returns a position frequency matrix (PFM) with columns summing to 1. The number of sequences used to build the PFM are available in the elementMetadata or mcols of the object. We use this information to obtain a PWM:

nseq <- as.integer(mcols(sox10_motifs["Hsapiens-jaspar2022-SOX10-MA0442.2"])$sequenceCount)
sox10pfm <- apply(round(nseq * Sox10JASP22,0), 2, as.integer)
rownames(sox10pfm) <- rownames(Sox10JASP22)
Sox10JASP22_PWM <- PWM(sox10pfm)

And then use this PWM to screen the whole genome. Again, this takes some time so we saved the rds object for you

## NOT RUN (too long...)
SOX10_hg38 <- matchPWM(Sox10JASP22_PWM, 
                       Hsapiens, 
                       min.score="90%")
saveRDS(SOX10_hg38, file.path("data", "SOX10_hg38.rds"))

We obtain a GRanges with the location of all the motifs:

SOX10_hg38 <- readRDS(file.path("data", "SOX10_hg38.rds"))
SOX10_hg38

## GRanges object with 3127286 ranges and 2 metadata columns:
##                        seqnames        ranges strand |     score         string
##                           <Rle>     <IRanges>  <Rle> | <numeric> <DNAStringSet>
##         [1]                chr1   11631-11641      + |  0.900030    TCCACAAAGTG
##         [2]                chr1   13572-13582      + |  0.911933    GCTACAAAGGT
##         [3]                chr1   16770-16780      + |  0.935700    AACACAAAGTG
##         [4]                chr1   19982-19992      + |  0.948358    GGGACAAAGGC
##         [5]                chr1   20010-20020      + |  0.930367    CCCACAAAGGA
##         ...                 ...           ...    ... .       ...            ...
##   [3127282] chrX_MU273397v1_alt 314005-314015      - |  0.929517    TGGACAAAGCT
##   [3127283] chrX_MU273397v1_alt 319703-319713      - |  0.952440    AATACAAAGGG
##   [3127284] chrX_MU273397v1_alt 319894-319904      - |  0.909040    AAAACAATGAA
##   [3127285] chrX_MU273397v1_alt 324082-324092      - |  0.955583    CCAACAAAGGC
##   [3127286] chrX_MU273397v1_alt 327059-327069      - |  0.960064    TACACAAAGAA
##   -------
##   seqinfo: 711 sequences (1 circular) from hg38 genome

We can now overlap these motifs with our ChIPseq peaks:

## We set a minimum overlap of 11bp (size of the motif) to count only
## the motifs that are fully within the peaks
count_sox10 <- countOverlaps(sox10peaks_300bp,
                             SOX10_hg38,
                             minoverlap = 11)
## Number of peaks that have 0, 1, 2, etc. SOX10 binding sites
table(count_sox10)

## count_sox10
##    0    1    2    3    4    5 
## 1734 1724  611  123   20    2

## Percentage of our peaks that have at least one SOX10 binding site:
100*mean(count_sox10>=1)

## [1] 58.85145

We see that more than 58.9% of the peaks contain at least one of these motifs.
Assessing whether our peaks are significantly enriched for SOX10 motifs would probably require to extract random sequences that share common characteristics with SOX10 ChIP-seq peaks in terms of size and genomic distribution and evaluate their enrichment in SOX10 binding sites.
Here as a substitute, let’s evaluate the frequency of SOX10 peaks in a random selection of promoters

set.seed(123) # for reproducibility
count_sox10_tss <- countOverlaps(promoters(genes(txdb), upstream=300, downstream=0)[sample.int(length(sox10peaks_300bp))],
                                 SOX10_hg38,
                                 minoverlap = 11)

##   2162 genes were dropped because they have exons located on both strands
##   of the same reference sequence or on more than one reference sequence,
##   so cannot be represented by a single genomic range.
##   Use 'single.strand.genes.only=FALSE' to get all the genes in a
##   GRangesList object, or use suppressMessages() to suppress this message.

## Percentage of 300bp promoters that contain SOX10 peaks:
100*mean(count_sox10_tss>=1)

## [1] 10.08543

This suggests that the enrichment of SOX10 binding sites in our ChIP-seq peaks is very significant.

8.1 Mining GEO ChIP-seq data with ChIPSeeker

There Gene Expression Omnibus (GEO) contains a vast number of ChIP-seq datasets.
We can compare our own dataset to those deposited in GEO to search for significant overlap data. Significant overlap of ChIP seq data by different binding proteins may be used to infer cooperative regulation and thus can be used to generate hypotheses.

The ChIPseeker package (Wang et al. 2022) package provides useful functions to download and use GEO datasets in R / Bioconductor.

We can collect >100,000 bed files deposited in GEO. You can use getGEOspecies to get the number of files available for different species.

getGEOspecies()[1:25,]

##                                                species Freq
## 1                      Actinobacillus pleuropneumoniae    1
## 2  Actinobacillus pleuropneumoniae serovar 3 str. JL03    1
## 3                                        Aedes aegypti   11
## 4                                             Anabaena    6
## 5                                  Anolis carolinensis    7
## 6                                    Anopheles gambiae    2
## 7                                       Apis mellifera   17
## 8                            Apis mellifera scutellata    1
## 9                                          Arabidopsis   17
## 10                                  Arabidopsis lyrata    4
## 11                                Arabidopsis thaliana 2473
## 12                                Atelerix albiventris    2
## 13                                   Bombus terrestris    8
## 14                                          Bos taurus  257
## 15                             Brachypodium distachyon   10
## 16                           Branchiostoma lanceolatum   93
## 17                                       Brassica rapa   12
## 18                              Caenorhabditis elegans  532
## 19                                  Callithrix jacchus   48
## 20                                    Candida albicans   34
## 21                                Candida dubliniensis   20
## 22                                         Canis lupus    1
## 23                              Canis lupus familiaris   35
## 24                               Capsaspora owczarzaki   21
## 25                                       Carica papaya    1

You can also use the getGEOgenomeVersion to get the same information but based on genome version used to analyze the data.

You can access the detailed information using the getGEOInfo function, for each genome version.

hg38 <- getGEOInfo(genome = "hg38", simplify = TRUE)
head(hg38)

##       series_id        gsm     organism                                 title
## 16488  GSE58207 GSM1403308 Homo sapiens         Lactimidomycin treated HCT116
## 16489  GSE58207 GSM1403308 Homo sapiens         Lactimidomycin treated HCT116
## 16490  GSE58207 GSM1403307 Homo sapiens          Cycloheximide treated HCT116
## 16491  GSE58207 GSM1403307 Homo sapiens          Cycloheximide treated HCT116
## 20345  GSE67978 GSM1660032 Homo sapiens   H3K27ac_Human_Brain_WhiteMatter_HS2
## 20351  GSE67978 GSM1660029 Homo sapiens H3K27ac_Human_Brain_OccipitalPole_HS2
##                                                                                                                                 supplementary_file
## 16488                                   ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403308/suppl/GSM1403308_HCT116_LTM_sense.bedGraph.gz
## 16489                              ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403308/suppl/GSM1403308_HCT116_LTM_anti-sense.bedGraph.gz
## 16490                              ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403307/suppl/GSM1403307_HCT116_CHX_anti-sense.bedGraph.gz
## 16491                                   ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403307/suppl/GSM1403307_HCT116_CHX_sense.bedGraph.gz
## 20345   ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1660nnn/GSM1660032/suppl/GSM1660032_H3K27ac_Human_Brain_WhiteMatter_HS2_hg38_peaks.narrowPeak.gz
## 20351 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1660nnn/GSM1660029/suppl/GSM1660029_H3K27ac_Human_Brain_OccipitalPole_HS2_hg38_peaks.narrowPeak.gz
##       genomeVersion pubmed_id
## 16488          hg38      <NA>
## 16489          hg38      <NA>
## 16490          hg38      <NA>
## 16491          hg38      <NA>
## 20345          hg38      <NA>
## 20351          hg38      <NA>

ChIPseeker provides the function downloadGEObedFiles to download all the bed files of a particular genome.

downloadGEObedFiles(genome = "hg38", destDir = "hg38")

Or a vector of GSM accession number by downloadGSMbedFiles.

gsm <- hg38$gsm[sample(nrow(hg38), 10)]
downloadGSMbedFiles(gsm, destDir = "hg38")

After the download of bed files from GEO, we can pass them to enrichPeakOverlap for testing the significance of the overlaps. Parameter targetPeak can be the folder, e.g. hg38, that containing bed files. enrichPeakOverlap will parse the folder and compare all the bed files. It is possible to test the overlap with bed files that are mapped to different genome or different genome versions: enrichPeakOverlap provides a parameter chainFile that can pass a chain file and liftOver the targetPeak to the genome version consistent with queryPeak. Significant overlaps can be used to generate hypotheses on cooperative regulation. By mining the data deposited in GEO, we can identify putative complexes or interacting regulators of gene expression or chromatin remodeling for further validation.

8.2 The AnnotationHub

The AnnotationHub is a web resource that provides a central location where genomic files in different formats (e.g. vcf, bed, wig, etc.) and other resources coming from different locations (e.g. UCSC, ENSEMBL, ENCODE, etc.) can be mined and downloaded to enrich your analyses. Importantly, the resource that are available include metadata with e.g. description of the experiments, tags and date of modification which are important to figure out where exactly the data comes from and if it fits with the aim of your analyses.
The AnnotationHub package is a client for R/Bioconductor to search these resources and download the data in appropriate formats for their analysis in the R/Bioconductor environment.
The client creates and manages a local cache of files that are retrieved by the user, which allows quick and reproducible access to the data.

First we need to create an AnnotationHub instance. This will create a local cache on your machine so that repeat queries for the same information (even in different R sessions) will be very fast. By default, the cache will be created in your home folder, which may not be appropriate in a shared server if you have a limited quota of space in your home directory. You can change the location of the cache with:

## NOT RUN: example of how to set the path to the cache
setAnnotationHubOption("CACHE") <- "/shared/projects/projectname/AnnotationHub"
## to check that the path was set correctly:
getAnnotationHubOption("CACHE")

Then we create the AnnotationHub instance with:

##NOT RUN (too long...)
ah <- AnnotationHub()
# |=============================================================================| 100%
#
# snapshotDate(): 2024-10-28

This will take some time the first time you run it since there are a lot of resources available. The ah object is quite big but it only contains pointers to online data, not the data itself which would be way too big to download all at once. The information content is changing frequently, which is why there is a snapshotDate.

Let’s take a look at the ah object:

##NOT RUN
ah
# AnnotationHub with 72100 records
# # snapshotDate(): 2024-10-28
# # $dataprovider: Ensembl, BroadInstitute, UCSC, ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, FANTOM5,DLRP,IUPHAR,HPRD,STRING,SWISSPROT,TREMBL,ENSEMBL,CELLPHONEDB,...
# # $species: Homo sapiens, Mus musculus, Drosophila melanogaster, Rattus norvegicus, Bos taurus, Pan troglodytes, Danio rerio, Gallus gallus, Sus scrofa, Mon...
# # $rdataclass: GRanges, TwoBitFile, BigWigFile, EnsDb, Rle, OrgDb, SQLiteFile, ChainFile, TxDb, Inparanoid8Db
# # additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags, rdatapath,
# #   sourceurl, sourcetype 
# # retrieve records with, e.g., 'object[["AH5012"]]' 
# 
#              title                                             
#   AH5012   | Chromosome Band                                   
#   AH5013   | STS Markers                                       
#   AH5014   | FISH Clones                                       
#   AH5015   | Recomb Rate                                       
#   AH5016   | ENCODE Pilot                                      
#   ...        ...                                               
#   AH119506 | Ensembl 113 EnsDb for Zonotrichia albicollis      
#   AH119507 | Ensembl 113 EnsDb for Zalophus californianus      
#   AH119508 | Ensembl 113 EnsDb for Zosterops lateralis melanops
#   AH119518 | MassBank CompDb for release 2024.06               
#   AH119519 | MassBank CompDb for release 2024.11             

As of 2024-10-28, the AnnotationHub contains 72,100 records.
The ah object is organized as a vector (so its length is 72100 in my case). You can get information about a single resource by indexing with a single [; using two brackets ([[) downloads the object:

##NOT RUN
ah[1]
# AnnotationHub with 1 record
# # snapshotDate(): 2024-10-28
# # names(): AH5012
# # $dataprovider: UCSC
# # $species: Homo sapiens
# # $rdataclass: GRanges
# # $rdatadateadded: 2013-03-26
# # $title: Chromosome Band
# # $description: GRanges object from UCSC track 'Chromosome Band'
# # $taxonomyid: 9606
# # $genome: hg19
# # $sourcetype: UCSC track
# # $sourceurl: rtracklayer://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/cytoBand
# # $sourcesize: NA
# # $tags: c("cytoBand", "UCSC", "track", "Gene", "Transcript", "Annotation") 
# # retrieve record with 'object[["AH5012"]]' 

Now, we can use various tools to narrow down the dataset or datasets that are actually of interest for our analyses.
Let’s say that we are looking for a dataset describing regions in the human genome that are potential distal enhancers. The ENCODE project has typically defined such regions using the data they have produced on e.g. histone marks H3K27ac, H3K4me3 or DNAseI accessibility.

Let’s take a look at the data providers:

##NOT RUN
unique(ah$dataprovider)
#  [1] "UCSC"                                                                                                      
#  [2] "Ensembl"                                                                                                   
#  [3] "RefNet"                                                                                                    
#  [4] "Inparanoid8"                                                                                               
# ...
# [72] "ENCODE" 
# ...
# [83] "MGI"                                                                                                       
# [84] "ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/" 

Let’s take a look at the available data from encode:

##NOT RUN
## Number of datasets available from ENCODE (note that there are more since in 
## "dataprovider" there is also e.g. "ENCODE Project" or "ENCODE SCREEN v3"
sum(ah$dataprovider == "ENCODE")
# [1] 20
# So there are 20 datasets from the dataprovider "ENCODE"

## Take a look at the titles of these datasets
ah[ah$dataprovider == "ENCODE"]$title
#  [1] "hg38.Kundaje.GRCh38_unified_Excludable"                   "hg38.Bernstein.Mint_Excludable_GRCh38"                   
#  [3] "hg38.Kundaje.GRCh38.Excludable"                           "hg38.Reddy.wgEncodeDacMapabilityConsensusExcludable.hg38"
#  [5] "hg38.Wold.hg38mitoExcludable"                             "hg38.Yeo.eCLIP_Excludableregions.hg38liftover.bed.fixed" 
#  [7] "hg19.Bernstein.Mint_Excludable_hg19"                      "hg19.Birney.wgEncodeDacMapabilityConsensusExcludable"    
#  [9] "hg19.Crawford.wgEncodeDukeMapabilityRegionsExcludable"    "hg19.Wold.hg19mitoExcludable"                            
# [11] "hg19.Yeo.eCLIP_Excludableregions.hg19"                    "mm10.Hardison.Excludable.full"                           
# [13] "mm10.Hardison.psuExcludable.mm10"                         "mm10.Kundaje.anshul.Excludable.mm10"                     
# [15] "mm10.Kundaje.mm10.Excludable"                             "mm10.Wold.mm10mitoExcludable"                            
# [17] "mm9.Wold.mm9mitoExcludable"                               "ENCODE_dELS_regions"                                     
# [19] "ENCODE_pELS_regions"                                      "ENCODE_PLS_regions"  

I wonder what “PLS_regions” are 🤔 and I want to know more. So let’s take a look at the object description:

##NOT RUN
ah[ah$dataprovider == "ENCODE" & ah$title == "ENCODE_PLS_regions"]$description
# [1] "A GRanges object containing regions of candidate cis-regulatory elements with promoter-like signatures 
# as identified by the ENCODE SCREEN project. These consist of regions with high H3K4me3 and DNase signal, and 
# located within 200 bp of a GENCODEtranscription start site."

OK, this looks very close to what we are looking for and I suppose that pELS and dELS refer to proximal and distal enhancer-like signatures. Let’s verify:

##NOT RUN
ah[ah$dataprovider == "ENCODE" & ah$title == "ENCODE_dELS_regions"]$description
# [1] "A GRanges object containing regions of candidate cis-regulatory elements with distal enhancer-like
# signatures as identified by the ENCODE SCREEN project. These consist of regions with high H3K27ac and 
# DNase signal, but low H3K4me3 signal, and located more than 2kb from GENCODE transcription start sites."

So now we can get the names of the object, download it and use it in R (it’s a GRanges object):

##NOT RUN
## Get the name / ID of the dataset
names(ah[ah$dataprovider == "ENCODE" & ah$title == "ENCODE_dELS_regions"])
# [1] "AH116727"

## Download the dataset (use double brackets `[[`)
dELS <- ah[["AH116727"]]
# downloading 1 resources
# retrieving 1 resource
#   |===========================================================================| 100%
# 
# loading from cache
# require(“GenomicRanges”)

## Take a look at the object:
dELS
# GRanges object with 786756 ranges and 0 metadata columns:
#            seqnames            ranges strand
#               <Rle>         <IRanges>  <Rle>
#        [1]     chr1     271227-271468      *
#        [2]     chr1     274330-274481      *
#        [3]     chr1     605331-605668      *
#        [4]     chr1     727122-727350      *
#        [5]     chr1     807737-807916      *
#        ...      ...               ...    ...
#   [786752]     chrY 26319537-26319816      *
#   [786753]     chrY 26319848-26320155      *
#   [786754]     chrY 26363541-26363721      *
#   [786755]     chrY 26670949-26671287      *
#   [786756]     chrY 26671293-26671478      *
#   -------
#   seqinfo: 24 sequences from an unspecified genome; no seqlengths

## Since seqinfo has not been set, you may wonder which genome version was used.
## You can get this info from the metadata of the ah object:
ah['AH116727']$genome
# [1] "hg38"

In another example, we will use some useful function to search the AnnotationHub:

• subset (or [) to do a specific subseting operation.
• query to do a command-line search over the metadata of the hub.
• display to get a Shiny interface in a browser, so you can browse the object.

We would like to know if there are any other dataset available for MITF in human. We first filter the hub to keep only the datasets obtained in human, or the datasets where the species is unspecified (NA), in case some metadata is missing.

##NOT RUN
ah <- subset(ah, species %in% c(NA, "Homo sapiens"))

Then we query the metadata for MITF:

##NOT RUN
query(ah, "MITF")
# AnnotationHub with 1 record
# # snapshotDate(): 2024-10-28
# # names(): AH22280
# # $dataprovider: Pazar
# # $species: NA
# # $rdataclass: GRanges
# # $rdatadateadded: 2015-03-09
# # $title: pazar_MITF_Strub_20120522.csv
# # $description: TF - Target Gene file from pazar_MITF_Strub_20120522
# # $taxonomyid: NA
# # $genome: NA
# # $sourcetype: CSV
# # $sourceurl: http://www.pazar.info/tftargets/pazar_MITF_Strub_20120522.csv
# # $sourcesize: 34964
# # $tags: c("Pazar", "Transcription Factors") 
# # retrieve record with 'object[["AH22280"]]' 

Indeed this dataset does not contain all the metadata but there is a link to the original source file and searching the web for “MITF, Strub” quickly gets you to this publication where all the info is available to evaluate if this dataset is relevant for us.

Finally, another way to mine the AnnotationHub is to use a browser with:

##NOT RUN
outdata <- display(ah)

Note how we assign the call to the display function to an R object (oudata) so that we can actually capture the selection of datasets that we will make from the browser into our R session.