\(^1\) Université de Bordeaux, INSERM U1218

\(^2\) GenomEast platform, IGBMC

1 Introduction

During this training session, we are going to use ChIP-seq and RNA-seq data from Laurette et al. (Laurette et al. 2015) which are available in GEO as GSE61967. Data have already been processed.

1.1 Processing of RNA-seq data

1.1.1 Preprocessing

Reads were preprocessed in order to remove adapter, polyA and low-quality sequences (Phred quality score below 20). After this preprocessing, reads shorter than 40 bases were discarded for further analysis. These preprocessing steps were performed using cutadapt (Martin 2011) version 1.10.

1.1.2 Mapping

Reads were mapped onto the hg38 assembly of Homo sapiens genome using STAR (Dobin et al. 2013) version 2.5.3a.

1.1.3 Quantification

Gene expression quantification was performed from uniquely aligned reads using htseq-count version 0.6.1.p1 (S. Anders, Pyl, and Huber 2015), with annotations from Ensembl version 103 and “union” mode. Only non-ambiguously assigned reads have been retained for further analyses.

1.1.4 Normalization

Read counts have been normalized across samples with the median-of-ratios method proposed by Anders and Huber (Anders and Huber 2010), to make these counts comparable between samples.

1.1.5 Differential expression analysis

Differential expressed genes between wild type ans shBRG1 cells was performed using the edgeR R package (Robinson, McCarthy, and Smyth 2010). We called significant changes when FDR < 0.01, absolute log fold change over 1 and minimum average log normalized count over 5.

1.2 Processing of ChIP-seq data

1.2.1 Mapping

Reads were mapped to Homo sapiens genome (assembly hg38) using Bowtie (Langmead et al. 2009) v1.0.0 with default parameters except for “-p 3 -m 1 –strata –best –chunkmbs 128.” The following table shows the number of reads aligned to the Homo sapiens genome.

Mapping statistics of ChIP-seq data analyzed during this training session. Column “Raw reads” corresponds to the number of input reads. Column “Aligned reads” corresponds to the number of reads aligned exactly 1 time. Column “Multimapped” corresponds to the number of reads aligned > 1 times. Column “Unmapped” corresponds to the number of reads aligned 0 time.
Sample ID Sample name raw reads aligned reads multimapped reads unmapped reads
SRR1594290 MITF 84856252 54369181 (64.07%) 12823841 (15.11%) 17663230 (20.82%)
SRR1594291 SOX10 38202152 29594694 (77.47%) 6128890 (16.04%) 2478568 (6.49%)
SRR1594292 BRG1 50190992 40134919 (79.96%) 6403809 (12.76%) 3652264 (7.28%)
SRR1594294 BRG1_CTRL 38148419 26273410 (68.87%) 5722823 (15.00%) 6152186 (16.13%)
SRR1594295 MITF_Input 29433042 19970925 (67.85%) 4049711 (13.76%) 5412406 (18.39%)
SRR1594296 SOX10_Input 35449561 27381422 (77.24%) 7045192 (19.87%) 1022947 (2.89%)

1.2.2 Peak Calling

Prior to peak calling, reads falling into Encode blacklisted regions (“(2014) Mod/Mouse/humanENCODE: Blacklisted Genomic Regions for Functional Genomics Analysis - Anshul Kundaje” n.d.) were removed using bedtools intersect v2.26.0 (Quinlan and Hall 2010). Then peak calling was done with Macs2 v2.1.1 with default parameters.

Number of peaks detected
IP sample Input sample No. of peaks
BRG1 BRG1_CTRL 72024
MITF MITF_Input 9702
SOX10 SOX10_Input 4538

1.2.3 Generation of BigWig files

Normalized BigWig files were generated using Homer (Heinz et al. 2010) makeUCSCfile v4.11.0 with the following parameter ’-norm 1e7’ meaning that data were normalized to 10M reads.

1.3 Define the paths to the data

We create a list containing the absolute paths to the data we will use along the practice and their file names.

params <- list(
  chipseq= "/shared/projects/ebai2021_n2/data/chip_seq/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= "/shared/projects/ebai2021_n2/data/chip_seq/RNAseq/RNAseq_diff_norm.RData"
  )

2 Description and statistics of BRG1 dataset

2.1 Load data

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

  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.
  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 no name 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.
## The package ChIPseeker provides a function to load peak files such as narrowPeaks as GRanges objects
## Here BRG1 peak set is loaded into a list of peaks
## this list can be extended if there are more datasets
library(ChIPseeker)
peaks <- list()
peaks[["BRG1"]] <- readPeakFile(paste0(params$chipseq, params$brg1), as="GRanges")
peaks
## $BRG1
## GRanges object with 72024 ranges and 7 metadata columns:
##           seqnames            ranges strand |                     V4        V5
##              <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
##                    V6        V7        V8        V9       V10
##           <character> <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 objects; this is an R format which look like the bed format, but is optimized in terms of memory requirements and speed of execution.

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 use these functions all along the course.

# sapply() function takes list, vector or data frame as input and gives output in vector or matrix
# sapply apply the same function (here length) to all elements
# of the list "peaks"
sapply(peaks,length)
##  BRG1 
## 72024

Make a simple barplot showing the number of BRG1 peaks.

Show code: barplot 
barplot(sapply(peaks,length))

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

Show code: load ggplot2 
# 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.
Show code: create the data.frame 
# create a table with the data to display
peak.lengths <- data.frame(IP=names(peaks),
                           NbPeaks=sapply(peaks,length))
  • use ggplot2 with the appropriate geometric object geom_*
Show code: barplot with ggplot2 
# make the barplot
ggplot(peak.lengths, aes(x=IP, y=NbPeaks)) +
         geom_col()

We can customize the plots by changing colors.

# Let's add colors to the barplot
# In R it exists some already defined colors palettes
# the most widely used palette is RColorBrewer.
# This R library offers several color palettes
# See:
library(RColorBrewer)
par(mar=c(3,4,2,2))
display.brewer.all()

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

Show code: color according to the chipped TF (column IP) 
ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
         geom_col()

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.

Show code: color according to the chipped TF (column IP) and set RColorBrewer palette to Set1 
ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
         geom_col()+
         scale_fill_brewer(palette="Set1")

2.3 How large are these peaks?

BRG1 peaks were loaded as a GRanges object.
- To manipulate it we load the library GenomicRanges - Use the appropriate function to retrieve the width of the peaks from the GRanges - Use a function to summarize the peak length (display the statistics (minimum, maximum, quartiles) of the distribution)

Show code: Display peak length statistics 
## 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
# or quantile(width(peaks$BRG1))

Create a simple boxplot of the peak sizes.

Show code: Display peak length statistics as a boxplot 
peak.width <- lapply(peaks,width)
boxplot(peak.width)

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:

# 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_*.

Show code: Display peak length statistics with ggplot2 
## 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 axis to log scale. Finally, we want to remove the x label, change the y label to ‘Peak sizes’ and the legend relative to the fill color from “L1” to “TF.”
Using the ggplot2 documentation (e.g. https://ggplot2.tidyverse.org/ or google it, try to enhance the boxplot.

Enhance it ! 
# - 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 >= 8. The qValue corresponds to the 9th column of narrowPeak files. So, we are going to set a threshold on this. How many peaks are selected ?

Show the code: select high quality peaks 
## Select high quality peaks
peaks$BRG1 <- peaks$BRG1[peaks$BRG1$V9 >= 8]
## Compute the number of remaining peaks
length(peaks$BRG1)
## [1] 30874

2.5 Where are the peaks located over the whole genome?

Sometime, peaks may occur more in some chromosomes than others. 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.

# genome wide BRG1 peak distribution
covplot(peaks$BRG1, weightCol="V5")

# chromosome wide BRG1 peak distribution
covplot(peaks$BRG1, chrs=c("chr1", "chr2"), weightCol="V5")

2.6 Functional annotation: genomic features enriched in BRG1 peaks

2.6.1 Generate annotation

We can assign peaks to the closest genes and genomic features (introns, exons, promoters, distal regions, etc…). We load ready-to-use annotation objects from Bioconductor: a OrgDB org.Hs.eg.db and a TxDB TxDb.Hsapiens.UCSC.hg38.knownGene objects.

## 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 annotation
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

This is done using the function annotatePeak which compares peak positions with the genomic feature positions of the reference genome. This function returns a complex object which contains all this information. As for the peaks coordinates, we store the peak annotations in a list that we first initialize as an empty list peakAnno <- list(). We define here the TSS regions as TSS positions -1000 bp +100 bp. Using the help of the function and annotation objects we have loaded, try to write the correct command line to annotate BRG1 peaks.

Show code: annotate BRG1 peaks 
## 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(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")
## >> preparing features information...      2021-05-25 19:09:26 
## >> identifying nearest features...        2021-05-25 19:09:28 
## >> calculating distance from peak to TSS...   2021-05-25 19:09:30 
## >> assigning genomic annotation...        2021-05-25 19:09:30 
## >> adding gene annotation...          2021-05-25 19:10:13 
## >> assigning chromosome lengths           2021-05-25 19:10:13 
## >> done...                    2021-05-25 19:10:13
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                          V4  V5 V6
## 1     chr1  983476  984526  1051      *  BRG1siCTRL_CHIP-seq_peak_2 108  .
## 2     chr1 1000729 1001179   451      *  BRG1siCTRL_CHIP-seq_peak_3  90  .
## 3     chr1 1020914 1021207   294      *  BRG1siCTRL_CHIP-seq_peak_5 114  .
## 4     chr1 1069156 1070192  1037      * BRG1siCTRL_CHIP-seq_peak_10  96  .
## 5     chr1 1079517 1080087   571      * BRG1siCTRL_CHIP-seq_peak_11 179  .
## 6     chr1 1304853 1306013  1161      * BRG1siCTRL_CHIP-seq_peak_15 108  .
##        V7       V8       V9 V10
## 1 5.97792 13.59661 10.83201 384
## 2 6.65677 11.62894  9.01201 199
## 3 8.11634 14.29797 11.46427 158
## 4 6.90466 12.34852  9.65604 410
## 5 9.95583 21.25171 17.96153 213
## 6 5.97792 13.59661 10.83201 176
##                                          annotation geneChr geneStart geneEnd
## 1                                 Distal Intergenic       1    975204  982093
## 2                                          Promoter       1   1001138 1014540
## 3 Intron (ENST00000379370.7/375790, intron 1 of 35)       1   1020123 1056118
## 4 Exon (ENST00000412397.2/100288175, exon 10 of 10)       1   1070967 1074306
## 5                                Downstream (1-2kb)       1   1082146 1084072
## 6                                          Promoter       1   1292390 1305929
##   geneLength geneStrand geneId      transcriptId distanceToTSS         ENSEMBL
## 1       6890          2  84808 ENST00000341290.6         -1383 ENSG00000187642
## 2      13403          1   9636 ENST00000624697.4             0 ENSG00000187608
## 3      35996          1 375790 ENST00000620552.4           791 ENSG00000188157
## 4       3340          2 401934 ENST00000453464.3          4114 ENSG00000237330
## 5       1927          2  54991 ENST00000464905.1          3985 ENSG00000131591
## 6      13540          2 116983 ENST00000492936.5             0 ENSG00000131584
##     SYMBOL                                                 GENENAME
## 1    PERM1             PPARGC1 and ESRR induced regulator, muscle 1
## 2    ISG15                            ISG15 ubiquitin like modifier
## 3     AGRN                                                    agrin
## 4   RNF223                                  ring finger protein 223
## 5 C1orf159                      chromosome 1 open reading frame 159
## 6    ACAP3 ArfGAP with coiled-coil, ankyrin repeat and PH domains 3

All peak information contained in the peak list will be retained in the output of annotatePeak. Positions and strand information of nearest genes are reported. The distance from peak to the TSS of its nearest gene is also reported. The genomic region of the peak is reported in annotation column. Since some annotation may overlap, ChIPseeker adopted the following priority in genomic annotation :

  • Promoter
  • 5’ UTR
  • 3’ UTR
  • Exon
  • Intron
  • Downstream
  • Intergenic
  • Downstream is defined as the downstream of gene end.

This hierachy can be customized using the parameter genomicAnnotationPriority.

The annotatePeak function report(s) detail(led) information when the annotation is Exon or Intron(.) (F)or tion when the annotation is Exon or Intron, for instance “Exon (uc002sbe.3/9736, exon 69 of 80),” means that the peak is overlap(ing) with an Exon of transcript uc002sbe.3, (whose) corresponding Entrez gene ID is 9736 (Transcripts that belong to the same gene ID may differ in splice events), and this (contains/harbors).

The “annoDb” parameter is optional. If it is provided, some extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added.

–

Reminder: The TxDb class is a container for storing transcript annotations.

  • Bioconductor provides several packages containing TxDb objects for model organisms sur as Human and mouse. For instance, TxDb.Hsapiens.UCSC.hg38.knownGene, TxDb.Hsapiens.UCSC.hg19.knownGene for human genome hg38 and hg19, TxDb.Mmusculus.UCSC.mm10.knownGene and TxDb.Mmusculus.UCSC.mm9.knownGene for mouse genome mm10 and mm9, etc.

  • User can also prepare their own TxDb by retrieving information from UCSC Genome Bioinformatics and BioMart data resources by R function makeTxDbFromBiomart and makeTxDbFromUCSC.

  • One can also create a TxDb objects for his favourite organism using an annotation file in GTF/GFF format using the function makeTxDbFromGFF or the package GenomicFeatures.

Expand to find Coturnix japonica example 
## 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')

–

2.6.2 Visualize genomic annotation

We can now perform a more detailed analysis of genomic features associated to our peaks (introns, exons, promoters, distal regions,…) visualizing the genomic distribution either as a pie chart with plotAnnoPie or as a bar chart plotAnnoBar.

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

# as a barplot
plotAnnoBar(peakAnno$BRG1)

However since some annotation overlap, ChIPseeker provides functions that help having a view of full annotation overlap.

library(UpSetR)
upsetplot(peakAnno$BRG1)

2.7 Heatmaps: visualisation of binding profiles

Heatmap(s) are widely used representation(s) of ChIP-seq data as they allow simultaneous visualization of read enrichment at various locations. For instance, one may want to represent reads related to a chipped protein in regions spanning +/-5Kb around all TSS of the reference genome. Another objective would be to compare read enrichment at the same locations in many chip-seq datasets.

2.7.1 BRG1 peak distribution

We want to know if BRG1 is binding large and/or narrow regions, unique and/or tandem etc.

First of all, for calculating the profile of ChIP peaks binding to the BRG1 center, we need to define peak centers and extend them each side. We will select the 10000 best q-value peaks (column V9) to lighten the matrix. Then (we need to) align the peaks that are mapped to these regions, and to generate the tagMatrix. The resulting matrix will be 2000*2+1 columns and 10000 lines which is quite big to display as a heatmap.

Step by step:
- Compute the center of BRG1 peaks using the GenomicRanges methods start and end. Assign the results to a new column in metadata of the BRG1 peaks object.

Show code: compute the peak centers 
# Compute the center of the peaks and attribute it to a new
# column in the metadata of the BRG1 peak GRanges
peaks$BRG1$center.peak <- (start(peaks$BRG1) + end(peaks$BRG1))/2
  • Sort the peaks in decreasing order (column V9) and select the top 10,000.
Show code: order and select peaks 
# For computation and memory efficiency reasons,
# we subset the top 10K peaks according to the FDR column (V9)
top.10000 <- peaks$BRG1[order(peaks$BRG1$V9,decreasing=T)][1:10000]
  • Create a GRanges object containing genomic positions of the 10K peak centers. You need to define seqnames and ranges in an IRanges object: IRanges(start,end).
Show code: Create GRanges object with peak center 
# Generate peak center GRanges for the 10K top peaks
centers.BRG1 <- GRanges(seqnames(top.10000),
                        IRanges(start=top.10000$center.peak,
                                end=top.10000$center.peak))
  • Extend the ranges each side of 2000 bp using the function resize. Fill the parameters width and fix.
Show code: resize the peak regions 
# Extend each side of 2000 bp
extended.2K.BRG1 <- resize(centers.BRG1, width = width(centers.BRG1)+4000, fix = "center")
  • You now have the regions you want to map the peaks to. Create the matrix using the getTagMatrix function and the parameter windows. Display it with tagHeatmap and its parameters xlim and color which correspond to the presence of peaks in the region.
Show code: create the matrix and display the heatmap 
## compute the density of peaks within the promoter regions
tagMatrix <- getTagMatrix(peaks$BRG1, windows=extended.2K.BRG1)
## plot the density
tagHeatmap(tagMatrix, xlim=c(-2000, 2000), color="red")

The regions are ordered relative to their peak enrichment.

We can display a summary of the binding profiles by looking at the corresponding average profiles. This kind of profiles is much less greedy, we can thus extend a bit more (e.g. +/- 5000) from the peak centers redoing the previous steps. Definition of the regions have to be redone. Try to create the GRanges and the tagMatrix.

Show code: Extend 5000bp from TSS and create the matrix 
# Extend from the peak center
extended.5K.BRG1 <- resize(centers.BRG1, width = width(centers.BRG1)+10000, fix = "center")
## compute the density of peaks within the promoter regions
tagMatrix <- getTagMatrix(peaks$BRG1, windows=extended.5K.BRG1)


Use the plotAvgProf function defining the xlim and setting the labels of the x and y axis.

Show code: Plot the average profile of BRG1 binding 
# Plot the profile
plotAvgProf(tagMatrix, xlim=c(-5000, 5000),
            xlab="Distance to peak center", ylab = "Peak Count Frequency")

It looks like BRG1 is having several binding patterns but the binary nature of the signal (presence/absence of peaks) and row ordering do not allow us to appreciate them.

2.7.2 Read enrichment in BRG1 peaks

For computation and memory efficiency reasons, we are not going to analyse read coverage at the nucleotide resolution. The strategy is rather to compute coverage in equally sized windows (e.g. 20nt). We thus need to build a matrix composed of rows that are all BRG1 peaks and columns that contain read enrichment in all bins.

This allows to consider a bigger set of peaks and covered region. We will now analyze the whole set of BRG1 peaks over 10Kb (+/- 5000bp). We will ask for a hundred bins each side of the center resulting in 200 windows of 50 bp.

2.7.2.1 Prepare the signal matrices

We need to load bigwig files for all datasets that we want to visualize. Data are imported using a function from the rtracklayer package.


First, create the GRanges object containing centers of all the peaks and extend them each side of 5000 bp.

Show code: Create all peak center GRanges object and extend them 
# # Generate peak center GRanges for all the peaks
centers.BRG1 <- GRanges(seqnames(peaks$BRG1),
                        IRanges(start=peaks$BRG1$center.peak,
                                end=peaks$BRG1$center.peak))
# Extend each side of 5000 bp
extended.5K.BRG1 <- resize(centers.BRG1, width = width(centers.BRG1)+10000, fix = "center")

We then load the rtracklayer library and import the read coverage from a bigWig file.

# load the library
library(rtracklayer)
# load the bw file for BRG1
cvg.BRG1 <- import(paste0(params$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 1333289 runs
##   Lengths: 979000      3     48     15     32 ...   1288     73     36  22594
##   Values :   0.00   0.32   0.16   0.49   0.32 ...   0.00   0.17   0.33   0.00
## 
## $chr10
## numeric-Rle of length 133797422 with 510080 runs
##   Lengths: 368603    154     85    154    304 ...     90     58    154 272220
##   Values :   0.00   0.16   0.00   0.16   0.00 ...   0.65   0.00   0.33   0.00
## 
## $chr11
## numeric-Rle of length 135086622 with 542307 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 564715 runs
##   Lengths: 194505     18    115    154    136 ...     22     68     19 189737
##   Values :   0.00   0.32   0.00   0.16   0.00 ...   0.16   0.65   0.49   0.00
## 
## $chr13
## numeric-Rle of length 114364328 with 338291 runs
##   Lengths: 19776001       19        5       12 ...       10       13    43427
##   Values :     0.00     1.29     1.13     0.97 ...     0.16     0.65     0.00
## 
## ...
## <20 more elements>

We can now create a matrix of read enrichment at the positions of interest with featureAlignedSignal from the ChIPpeakAnno package.


Load the library and find the required parameters. Be careful to the type of the required objects.

Show code: Create the read enrichment matrix 
# load library
library(ChIPpeakAnno)
# featureAlignedSignal needs the coverage (cvg) stored in a list.
cvglist <- list(BRG1=cvg.BRG1)
sig <- featureAlignedSignal(cvglists = cvglist, feature.gr = centers.BRG1,n.tile=200,
                           upstream=5000, downstream=5000)
dim(sig$BRG1)
## [1] 30874   200

2.7.2.2 Create heatmaps

Let’s draw the heatmaps using the EnrichedHeatmap library. We need first to transform our enrichment matrix in an object readable for the EnrichedHeatmap function. We use the function as.normalizedMatrix for which many parameters need to be set. What do they control ?

## Load the library
library(EnrichedHeatmap)
## Create a list of normalizedMatrix that is the input format
## for EnrichedHeatmap
mat1 <- list()
mat1[["BRG1"]] <- as.normalizedMatrix(as.matrix(sig[["BRG1"]]),
    k_upstream = 100,
    k_downstream = 100,
    k_target = 0,
    extend = c(5000, 5000),
    signal_name = names(sig[["BRG1"]]),
    target_name = "Peak centers"
)
## Create the Heatmap with default parameters
EnrichedHeatmap(mat1$BRG1, name = "BRG1")

EnrichedHeatmap combines the average profile and the density heatmap. We can observe a greater precision of the signal around the peak centers. As with ChIPpeakAnno, by default, heatmaps are sorted by read enrichment. However, it would be worth grouping together regions that have similar read enrichment patterns. This can be done using a clustering method such as k-means. This type of clustering requires the number of expected clusters to be set. Moreover, to obtain reproducible clustering results, we need to set a seed.

## define a seed value to get the same results when re-running the analysis
set.seed(123)
## Create Heatmaps with k-means clustering on BRG1 data
## We keep the generated object in order to use the clustering
## information.
heatmap.kmeans <- EnrichedHeatmap(mat1$BRG1, name = "BRG1", row_km = 8,
    column_title = "BRG1", row_title_rot = 0)
## draw the heatmap
htlist <- draw(heatmap.kmeans)
## Warning: did not converge in 10 iterations

Let’s enhance it!

  • We first retrieve the peak clusters from htlist object to use them as a pre-defined partition in our enhanced heatmap.
# Use the row_order function to retrieve peaks index
# belonging to each cluster
clusters <- row_order(htlist)
# rename the clusters
names(clusters) <- paste0("cluster",names(clusters))
# transform to a vector
# Check the class and length of the resulting object
class(clusters)
## [1] "list"
length(clusters)
## [1] 8
# Check what is actually in the list elements
head(clusters[[1]])
## [1]  1353  1815 22265  1564 19335 19294
# Each element of the list contain the indexes of the peaks in the
# original object
# Check cluster sizes
lapply(clusters,length)
## $cluster1
## [1] 1214
## 
## $cluster2
## [1] 1462
## 
## $cluster3
## [1] 1691
## 
## $cluster4
## [1] 3055
## 
## $cluster5
## [1] 2035
## 
## $cluster6
## [1] 5363
## 
## $cluster7
## [1] 15114
## 
## $cluster8
## [1] 940
# create the partition by transforming the list in vector
partition <- unlist(clusters)
head(partition)
## cluster11 cluster12 cluster13 cluster14 cluster15 cluster16 
##      1353      1815     22265      1564     19335     19294
# names of the elements were extended with a number
# We thus trim them to retrieve the cluster name
names(partition) <- substring(first = 1,last = 8,text=names(partition))
# The numbers are the indexes of the rows, we need to
# sort the indexes to get the right order of cluster
# labels
partition <- names(partition)[order(partition)]

# keep this information with you BRG1 peaks
peaks$BRG1$cluster <- partition
head(peaks$BRG1)
## GRanges object with 6 ranges and 9 metadata columns:
##       seqnames          ranges strand |                     V4        V5
##          <Rle>       <IRanges>  <Rle> |            <character> <integer>
##   [1]     chr1   983476-984526      * | BRG1siCTRL_CHIP-seq_..       108
##   [2]     chr1 1000729-1001179      * | BRG1siCTRL_CHIP-seq_..        90
##   [3]     chr1 1020914-1021207      * | BRG1siCTRL_CHIP-seq_..       114
##   [4]     chr1 1069156-1070192      * | BRG1siCTRL_CHIP-seq_..        96
##   [5]     chr1 1079517-1080087      * | BRG1siCTRL_CHIP-seq_..       179
##   [6]     chr1 1304853-1306013      * | BRG1siCTRL_CHIP-seq_..       108
##                V6        V7        V8        V9       V10 center.peak
##       <character> <numeric> <numeric> <numeric> <integer>   <numeric>
##   [1]           .   5.97792   13.5966  10.83201       384      984001
##   [2]           .   6.65677   11.6289   9.01201       199     1000954
##   [3]           .   8.11634   14.2980  11.46427       158     1021060
##   [4]           .   6.90466   12.3485   9.65604       410     1069674
##   [5]           .   9.95583   21.2517  17.96153       213     1079802
##   [6]           .   5.97792   13.5966  10.83201       176     1305433
##           cluster
##       <character>
##   [1]    cluster7
##   [2]    cluster7
##   [3]    cluster7
##   [4]    cluster7
##   [5]    cluster7
##   [6]    cluster6
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths
  • Now, customize the enrichment colors and differentiate clusters with colors.
library(circlize)
# Define new colors for each heatmap
col_brg1 <- colorRamp2(c(0,10,15), c("white", "blue","black"))
# create a legend for the cluster labels
lgd <- Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5","cluster6","cluster7","cluster8"),
    title = "Clusters",
    type = "lines", legend_gp = gpar(col = 2:9))
# Add a first column containing the cluster assignment
ht_list <- Heatmap(partition, col = structure(2:9, names = paste0("cluster", 1:8)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(mat1$BRG1, name = "BRG1", col=col_brg1, # specify the heat colors
  # color per cluster
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
  column_title = "BRG1")
draw(ht_list, split = partition, annotation_legend_list = list(lgd),
    ht_gap = unit(c(2, 8), "mm"))

2.7.2.3 Are BRG1 binding profiles associated with particular genomic regions ?

We will use here the clusters obtained with EnrichedHeatmap and look at the genomic distribution of peaks using TxDb.Hsapiens.UCSC.hg38.knownGene and plotAnnoBar.

# The following code uses each element of the list (peaks indexes)
# to select corresponding peaks in the original data and return them
BRG1.clusters <- lapply(clusters,function(x,peaks){
  # x represents one element of the list, one set of
  # indexes
  cl <- peaks[x]
  return(cl)
},peaks=peaks$BRG1)
# check names of the resulting list
names(BRG1.clusters)
## [1] "cluster1" "cluster2" "cluster3" "cluster4" "cluster5" "cluster6" "cluster7"
## [8] "cluster8"
# Transform the list as GRangesList to be able to use it with the annotatePeak
# function
BRG1.clusters <- as(BRG1.clusters,"GRangesList")
# apply annotatePeak function to each element of the BRG1.clusters GRangesList
peakAnnoList <- lapply(BRG1.clusters, annotatePeak, TxDb=txdb,
                       tssRegion=c(-1000, 100), annoDb="org.Hs.eg.db",
                       verbose=FALSE)
plotAnnoBar(peakAnnoList)

3 Data integration: let’s compare BRG1 to MITF and SOX10

3.1 Dataset description

3.1.1 Load MITF and SOX10 datasets

Using the readPeakFile function, add MITF (MITF_CHIP-seq_peaks.narrowPeak) and SOX10 (SOX10_CHIP-seq_peaks.narrowPeak) peaks to the peaks list. The path to these files is stored in the params list as the chipseq element.

Show code: Read MITF and SOX10 peak files 
## Let's load MITF and SOX10 peak sets
peaks[["MITF"]] <- readPeakFile(paste0(params$chipseq, params$mitf), as="GRanges")
peaks[["SOX10"]] <- readPeakFile(paste0(params$chipseq, params$sox10), as="GRanges")

The list peaks now contains 3 elements.

length(peaks)
## [1] 3
names(peaks)
## [1] "BRG1"  "MITF"  "SOX10"

3.1.2 How many peaks were called?

We want to answer this question using a bar plot with ggplot2. Here we are basically plotting the same graphs as in section 2.2 with more data sets.

Try to create the command lines enabling to get a data.frame containing the number of peaks per chipped protein (use the sapply function to compute the length of each peak set). Then, plot the results as a bar plot filled with colors from the brewer palette “Set1” according to the TF name.

Show code: Plot peak set size per TF 
# The `sapply` function makes now completely sense as we have a list with
# several elements to which we want to apply functions, here the length function.
peak.lengths <- data.frame(IP=names(peaks),
                           NbPeaks=sapply(peaks,length))
# check the object
peak.lengths
##          IP NbPeaks
## BRG1   BRG1   30874
## MITF   MITF    9702
## SOX10 SOX10    4538
ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
         geom_col()+
         scale_fill_brewer(palette="Set1")+
        theme_classic()

3.1.3 How large are these peaks?

Compute statistics on all peak sizes as in section 2.3.

  • Compute the peak width within each peak set
  • Look at the structure of the return object (str())
  • Use the melt function from the reshape2 package to obtain a data.frame with one column with all the peak width and one column containing the name of the peak set they come from.
Show code: create a data.frame containing all peak widths 
peak.width = lapply(peaks,width)
str(peak.width)
## List of 3
##  $ BRG1 : int [1:30874] 1051 451 294 1037 571 1161 360 285 897 905 ...
##  $ MITF : int [1:9702] 260 190 262 175 195 257 233 188 1251 317 ...
##  $ SOX10: int [1:4538] 252 272 331 517 253 362 246 306 254 400 ...
peak.width.table <- melt(peak.width)
head(peak.width.table)
##   value   L1
## 1  1051 BRG1
## 2   451 BRG1
## 3   294 BRG1
## 4  1037 BRG1
## 5   571 BRG1
## 6  1161 BRG1

In this table, there are as many rows as the total number of peaks in all peak sets. It contains all possible pairs IP <-> number of peaks possible.

# Number of peaks per chipped protein
sapply(peaks, length)
##  BRG1  MITF SOX10 
## 30874  9702  4538
# total number of peaks
sum(sapply(peaks, length))
## [1] 45114
# size of the table we've just generated
dim(peak.width.table)
## [1] 45114     2

We now want to represent the log transformed peak length distribution as boxplots filled with the corresponding colors. See bellow:

Show code: plot the distribution on peak length per experiment 
ggplot(peak.width.table, aes(x=L1, y=value, fill=L1)) +
         geom_boxplot()+
         theme_classic()+
         scale_fill_brewer(palette="Set1")+
         labs(x = "", y = "log10(Peak sizes)", fill = "")+
         scale_y_log10()

3.2 Peak filtering

In order to discard poor quality peaks, a threshold is set on qValue. The following code aims at selecting peaks with a qValue greater than 4. The qValue corresponds to the 9th column of narrowPeak files. So, we are going to set a threshold on this.

## We will filter the peaks in each set using lapply. We first define
## a function to be applied to each set

filter_V9 <- function(x){
    res <- x[x$V9 >= 4,]
    return(res)
}

## Select high quality peaks in each element of the list peaks
peaks <- lapply(peaks, filter_V9)
## Compute the number of remaining peaks
sapply(peaks, length)
##  BRG1  MITF SOX10 
## 30874  9279  4214

3.2.1 Peak annotation

Load and add MITF and SOX10 annotations to the peakAnno list that already contains BRG1 annotations using annotatePeak as in 2.6.1.

Show code: annotate the peaks 
peakAnno[["MITF"]] = annotatePeak(peaks$MITF, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")
## >> preparing features information...      2021-05-25 19:13:12 
## >> identifying nearest features...        2021-05-25 19:13:12 
## >> calculating distance from peak to TSS...   2021-05-25 19:13:12 
## >> assigning genomic annotation...        2021-05-25 19:13:12 
## >> adding gene annotation...          2021-05-25 19:13:17 
## >> assigning chromosome lengths           2021-05-25 19:13:18 
## >> done...                    2021-05-25 19:13:18
peakAnno[["SOX10"]] = annotatePeak(peaks$SOX10, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")
## >> preparing features information...      2021-05-25 19:13:18 
## >> identifying nearest features...        2021-05-25 19:13:18 
## >> calculating distance from peak to TSS...   2021-05-25 19:13:18 
## >> assigning genomic annotation...        2021-05-25 19:13:18 
## >> adding gene annotation...          2021-05-25 19:13:23 
## >> assigning chromosome lengths           2021-05-25 19:13:23 
## >> done...                    2021-05-25 19:13:23

And plot the distribution of genomic feature of the peaks as a barplot.

Show code: Plot genomic feature distribution per experiment 
plotAnnoBar(peakAnno)

3.3 Are BRG1, MITF and SOX10 co-localizing ?

3.3.1 Compare BRG1, MITF and SOX10 peak positions (Venn Diagram)

Overlaps between peak datasets can be evaluated using a Venn diagram. Such approach is implemented in the ChIPpeakAnno package.

library(ChIPpeakAnno)
# We first compute the overlap between peak sets, keeping the information
# of all peaks overlapping in each set (see ?findOverlapsOfPeaks for help)
ovl <- findOverlapsOfPeaks(peaks, connectedPeaks="keepAll")

ChIPpeakAnno imposes, while plotting the Venn diagram, to compute the significance of the pairwise associations using a hypergeometric test. To this end, we need to estimate the number of all potential binding events which is used by the makeVennDiagram function through the totalTest number. It is used for the hypergeometric sampling that is used to determine if the overlap between two datasets is more than would be expected by chance. This is not a trivial question, the answer is driven by what you know about the binding properties of your factors (eg. sequence specific, mainly intergenic etc). You can find an interesting discussion here. In our case we can refer to the genomic distribution of the peaks that we have plotted previously. We can assume that our TFs have a gene body binding preference. Genes cover roughly 10% of the genome.

# 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)

TIPS: We can define the colors attributed to each set using the function colours(). In any case you want to set a color you can use this function. Please, have a look to this page or more generally this page if you are interested in finely customizing colors of your plots.

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"))

## $p.value
##      BRG1 MITF SOX10          pval
## [1,]    0    1     1 7.824528e-215
## [2,]    1    0     1  0.000000e+00
## [3,]    1    1     0  0.000000e+00
## 
## $vennCounts
##      BRG1 MITF SOX10 Counts count.BRG1 count.MITF count.SOX10
## [1,]    0    0     0 411865          0          0           0
## [2,]    0    0     1   1127          0          0        1127
## [3,]    0    1     0   6930          0       6930           0
## [4,]    0    1     1     40          0         40          40
## [5,]    1    0     0  24781      24781          0           0
## [6,]    1    0     1   2547       3491          0        2572
## [7,]    1    1     0   1784       1907       1822           0
## [8,]    1    1     1    454        695        487         475
## attr(,"class")
## [1] "VennCounts"

According to the hypergeometric test p-values all pairwise comparisons are highly significant (=0). This has to be taken very carefully as it depends largely on the background estimation that we may have over-estimated. If you are not sure about your estimation, you will prefer to use a non-parametric approach based on your peak genomic distribution to estimate randomness. To this end the TxDb.Hsapiens.UCSC.hg38.knownGene is used.

3.3.1.1 Is the overlap significant?

ChIPpeakAnno provides the preparePool and peakPermTest functions to compute this test. These tests are made by pair, we will look at SOX10/MITF overlap significance as an example.

# 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=1, force.parallel=FALSE)
SOX10.MITF
## $cntOverlaps
## P-value: 0.0099009900990099
## Z-score: 139.5137
## Number of iterations: 100
## Alternative: greater
## Evaluation of the original region set: 437
## Evaluation function: cntOverlaps
## Randomization function: randPeaks
## 
## attr(,"class")
## [1] "permTestResultsList"
plot(SOX10.MITF)

Venn diagrams are widely used to represent overlaps, intersections. However, in ChIP-seq analysis, the definition of the peaks is dependent on the peak caller, the FDR thresholds (what about peaks just bellow the threshold). The overlap is also difficult to assess, indeed do we call an overlap a 1 nucleotide an intersection ? It is one of the numerous parameters that can be tuned …

Are different combinations of TFs bind specific genomic regions ?

coocs <- as(ovl$peaklist,"GRangesList")
# apply annotatePeak to each set of peaks in the list coocs
peakAnnoList <- lapply(coocs, annotatePeak, TxDb=txdb,
                        tssRegion=c(-1000, 100), annoDb="org.Hs.eg.db",
                        verbose=FALSE)
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
# Plot peak distribution relatively to gene features
plotAnnoBar(peakAnnoList)

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

3.3.2 Heatmaps / Profiles

One way to circumvent hard thresholds and relatively arbitrary choices, we can choose to use heatmap and average profile representations. To this end, we need to define the reference point of view. We use the BRG1 peak centers.

Load bigwig file for MITF and SOX10 at BRG1 peaks.

# Load and add bigwig profiles for MITF and SOX10 in cvglist already containing
# BRG1
names(cvglist)
## [1] "BRG1"
cvglist$MITF <- import(file.path(params$chipseq,params$mitf.bw),format="BigWig",
                      which=extended.5K.BRG1,
                      as="RleList")
cvglist$SOX10 <- import(file.path(params$chipseq,params$sox10.bw),format="BigWig",
                      which=extended.5K.BRG1,
                      as="RleList")
names(cvglist)
## [1] "BRG1"  "MITF"  "SOX10"

Prepare the matrices binned in 50bp windows.

sig <- featureAlignedSignal(cvglist, centers.BRG1,
                           upstream=5000, downstream=5000,n.tile=200)
lapply(sig, dim)
## $BRG1
## [1] 30874   200
## 
## $MITF
## [1] 30874   200
## 
## $SOX10
## [1] 30874   200

Create normalized matrices.

## Create a list of normalizedMatrix that is the input format
## for EnrichedHeatmap
BRG1.mat <- lapply(sig, function(x){
  # x represent each element of the sig list
  mat <- as.normalizedMatrix(as.matrix(x),
                      k_upstream = 100,
                      k_downstream = 100,
                      k_target = 0,
                      extend = c(5000, 5000),
                      #signal_name = names(sig[["MITF"]]),
                      target_name = "Peak center"
                      )
  return(mat)
})

We use the partition we computed in section 2.7.2.2

# Define new colors for each heatmap
col_sox10 <-colorRamp2(c(0,3,4), c("white", "blue","black"))
col_mitf <- colorRamp2(c(0,4,5), c("white", "blue","black"))
col_brg1 <- colorRamp2(c(0,5,6), c("white", "blue","black"))
# create a legend for the cluster labels
lgd <- Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5","cluster6","cluster7","cluster8"),
    title = "Clusters",
    type = "lines", legend_gp = gpar(col = 2:9))
# Add a first column containing the cluster assignment
ht_list <- Heatmap(partition, col = structure(2:9, names = paste0("cluster", 1:8)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(BRG1.mat$BRG1, name = "BRG1", col=col_brg1,
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
  column_title = "BRG1") +
EnrichedHeatmap(BRG1.mat$MITF, name = "MITF",
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
  column_title = "MITF", row_title_rot = 0, col=col_mitf) +
EnrichedHeatmap(BRG1.mat$SOX10, name = "SOX10",
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
  column_title = "SOX10", row_title_rot = 0, col=col_sox10)
draw(ht_list, split = partition, annotation_legend_list = list(lgd),
    ht_gap = unit(c(2, 8, 8,8), "mm"))

Lets have a look at the TF binding profiles in cluster 5.

# For each matrix, filter peaks from cluster 5
sig.tmp <- lapply(sig,function(x,part){
    mat <- x[part=="cluster5",]
    return(mat)
  },part=partition)

# plot the average profiles
featureAlignedDistribution(sig.tmp, centers.BRG1,
                           upstream=5000, downstream=5000,
                           type="l",n.tile=200)

3.4 Is gene expression influenced by TF binding to their promoters ?

  • We will concentrate on genes having BRG1 binding close to their TSS. To do that we will plot a heatmap centered on TSS.

First retrieve the TSS positions from the TxDB object

# Generate TSS and promoter GRanges, the function promoters allows to retrieve the
# gene ID (entrez)
TSS <- GenomicFeatures::promoters(txdb, upstream=0, downstream=0, columns="gene_id")

Clean the TSS list removing TSS with no gene ID assigned (column gene_id in the metadata) and on scaffolds (seqnames).
The column gene_id is of type characterList, this type is difficult to manipulate. We then decide to transform it to a character vector. We can’t directly unlist the object because, to save memory, not annotated TSS are assigned empty elements which disappear when unlisting.

# For each element of the characterList we check it's length, if it is 0
# we fill the element with NA. We the unlist.
gene_ids <- unlist(lapply(TSS$gene_id,function(x){
  if (length(x)==0){
     NA
  }
  else {
    x
  }
} ))
# Keep TSS on selected chromosomes annotated with a gene ID and remove
# duplicate IDs due to transcript isoforms. Here, the selection of the
# isoform is random.
TSS <- TSS[as.vector(seqnames(TSS))%in%paste0("chr",c(1:22,"X","Y")) & !is.na(gene_ids) & !duplicated(gene_ids)]

# extend TSS both direction
TSS.extended <- resize(TSS, width = width(TSS)+4000, fix = "center")

Load bigwig profiles for TSS extended regions, create matrices and transform in the enrichedHeatmap recognized format.

# load the bw file for all TFs within the promoter
cvglist <- list()
cvglist$BRG1 <- import(file.path(params$chipseq,params$brg1.bw),format="BigWig",
                      which=TSS.extended,
                      as="RleList")
cvglist$MITF <- import(file.path(params$chipseq,params$mitf.bw),format="BigWig",
                      which=TSS.extended,
                      as="RleList")
cvglist$SOX10 <- import(file.path(params$chipseq,params$sox10.bw),format="BigWig",
                      which=TSS.extended,
                      as="RleList")
# Produce the signal matrices
sig <- featureAlignedSignal(cvglist, TSS,n.tile=400,
                           upstream=2000, downstream=2000)
# Transform the signal matrices as normalizedMatrix
TSS.mat <- lapply(sig, function(x){
  # x represent each element of the sig list
  mat <- as.normalizedMatrix(as.matrix(x),
                      k_upstream = 200,
                      k_downstream = 200,
                      k_target = 0,
                      extend = c(2000, 2000),
                      target_name = "TSS"
                      )
  return(mat)
})

We will compute K-means clustering with kmeans function to define the partitions of the TSS based on BRG1 signal.

# Compute a partition using the kmeans function, asking for 5 clusters
set.seed(20210526)
partition.TSS = paste0("cluster", kmeans(TSS.mat$BRG1, centers = 5)$cluster)
# Specify colors for the position enrichment for each matrix
col_sox10 = colorRamp2(c(0,3,4), c("white", "blue","black"))
col_mitf = colorRamp2(c(0,4,5), c("white", "blue","black"))
col_brg1 = colorRamp2(c(0,5,6), c("white", "blue","black"))
lgd = Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5"),
    title = "Clusters",
    type = "lines", legend_gp = gpar(col = 2:6))
ht_list = Heatmap(partition.TSS, col = structure(2:7, names = paste0("cluster", 1:5)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(TSS.mat$BRG1, name = "BRG1", col=col_brg1,
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
  column_title = "BRG1") +
EnrichedHeatmap(TSS.mat$MITF, name = "MITF",
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
  column_title = "MITF", row_title_rot = 0, col=col_mitf) +
EnrichedHeatmap(TSS.mat$SOX10, name = "SOX10",
  top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
  column_title = "SOX10", row_title_rot = 0, col=col_sox10)
draw(ht_list, split = partition.TSS, annotation_legend_list = list(lgd),
    ht_gap = unit(c(2, 8, 8,8), "mm"))

Lets have a look at the TF binding profiles in cluster 4.

# For each matrix, filter peaks from cluster 5
sig.tmp <- lapply(sig,function(x,part){
    mat <- x[part=="cluster4",]
    return(mat)
  },part=partition.TSS)

# plot the average profiles
featureAlignedDistribution(sig.tmp, TSS,
                           upstream=2000, downstream=2000,
                           type="l",n.tile=400)
## Warning in featureAlignedDistribution(sig.tmp, TSS, upstream = 2000, downstream
## = 2000, : feature.gr is set to the center of feature.gr

  • Let’s integrate expression data from RNA-seq experiment.

The goal is to visualize the gene expression according to the cluster the TSS belong.

# load RNA-seq normalized expression data
load(params$rnaseq) # RNAseq
head(RNAseq)
##   Ensembl.Gene.ID    WT.norm shBGR1.norm   logCPM      logFC          FDR
## 1 ENSG00000000003  722.38046   119.10483 4.387012  2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280  2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457  259.21973   314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460  333.22470   380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971   13.42921    13.59959 0.341927  0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492  1521.75139 5.567526 -0.3587066 3.117933e-06
##   signif
## 1 stable
## 2 stable
## 3 stable
## 4 stable
## 5 stable
## 6 stable

RNA-seq columns: - Ensembl.Gene.ID: Gene idenfier from ensEMBL (ENSG*****) - WT.norm: normalized expression for wild type cell lines - shBGR1.norm: normalized expression for shBRG1 cell lines - logCPM: average log transformed overall expression - logFC: log fol change between the two conditions - FDR: multi-testing corrected p-value from edgeR statistical test - signif (up, down, stable): significance was set as: FDR < 0.01, logCPM > 5, abs(logFC) > 1

Reads from RNA-seq were assigned to genes modeled from ensEMBL and thus are identified with ensEMBL ID. However, the TSS are annotated with ENTREZ IDs, we thus need to map together the identifiers. We use an OrgDB object for the human genome the available information are displayed using columns().

columns(org.Hs.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GO"           "GOALL"        "IPI"          "MAP"          "OMIM"        
## [16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"       "UNIGENE"     
## [26] "UNIPROT"

We then use the function mapIds from AnnotationDbi package with the keys ENSEMBL and ENTREZID. You can see that other identifiers are available from this object (e.g: SYMBOL, REFSEQ, UNIPROT etc).

conv <- mapIds(x=org.Hs.eg.db,
              keys=RNAseq$Ensembl.Gene.ID, # what do we want to be mapped
              column="ENTREZID", # which type of ID we want
              keytype="ENSEMBL") # what type of ID we give
## 'select()' returned 1:many mapping between keys and columns
head(conv)
## ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000971 
##          "7105"          "8813"         "57147"         "55732"          "3075" 
## ENSG00000001036 
##          "2519"
# names of the vector are ENSEMBL IDs and elements are ENTREZ IDs

We use the match function that allows to return match between the IDs in RNAseq and the ensEMBL IDs from the conversion vector

m <- match(RNAseq$Ensembl.Gene.ID,names(conv))
head(m)
## [1] 1 2 3 4 5 6
# m contains NA for element in the first vector not present in the 2nd
# or the index of the element of the 2nd vector corresponding the those
# in the first
# Attribute the matched ENTREZ ID
RNAseq$ENTREZID <- conv[m]
head(RNAseq)
##   Ensembl.Gene.ID    WT.norm shBGR1.norm   logCPM      logFC          FDR
## 1 ENSG00000000003  722.38046   119.10483 4.387012  2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280  2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457  259.21973   314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460  333.22470   380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971   13.42921    13.59959 0.341927  0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492  1521.75139 5.567526 -0.3587066 3.117933e-06
##   signif ENTREZID
## 1 stable     7105
## 2 stable     8813
## 3 stable    57147
## 4 stable    55732
## 5 stable     3075
## 6 stable     2519

We can now associate RNA-seq and TSS information using the ENTREZ IDs and the match function.

# We add the partition (clusters) to the TSS object in order to add the
# clusters to the RNA-seq object
TSS$cluster <- partition.TSS
# The common IDs this time are the ENTREZ IDs
m <- match(RNAseq$ENTREZID,unlist(TSS$gene_id))
RNAseq$cluster <- TSS$cluster[m]
head(RNAseq)
##   Ensembl.Gene.ID    WT.norm shBGR1.norm   logCPM      logFC          FDR
## 1 ENSG00000000003  722.38046   119.10483 4.387012  2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280  2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457  259.21973   314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460  333.22470   380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971   13.42921    13.59959 0.341927  0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492  1521.75139 5.567526 -0.3587066 3.117933e-06
##   signif ENTREZID  cluster
## 1 stable     7105 cluster3
## 2 stable     8813 cluster4
## 3 stable    57147 cluster4
## 4 stable    55732 cluster3
## 5 stable     3075 cluster3
## 6 stable     2519 cluster4

We first want to compare the expression between this to condition in the different clusters.
We need to melt the RNAseq object keeping only the expression values in the “WT.norm” and the “shBGR1.norm” columns and the cluster assignment. We want to evaluate the differences using boxplots as bellow.

##    cluster variable      value
## 1 cluster3  WT.norm  722.38046
## 2 cluster4  WT.norm 1481.30280
## 3 cluster4  WT.norm  259.21973
## 4 cluster3  WT.norm  333.22470
## 5 cluster3  WT.norm   13.42921
## 6 cluster4  WT.norm 1157.52492

Show code: plot expression distribution in both condition 
# melt the RNAseq table to get a data.frame with a column indicating
# the experimental condition, a column with the cluster assignment
# and the expression values

expr.melt <- melt(RNAseq[,c("WT.norm","shBGR1.norm","cluster")])
head(expr.melt)

ggplot(expr.melt,aes(x=cluster,y=log2(value+1))) +
  geom_boxplot(aes(color=variable)) +
  theme_bw() + theme(legend.position = "top") +
  scale_color_manual(values=c(WT.norm="black",shBGR1.norm=colours()[613]))

We can see that there are no visible difference between the conditions even if cluster 3, which has the less BRG1 signal, is much less expressed than the others. We can check if the distribution of significantly differentially expressed genes differs among the clusters using barplots.

Show code: plot significantly differentially expressed genes in clusters 
# Distribution of differentially expressed genes as barplots
ggplot(RNAseq,aes(x=cluster)) +
  geom_bar(aes(fill=signif),position = "fill",color="black") +
  theme_bw() +
  scale_fill_manual(values=c(up="red",down="blue",stable="white"))

3.5 Are there particular biological function associated with the differentially expressed genes having different BRG1 binding patterns ?

To answer this question, we use the R package clusterProfiler to compute hypergeometric enrichment of biological function in each cluster using KEGG annotation database.
clusterProfiler provide the function compareCluster that allows to analyze and compare enrichment in different group of genes. This function recognize ENTREZ identifiers.


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

# load the library
library(clusterProfiler)
# Select genes whose gene expression is significantly changing
# between conditions
geneList <- RNAseq[RNAseq$signif!="stable",]
# The function split help at splitting a vector to a list following categories
# in an other vector.
geneList <- split(geneList$ENTREZID,geneList$cluster)
names(geneList)
## [1] "cluster1" "cluster2" "cluster3" "cluster4" "cluster5"
head(geneList$cluster1)
## [1] "27072" "54443" "57514" "57609" "998"   "58476"
compKEGG <- compareCluster(geneCluster   = geneList,
                         fun           = "enrichKEGG",
                         pvalueCutoff  = 0.05,
                         pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")

Let’s focus on cluster 4 and visualize which genes of the the KEGG pathway “Cell cycle” (hsa04110) are both differentially expressed and bound by BRG1 and MITF. We use the pathview package which maps differential gene expression values to KEGG maps and create a PNG file in your current directory that you get with getwd().

# load the library
library(pathview)
# Retrieve gene id from cluster 4
geneList.cl4 <- geneList$cluster4
# use match function to select correspondant logFC from RNA-seq data
m <- match(geneList.cl4,RNAseq$ENTREZID)
logFC.cl4 <- RNAseq$logFC[m]
names(logFC.cl4) <- geneList.cl4
# call the function
pathview(logFC.cl4,pathway.id = "hsa04110",species = "hsa")

Pathview

3.5.1 Are there published data I could use to enrich my analysis ?

  • Data Mining with ChIP seq data deposited in GEO with ChIPSeeker

There are many ChIP seq data sets that have been published and deposited in GEO database. 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.

We collect about 17,000 bed files deposited in GEO, user can use getGEOspecies to get a summary based on spieces.

getGEOspecies()
##                                               species  Freq
## 1                                       Aedes aegypti    11
## 2                                            Anabaena     6
## 3                                 Anolis carolinensis     5
## 4                                   Anopheles gambiae     2
## 5                                      Apis mellifera     5
## 6                           Apis mellifera scutellata     1
## 7                                  Arabidopsis lyrata     4
## 8                                Arabidopsis thaliana   288
## 9                                Atelerix albiventris     2
## 10                                         Bos taurus    37
## 11                          Branchiostoma lanceolatum    62
## 12                                      Brassica rapa    12
## 13                             Caenorhabditis elegans   189
## 14                                   Candida albicans    25
## 15                               Candida dubliniensis    20
## 16                             Canis lupus familiaris     7
## 17                          Chlamydomonas reinhardtii    51
## 18                               Chlorocebus aethiops     2
## 19                                 Cleome hassleriana     1
## 20                                      Columba livia     6
## 21                                  Crassostrea gigas     1
## 22                            Cryptococcus neoformans    51
## 23                                    Cyprinus carpio    40
## 24                                        Danio rerio   308
## 25                                 Drosophila busckii     1
## 26                            Drosophila melanogaster  1069
## 27            Drosophila melanogaster;\tSindbis virus     3
## 28                                 Drosophila miranda     2
## 29                           Drosophila pseudoobscura     7
## 30                                Drosophila simulans    12
## 31                                 Drosophila virilis    26
## 32                              Drosophila willistoni     1
## 33                                  Drosophila yakuba     8
## 34                                     Equus caballus     1
## 35                                   Escherichia coli    15
## 36                           Escherichia coli BW25113     4
## 37                              Escherichia coli K-12     2
## 38          Escherichia coli str. K-12 substr. MG1655    16
## 39                                      Gallus gallus    58
## 40                       Geobacter sulfurreducens PCA     3
## 41                                    Gorilla gorilla     2
## 42                                  Histophilus somni     1
## 43                                       Homo sapiens 29978
## 44                 Homo sapiens;\tHuman herpesvirus 8     6
## 45                               Human herpesvirus 6B     2
## 46                                Human herpesvirus 8     6
## 47                                Larimichthys crocea     7
## 48                             Legionella pneumophila     5
## 49                             Leishmania amazonensis     4
## 50                                   Leishmania major     2
## 51                   Leishmania major strain Friedlin     4
## 52                              Leishmania tarentolae    15
## 53                                     Macaca mulatta   120
## 54                              Monodelphis domestica     8
## 55                         Moraxella catarrhalis O35E     6
## 56                                                Mus     2
## 57                                       Mus musculus 16748
## 58                         Mus musculus x Mus spretus     1
## 59                         Mycobacterium tuberculosis     2
## 60                                    Myotis brandtii    15
## 61                              Naumovozyma castellii     1
## 62                             Nematostella vectensis    23
## 63                              Oreochromis niloticus     1
## 64                           Ornithorhynchus anatinus     5
## 65                                       Oryza sativa    30
## 66                                    Oryzias latipes     2
## 67                                    Pan troglodytes    93
## 68                                       Papio anubis     1
## 69                              Plasmodium falciparum   129
## 70                          Plasmodium falciparum 3D7    29
## 71                          Pseudomonas putida KT2440     2
## 72                                 Pseudozyma aphidis    11
## 73                                Pyrococcus furiosus     4
## 74                                  Rattus norvegicus   108
## 75                         Rhodopseudomonas palustris     6
## 76                  Rhodopseudomonas palustris CGA009     3
## 77                           Saccharomyces cerevisiae   813
## 78 Saccharomyces cerevisiae x Saccharomyces paradoxus    16
## 79            Saccharomyces cerevisiae;\tMus musculus    12
## 80                         Saccharomyces kudriavzevii     1
## 81                            Saccharomyces paradoxus     8
## 82                               Saccharomyces uvarum     1
## 83                      Schizosaccharomyces japonicus     2
## 84                          Schizosaccharomyces pombe   179
## 85                             Schmidtea mediterranea     7
## 86                               Solanum lycopersicum     2
## 87                                    Sorghum bicolor     2
## 88                              Spodoptera frugiperda    16
## 89                      Streptomyces coelicolor A3(2)     6
## 90                                         Sus scrofa    41
## 91                                Taeniopygia guttata     1
## 92                                   Tupaia chinensis     7
## 93                                    Vibrio cholerae     8
## 94                      Xenopus (Silurana) tropicalis    62
## 95                                     Xenopus laevis    10
## 96                                 Xenopus tropicalis    74
## 97                                           Zea mays    65

User can access the detail information by getGEOInfo, 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 significant of overlap. Parameter targetPeak can be the folder, e.g. hg19, 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. Signifcant overlap can be used to generate hypothesis of cooperative regulation. By mining the data deposited in GEO, we can identify some putative complex or interacted regulators in gene expression regulation or chromosome remodelling for further validation.

  • The (mod)ENCODE project link

  • The UCSC genome browser link

4 Importance of the annotation databases

In previous annotation steps, we have annotated peaks using annotation from UCSC’s knownGene. The database chosen for annotation can have an impact on subsequent results including data integration such as comparison between genes associated to peaks and gene expression using RNA-seq data for instance.

Let’s look at the overlap between genes expressed in RNA-seq data with genes associated to peaks (let’s remind that peaks are associated to closest genes if no other evidences are used).

## Download GTF files
library(GenomicFeatures)

## Import BRG1 peaks annotated with Ensembl and Refseq annotation
load(params$annot)

## Here is how annotation data were downloaded
# txdb.ensembl = makeTxDbFromGFF("data/ChIPseq/Homo_sapiens.GRCh38.103_UCSC_chr.gtf.gz")
# txdb.refseq = makeTxDbFromUCSC(genome="hg38", tablename="refGene")

## Here is how annotation was performed
# annot <- list()
# annot[["ensembl"]] <- annotatePeak(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb.ensembl, annoDb="org.Hs.eg.db")
# annot[["refseq"]] <- annotatePeak(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb.refseq, annoDb="org.Hs.eg.db")

## Add Gene symbols to RNAseq data
conv.symbol <- mapIds(x=org.Hs.eg.db,
              keys=RNAseq$Ensembl.Gene.ID, # what do we want to be mapped
              column=c("SYMBOL"), # which type of ID we want
              keytype="ENSEMBL") # what type of ID we give

m <- match(RNAseq$Ensembl.Gene.ID,names(conv.symbol))

RNAseq$SYMBOL <- conv.symbol[m]

## Now let's create a venn diagram that compares peak annotation with
## refseq or Ensembl
library(ggVennDiagram)
x <- list(RNAseq=RNAseq[RNAseq$signif=="up",]$SYMBOL,
          Ensembl=as.data.frame(annot$ensembl)$SYMBOL,
          Refseq=as.data.frame(annot$refseq)$SYMBOL
)

ggVennDiagram(x) + scale_fill_gradient(low="white",high = "blue")

5 Session info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] ggVennDiagram_1.1.1                     
##  [2] pathview_1.30.1                         
##  [3] clusterProfiler_3.18.1                  
##  [4] UpSetR_1.4.0                            
##  [5] reshape2_1.4.4                          
##  [6] RColorBrewer_1.1-2                      
##  [7] knitr_1.33                              
##  [8] org.Hs.eg.db_3.12.0                     
##  [9] TxDb.Hsapiens.UCSC.hg38.knownGene_3.10.0
## [10] GenomicFeatures_1.42.3                  
## [11] AnnotationDbi_1.52.0                    
## [12] Biobase_2.50.0                          
## [13] rtracklayer_1.50.0                      
## [14] circlize_0.4.12                         
## [15] ChIPseeker_1.26.2                       
## [16] ChIPpeakAnno_3.24.2                     
## [17] ggplot2_3.3.3                           
## [18] EnrichedHeatmap_1.20.0                  
## [19] GenomicRanges_1.42.0                    
## [20] GenomeInfoDb_1.26.7                     
## [21] IRanges_2.24.1                          
## [22] S4Vectors_0.28.1                        
## [23] BiocGenerics_0.36.1                     
## [24] ComplexHeatmap_2.6.2                    
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.1                             
##   [2] tidyselect_1.1.1                       
##   [3] RSQLite_2.2.7                          
##   [4] BiocParallel_1.24.1                    
##   [5] scatterpie_0.1.6                       
##   [6] munsell_0.5.0                          
##   [7] units_0.7-1                            
##   [8] withr_2.4.2                            
##   [9] colorspace_2.0-1                       
##  [10] GOSemSim_2.16.1                        
##  [11] highr_0.9                              
##  [12] DOSE_3.16.0                            
##  [13] MatrixGenerics_1.2.1                   
##  [14] labeling_0.4.2                         
##  [15] KEGGgraph_1.50.0                       
##  [16] GenomeInfoDbData_1.2.4                 
##  [17] polyclip_1.10-0                        
##  [18] bit64_4.0.5                            
##  [19] farver_2.1.0                           
##  [20] downloader_0.4                         
##  [21] vctrs_0.3.8                            
##  [22] generics_0.1.0                         
##  [23] lambda.r_1.2.4                         
##  [24] xfun_0.22                              
##  [25] BiocFileCache_1.14.0                   
##  [26] regioneR_1.22.0                        
##  [27] R6_2.5.0                               
##  [28] clue_0.3-59                            
##  [29] graphlayouts_0.7.1                     
##  [30] RVenn_1.1.0                            
##  [31] locfit_1.5-9.4                         
##  [32] AnnotationFilter_1.14.0                
##  [33] bitops_1.0-7                           
##  [34] cachem_1.0.4                           
##  [35] fgsea_1.16.0                           
##  [36] DelayedArray_0.16.3                    
##  [37] assertthat_0.2.1                       
##  [38] scales_1.1.1                           
##  [39] ggraph_2.0.5                           
##  [40] enrichplot_1.10.2                      
##  [41] gtable_0.3.0                           
##  [42] Cairo_1.5-12.2                         
##  [43] tidygraph_1.2.0                        
##  [44] ensembldb_2.14.1                       
##  [45] rlang_0.4.11                           
##  [46] GlobalOptions_0.1.2                    
##  [47] splines_4.0.5                          
##  [48] lazyeval_0.2.2                         
##  [49] BiocManager_1.30.14                    
##  [50] yaml_2.2.1                             
##  [51] qvalue_2.22.0                          
##  [52] RBGL_1.66.0                            
##  [53] tools_4.0.5                            
##  [54] ellipsis_0.3.2                         
##  [55] gplots_3.1.1                           
##  [56] jquerylib_0.1.4                        
##  [57] proxy_0.4-25                           
##  [58] Rcpp_1.0.6                             
##  [59] plyr_1.8.6                             
##  [60] progress_1.2.2                         
##  [61] zlibbioc_1.36.0                        
##  [62] classInt_0.4-3                         
##  [63] purrr_0.3.4                            
##  [64] RCurl_1.98-1.3                         
##  [65] prettyunits_1.1.1                      
##  [66] openssl_1.4.4                          
##  [67] GetoptLong_1.0.5                       
##  [68] viridis_0.6.0                          
##  [69] cowplot_1.1.1                          
##  [70] SummarizedExperiment_1.20.0            
##  [71] ggrepel_0.9.1                          
##  [72] cluster_2.1.2                          
##  [73] magrittr_2.0.1                         
##  [74] data.table_1.14.0                      
##  [75] futile.options_1.0.1                   
##  [76] DO.db_2.9                              
##  [77] ProtGenerics_1.22.0                    
##  [78] matrixStats_0.58.0                     
##  [79] hms_1.0.0                              
##  [80] evaluate_0.14                          
##  [81] XML_3.99-0.6                           
##  [82] VennDiagram_1.6.20                     
##  [83] gridExtra_2.3                          
##  [84] shape_1.4.5                            
##  [85] ggupset_0.3.0                          
##  [86] compiler_4.0.5                         
##  [87] biomaRt_2.46.3                         
##  [88] tibble_3.1.1                           
##  [89] KernSmooth_2.23-20                     
##  [90] crayon_1.4.1                           
##  [91] shadowtext_0.0.8                       
##  [92] htmltools_0.5.1.1                      
##  [93] tidyr_1.1.3                            
##  [94] DBI_1.1.1                              
##  [95] tweenr_1.0.2                           
##  [96] formatR_1.9                            
##  [97] dbplyr_2.1.1                           
##  [98] MASS_7.3-54                            
##  [99] rappdirs_0.3.3                         
## [100] sf_0.9-8                               
## [101] boot_1.3-28                            
## [102] Matrix_1.3-3                           
## [103] igraph_1.2.6                           
## [104] pkgconfig_2.0.3                        
## [105] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [106] rvcheck_0.1.8                          
## [107] GenomicAlignments_1.26.0               
## [108] xml2_1.3.2                             
## [109] bslib_0.2.4                            
## [110] multtest_2.46.0                        
## [111] XVector_0.30.0                         
## [112] stringr_1.4.0                          
## [113] digest_0.6.27                          
## [114] graph_1.68.0                           
## [115] Biostrings_2.58.0                      
## [116] rmarkdown_2.8                          
## [117] fastmatch_1.1-0                        
## [118] curl_4.3.1                             
## [119] Rsamtools_2.6.0                        
## [120] gtools_3.8.2                           
## [121] rjson_0.2.20                           
## [122] lifecycle_1.0.0                        
## [123] jsonlite_1.7.2                         
## [124] futile.logger_1.4.3                    
## [125] viridisLite_0.4.0                      
## [126] askpass_1.1                            
## [127] BSgenome_1.58.0                        
## [128] fansi_0.4.2                            
## [129] pillar_1.6.0                           
## [130] lattice_0.20-44                        
## [131] KEGGREST_1.30.1                        
## [132] fastmap_1.1.0                          
## [133] httr_1.4.2                             
## [134] plotrix_3.8-1                          
## [135] survival_3.2-11                        
## [136] GO.db_3.12.1                           
## [137] glue_1.4.2                             
## [138] png_0.1-7                              
## [139] Rgraphviz_2.34.0                       
## [140] bit_4.0.4                              
## [141] class_7.3-19                           
## [142] ggforce_0.3.3                          
## [143] stringi_1.5.3                          
## [144] sass_0.3.1                             
## [145] blob_1.2.1                             
## [146] caTools_1.18.2                         
## [147] memoise_2.0.0                          
## [148] dplyr_1.0.6                            
## [149] e1071_1.7-6

References

“(2014) Mod/Mouse/humanENCODE: Blacklisted Genomic Regions for Functional Genomics Analysis - Anshul Kundaje.” n.d. Accessed August 26, 2016. https://sites.google.com/site/anshulkundaje/projects/blacklists.
Anders, Simon, Paul Theodor Pyl, and Wolfgang Huber. 2015. “HTSeq—a Python Framework to Work with High-Throughput Sequencing Data.” Bioinformatics 31 (2): 166–69. https://doi.org/10.1093/bioinformatics/btu638.
Anders, and Huber. 2010. “Differential Expression Analysis for Sequence Count Data.” Genome Biology 11.
Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.
Heinz, Sven, Christopher Benner, Nathanael Spann, Eric Bertolino, Yin C. Lin, Peter Laslo, Jason X. Cheng, Cornelis Murre, Harinder Singh, and Christopher K. Glass. 2010. “Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and b Cell Identities.” Molecular Cell 38 (4): 576–89. https://doi.org/10.1016/j.molcel.2010.05.004.
Langmead, Ben, Cole Trapnell, Mihai Pop, and Steven L. Salzberg. 2009. “Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the Human Genome.” Genome Biology 10 (3): R25. https://doi.org/10.1186/gb-2009-10-3-r25.
Laurette, Patrick, Thomas Strub, Dana Koludrovic, Céline Keime, Stéphanie Le Gras, Hannah Seberg, Eric Van Otterloo, et al. 2015. “Transcription Factor MITF and Remodeller Brg1 Define Chromatin Organisation at Regulatory Elements in Melanoma Cells.” Edited by Michael R Green. eLife 4 (March): e06857. https://doi.org/10.7554/eLife.06857.
Martin, Marcel. 2011. “Cutadapt Removes Adapter Sequences from High-Throughput Sequencing Reads.” EMBnet.journal 17 (1): pp. 10–12. http://journal.embnet.org/index.php/embnetjournal/article/view/200.
Quinlan, Aaron R., and Ira M. Hall. 2010. “BEDTools: A Flexible Suite of Utilities for Comparing Genomic Features.” Bioinformatics 26 (6): 841–42. https://doi.org/10.1093/bioinformatics/btq033.
Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.