Convert peaks to GRanges

The ChIPpeakAnno package provides an example peak file in narrowPeak format, which is generated by MACS. To work with this example file, we need to locate it with system.file and convert it into a GRanges object using toGRanges function.

library(ChIPpeakAnno)

macs_peak <- system.file("extdata", "MACS_peaks.xls", 
                         package = "ChIPpeakAnno")
macs_peak_gr <- toGRanges(macs_peak, format = "MACS")
head(macs_peak_gr, n = 2)

## GRanges object with 2 ranges and 6 metadata columns:
##       seqnames            ranges strand |    length    summit      tags
##          <Rle>         <IRanges>  <Rle> | <integer> <integer> <integer>
##   X01     chr1 25323511-25324015      * |       505       252        45
##   X02     chr1 25362685-25362997      * |       313       211        33
##       -10*log10(pvalue) fold_enrichment       FDR
##               <numeric>       <numeric> <numeric>
##   X01             59.17           17.01       5.8
##   X02             60.63           22.41       4.2
##   -------
##   seqinfo: 12 sequences from an unspecified genome; no seqlengths

Prepare annotations

This section demonstrates how to prepare annotation file from Ensembl package. The Ensembl package is generally favored, as it offers more comprehensive and superior annotation. For more information, please refer to Section @ref(prepanno).

# Use Ensembl annotation package:
library(EnsDb.Hsapiens.v86)

ensembl.hs86.transcript <- transcripts(EnsDb.Hsapiens.v86)

Annotate peaks

# Use Ensembl annotation package:
macs_peak_ensembl <- annotatePeakInBatch(macs_peak_gr, 
                                         AnnotationData = ensembl.hs86.transcript)
head(macs_peak_ensembl, n = 2)

## GRanges object with 2 ranges and 15 metadata columns:
##                       seqnames            ranges strand |    length    summit
##                          <Rle>         <IRanges>  <Rle> | <integer> <integer>
##   X01.ENST00000564413     chr1 25323511-25324015      * |       505       252
##   X02.ENST00000487234     chr1 25362685-25362997      * |       313       211
##                            tags -10*log10(pvalue) fold_enrichment       FDR
##                       <integer>         <numeric>       <numeric> <numeric>
##   X01.ENST00000564413        45             59.17           17.01       5.8
##   X02.ENST00000487234        33             60.63           22.41       4.2
##                              peak         feature start_position end_position
##                       <character>     <character>      <integer>    <integer>
##   X01.ENST00000564413         X01 ENST00000564413       25335846     25337218
##   X02.ENST00000487234         X02 ENST00000487234       25342592     25351654
##                       feature_strand insideFeature distancetoFeature
##                          <character>   <character>         <numeric>
##   X01.ENST00000564413              -    downstream             13707
##   X02.ENST00000487234              +    downstream             20093
##                       shortestDistance fromOverlappingOrNearest
##                              <integer>              <character>
##   X01.ENST00000564413            11831          NearestLocation
##   X02.ENST00000487234            11031          NearestLocation
##   -------
##   seqinfo: 12 sequences from an unspecified genome; no seqlengths

For more regarding the metadata columns, please refer to the Section @ref(findnearest). Of note, the metadata does not come with “gene symbol” column. We can add it using the addGeneIDs function as below.

Add gene symbols

To add gene symbols, we can leverage either the organism annotation package, namely org.Hs.eg.db, or the bioRmart package. The latter offers the most current information but may have a slower performance as it requires querying the BioMart databases in real time. Once the required package is loaded, you can use the addGeneIDs function to add gene symbols. The following example utilizes the bioRmart method.

library(biomaRt)
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                dataset = "hsapiens_gene_ensembl")

# Use Ensembl annotation package:
macs_peak_ensembl <- addGeneIDs(annotatedPeak = macs_peak_ensembl, 
                                mart = mart,
                                feature_id_type = "ensembl_transcript_id",
                                IDs2Add = "hgnc_symbol")
head(macs_peak_ensembl, n = 2)

## GRanges object with 2 ranges and 16 metadata columns:
##                       seqnames            ranges strand |    length    summit
##                          <Rle>         <IRanges>  <Rle> | <integer> <integer>
##   X01.ENST00000564413     chr1 25323511-25324015      * |       505       252
##   X02.ENST00000487234     chr1 25362685-25362997      * |       313       211
##                            tags -10*log10(pvalue) fold_enrichment       FDR
##                       <integer>         <numeric>       <numeric> <numeric>
##   X01.ENST00000564413        45             59.17           17.01       5.8
##   X02.ENST00000487234        33             60.63           22.41       4.2
##                              peak         feature start_position end_position
##                       <character>     <character>      <integer>    <integer>
##   X01.ENST00000564413         X01 ENST00000564413       25335846     25337218
##   X02.ENST00000487234         X02 ENST00000487234       25342592     25351654
##                       feature_strand insideFeature distancetoFeature
##                          <character>   <character>         <numeric>
##   X01.ENST00000564413              -    downstream             13707
##   X02.ENST00000487234              +    downstream             20093
##                       shortestDistance fromOverlappingOrNearest hgnc_symbol
##                              <integer>              <character> <character>
##   X01.ENST00000564413            11831          NearestLocation       RSRP1
##   X02.ENST00000487234            11031          NearestLocation     TMEM50A
##   -------
##   seqinfo: 12 sequences from an unspecified genome; no seqlengths

Now, a hgnc_symbol metadata column has been successfully inserted to the resulting GRanges object. It is crucial to specify the correct feature_id_type and IDs2Add for the function to work properly. By default, feature_id_type = "ensembl_gene_id".

The following sections demonstration each topic in more details.

Obtain final peak set from replicates

For experiments with two or more biological replicates, there are several strategies² to retrieve a final peak set:

• Take the intersection of peaks:
  • Most conservative, reduces noise by considering the most robust and consistently enriched peaks
  • May result in small number of peaks with narrower widths
• Take the union of all peaks: may lead to large number of peaks
  • Least conservative, provides a more comprehensive view of the binding landscape
  • Useful for identifying condition-specific or rare binding events
• Merge the overlapping peaks: may lead to peaks with broader widths
  • Can be achieved with ChIPpeakAnno::findOverlapsOfPeaks
• Use peak coverage voting to determine a final set:
  • Can be achieved with consensusSeekeR
• Statistical model based approaches that leverage peak coverage and/or read depth information

Investigators may choose one of these strategies based on their research questions and the quality of their data.

Assess the concordance of peak replicates

Concordant peaks refer to the peaks that come from different biological replicates and exhibit a high degree of overlap in their genomic coordinates and signal intensity. The presence of concordant peaks indicates that the observed signal is likely to be true positive rather than arising from random noise or technical artifacts. You can evaluate the concordance of replicates by visualizing the peak overlappings.

# For two samples:
data(peaks1)
data(peaks2)
head(peaks1, n = 2)

## GRanges object with 2 ranges and 0 metadata columns:
##         seqnames          ranges strand
##            <Rle>       <IRanges>  <Rle>
##   Site1        1   967654-967754      +
##   Site2        2 2010897-2010997      +
##   -------
##   seqinfo: 6 sequences from an unspecified genome; no seqlengths

head(peaks2, n = 2)

## GRanges object with 2 ranges and 0 metadata columns:
##      seqnames          ranges strand
##         <Rle>       <IRanges>  <Rle>
##   t1        1   967659-967869      +
##   t2        2 2010898-2011108      +
##   -------
##   seqinfo: 6 sequences from an unspecified genome; no seqlengths

ol <- findOverlapsOfPeaks(peaks1, peaks2)
names(ol)

## [1] "venn_cnt"           "peaklist"           "uniquePeaks"       
## [4] "mergedPeaks"        "peaksInMergedPeaks" "overlappingPeaks"  
## [7] "all.peaks"

# For three samples:
data(peaks3)
ol2 <- findOverlapsOfPeaks(peaks1, peaks2, peaks3,
                           connectedPeaks = "min")

# Find peaks with at least 90% overlap:
ol3 <- findOverlapsOfPeaks(peaks1, peaks2, minoverlap = 0.9)

# For peaks1:
length(ol$peaklist[["peaks1"]])

## [1] 1

# For peaks2:
length(ol$peaklist[["peaks2"]])

## [1] 6

Visualize the overlaps using Venn diagram

The output from findOverlapsOfPeaks can be directly fed to makeVennDiagram to draw a Venn diagram and evaluate the degree of overlap between peak sets. Additionally, the makeVennDiagram function also calculates a P-value, which indicates whether the observed overlap is statistically significant or not.

# For two samples:
venn <- makeVennDiagram(ol, totalTest = 100,
                        fill = c("#009E73", "#F0E442"),
                        col = c("#D55E00", "#0072B2"),
                        cat.col = c("#D55E00", "#0072B2"))

Venn diagram showing the number of overlapping peaks for two samples

# For three samples:
venn2 <- makeVennDiagram(ol2, totalTest = 100,
                         fill = c("#CC79A7", "#56B4E9", "#F0E442"),
                         col = c("#D55E00", "#0072B2", "#E69F00"),
                         cat.col = c("#D55E00", "#0072B2", "#E69F00"))

Venn diagram showing the number of overlapping peaks for two samples

# For peaks overlap at least 90%:
venn3 <- makeVennDiagram(ol3, totalTest = 100,
                         fill = c("#009E73", "#F0E442"),
                         col = c("#D55E00", "#0072B2"),
                         cat.col = c("#D55E00", "#0072B2"))

Venn diagram showing the number of overlapping peaks for two samples

The parameter totalTest refers to the total number of potential peaks that are considered in the hypergeometric test. It should be set to a value larger than the largest number of peaks in the replicates. Refer to Section @ref(hypertest) on how to choose totalTest. Since P-value is sensitive to the choice of totalTest, we recommend using permutation test, implemented as peakPermTest in Biocpkg("ChIPpeakAnno"). For details, go to Section @ref(permtest).

if (!library(Vennerable)) {
  install.packages("Vennerable", repos="http://R-Forge.R-project.org")
  library(Vennerable)
}

n <- which(colnames(ol$vennCounts) == "Counts") - 1
set_names <- colnames(ol$vennCounts)[1:n]
weight <- ol$vennCounts[, "Counts"]
names(weight) <- apply(ol$vennCounts[, 1:n], 1, base::paste, collapse = "")
v <- Venn(SetNames = set_names, Weight = weight)
plot(v)

names(ol$overlappingPeaks)

## [1] "peaks1///peaks2"

dim(ol$overlappingPeaks[["peaks1///peaks2"]])

## [1] 12 14

ol$overlappingPeaks[["peaks1///peaks2"]][1:2, ]

##          peaks1 seqnames   start     end width strand peaks2 seqnames   start
## Site1_t1  Site1        1  967654  967754   101      +     t1        1  967659
## Site2_t2  Site2        2 2010897 2010997   101      +     t2        2 2010898
##              end width strand overlapFeature shortestDistance
## Site1_t1  967869   211      +   overlapStart                5
## Site2_t2 2011108   211      +   overlapStart                1

unique(ol$overlappingPeaks[["peaks1///peaks2"]][["overlapFeature"]])

## [1] "overlapStart"   "inside"         "overlapEnd"     "includeFeature"

pie1(table(ol$overlappingPeaks[["peaks1///peaks2"]]$overlapFeature))

Commonly used annotations

Popular annotation files come from two sources and are typically stored in tab-delimited formats, such as GTF or BED.

In addition, the UCSC Genome Browser team provides processed annotations based on the above resources that can be visualized easily with the browser. For human and mouse, there are four GTF files provided respectively.

We would recommend using EnsDb annotations because they are usually more comprehensive and consistent.

Bioconductor supported TxDb and EnsDb

In Bioconductor, the tab-delimited files are often converted into TxDb or EnsDb object. Both can be created with the tab delimited files mentioned above. While EnsDb is tailored to Ensembl annotations and contains additional information such as gene_name, symbol, and gene_biotype, the TxDb counterpart is typically created from RefSeq or UCSC Genome Browser hosted annotations.

There is a comprehensive list of pre-built Bioconductor maintained TxDb and EnsDb packages such as TxDb.Hsapiens.UCSC.hg38.knownGene and EnsDb.Hsapiens.v86. The list is regularly updated and can be found here.

Obtain EnsDb and TxDb from AnnotationHub

Users can also retrieve annotations using the AnnotationHub package. By default, AnnotationHub uses a snapshot that matches the version of Bioconductor being used, so it may be slightly more up-to-date compared to the pre-built packages. Additionally, users have the option to switch to an earlier version³. Below are examples of how to obtain annotations using AnnotationHub.

library(AnnotationHub)

ah <- AnnotationHub()

# Obtain EnsDb:
EnsDb_Mmusculus_all <- query(ah, pattern = c("Mus musculus", "EnsDb"))
head(EnsDb_Mmusculus_all, n = 2)

## AnnotationHub with 2 records
## # snapshotDate(): 2024-10-28
## # $dataprovider: Ensembl
## # $species: Mus musculus
## # $rdataclass: EnsDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH53222"]]' 
## 
##             title                            
##   AH53222 | Ensembl 87 EnsDb for Mus Musculus
##   AH53726 | Ensembl 88 EnsDb for Mus Musculus

EnsDb_Mmusculus <- EnsDb_Mmusculus_all[["AH53222"]]
class(EnsDb_Mmusculus)

## [1] "EnsDb"
## attr(,"package")
## [1] "ensembldb"

# Obtain TxDb:
TxDb_Mmusculus_all <- query(ah, pattern = c("Mus musculus", "TxDb"))
head(TxDb_Mmusculus_all, n = 2)

## AnnotationHub with 2 records
## # snapshotDate(): 2024-10-28
## # $dataprovider: UCSC
## # $species: Mus musculus
## # $rdataclass: TxDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH52263"]]' 
## 
##             title                                    
##   AH52263 | TxDb.Mmusculus.UCSC.mm10.ensGene.sqlite  
##   AH52264 | TxDb.Mmusculus.UCSC.mm10.knownGene.sqlite

TxDb_Mmusculus <- TxDb_Mmusculus_all[["AH52264"]]
class(TxDb_Mmusculus)

## [1] "TxDb"
## attr(,"package")
## [1] "GenomicFeatures"

Build custom EnsDb and TxDb

To create the most up-to-date or custom TxDb and EnsDb objects, you can utilize functions such as makeTxDbFromUCSC, makeTxDbFromEnsembl, and ensDbFromGRanges from the GenomicFeatures⁴ package and the Ensembldb⁵ package.

Use biomaRt

The biomaRt package offers a convenient interface to the BioMart databases that are prominently maintained by Ensembl. By querying biomaRt, you can access the latest available annotations from Ensembl. Check this vignette for details.

ChIPpeakAnno provides a helper function called getAnnotation, which simplifies the retrieval of desired annotations by leveraging biomaRt. Here is an example:

library(biomaRt)

listMarts()
head(listDatasets(useMart("ENSEMBL_MART_ENSEMBL")), n = 2)
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                dataset = "mmusculus_gene_ensembl")
anno_from_biomart <- getAnnotation(mart, 
                                   featureType = "transcript")
head(anno_from_biomart, n = 2)

Use toGRanges function

ChIPpeakAnno offers a helper function called toGRanges, which can convert annotations from various formats including GFF, BED, CSV, TxDb, and EnsDb into GRanges. Below is an example.

library(EnsDb.Hsapiens.v86)

anno_ensdb_transcript <- toGRanges(EnsDb.Hsapiens.v86, 
                                   feature = "transcript")
head(anno_ensdb_transcript, n = 2)

## GRanges object with 2 ranges and 3 metadata columns:
##                   seqnames      ranges strand |           tx_id         gene_id
##                      <Rle>   <IRanges>  <Rle> |     <character>     <character>
##   ENST00000456328     chr1 11869-14409      + | ENST00000456328 ENSG00000223972
##   ENST00000450305     chr1 12010-13670      + | ENST00000450305 ENSG00000223972
##                     gene_name
##                   <character>
##   ENST00000456328     DDX11L1
##   ENST00000450305     DDX11L1
##   -------
##   seqinfo: 357 sequences (1 circular) from 2 genomes (hg38, GRCh38)

Similarly, set feature = "gene" to retrieve gene annotations in GRanges format.

anno_ensdb_transcript <- transcripts(EnsDb.Hsapiens.v86)

head(anno_ensdb_transcript, n = 2)

## GRanges object with 2 ranges and 6 metadata columns:
##                   seqnames      ranges strand |           tx_id
##                      <Rle>   <IRanges>  <Rle> |     <character>
##   ENST00000456328        1 11869-14409      + | ENST00000456328
##   ENST00000450305        1 12010-13670      + | ENST00000450305
##                               tx_biotype tx_cds_seq_start tx_cds_seq_end
##                              <character>        <integer>      <integer>
##   ENST00000456328   processed_transcript             <NA>           <NA>
##   ENST00000450305 transcribed_unproces..             <NA>           <NA>
##                           gene_id         tx_name
##                       <character>     <character>
##   ENST00000456328 ENSG00000223972 ENST00000456328
##   ENST00000450305 ENSG00000223972 ENST00000450305
##   -------
##   seqinfo: 357 sequences (1 circular) from GRCh38 genome

Plot peak distributions relative to genomic features

The binOverFeature function can be used to plot the distribution of peak counts relative to a specific genomic feature. The following example shows the distribution of peaks in the macs_peak_gr2 dataset relative to the gene feature (transcription start site (TSS)).

library(TxDb.Hsapiens.UCSC.hg38.knownGene)

annotation_data <- transcripts(TxDb.Hsapiens.UCSC.hg38.knownGene)
binOverFeature(macs_peak_gr2, 
               featureSite = "FeatureStart",
               nbins = 20,
               annotationData = annotation_data,
               xlab = "peak distance from TSS (bp)", 
               ylab = "peak count", 
               main = "Distribution of aggregated peak numbers around TSS")

Peak count distribution around TSSs

By default, featureSite = "FeatureStart" meaning that distance is calculated as peak to feature start (i.e. TSS for transcript). The plot above, describing the distribution of peaks around TSS, exhibits a signal summit around TSS, a characteristic pattern observed in peaks obtained from transcription factor binding experiments.

You can also plot peak distribution over multiple genomic features including exon, intron, enhancer, proximal promoter, 5’ UTR and 3’ UTR in a single bar graph using assignChromosomeRegion.

chromosome_region <- assignChromosomeRegion(macs_peak_gr2,
                                            TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene,
                                            nucleotideLevel = FALSE,
                                            precedence=c("Promoters",
                                                         "immediateDownstream", 
                                                         "fiveUTRs", 
                                                         "threeUTRs",
                                                         "Exons", 
                                                         "Introns"))

# optional helper function to rotate x-axis labels for barplot(): 
# ref: https://stackoverflow.com/questions/10286473/rotating-x-axis-labels-in-r-for-barplot
rotate_x <- function(data, rot_angle) {
  plt <- barplot(data, xaxt = "n")
  text(plt, par("usr")[3], 
       labels = names(data), 
       srt = rot_angle, adj = c(1.1,1.1), 
       xpd = TRUE, cex = 0.6)
}

rotate_x(chromosome_region[["percentage"]], 45)

Bar graph showing peak distributions over different genomic features

By default, nucleotideLevel = FALSE meaning that peaks are treated as ranges when determining overlaps with genomic features. If a peak intersects with multiple features, the feature assignment is determined by the order specified in the precedence argument. If precedence is not set, counts for each overlapping feature will be incremented. Otherwise, if nucleotideLevel = TRUE, the summit of the peak (a single position, not suitable for broad peaks) will be used when determining overlaps.

Plot peak distributions over different feature levels

In addition to inspecting the peak enrichment pattern by plotting the distribution against genomic features, users can plot distributions over different feature levels, each containing multiple categories, using the genomicElementDistribution function.

• Transcript Level: “geneLevel”
  • “Promoter”
  • “Gene body”
  • “Downstream”
  • “Distal Intergenic”
• Exon Level: “Exons”
  • “5’ UTR”
  • “CDS”
  • “3’ URT”
  • “Other exon”
• Exon/Intron/Intergenic: “ExonIntron”
  • “Exon”
  • “Intron”
  • “Intergenic”

Please note that peaks can be classified into multiple categories from different levels, leading to the total percentage of annotated features being greater than 100%. At each level, since a peak spans a genomic range, it may overlap with multiple categories of features. In such cases, by default nucleotideLevel = FALSE, which means that the precedence is determined by the order listed in the labels argument.

The genomicElementDistribution function accepts either a single peak object or a list of peak objects as input. If a single peak object if provided, a pie chart will be created; if a list of peak objects is provided, a bar graph will be created.

Pie graph with one peak set

genomicElementDistribution(macs_peak_gr1, 
                           TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)

Pie graph showing peak distributions over features at different levels

As can be seen, a significant number of peaks originate from promoter regions.

Bar graph with a list of peak sets

macs_peaks <- GRangesList(rep1 = macs_peak_gr1,
                          rep2 = macs_peak_gr2)
genomicElementDistribution(macs_peaks, 
                           TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)

Bar graph showing peak distributions over features at different levels

The consistent patterns from “rep1” and “rep2” indicate a high correlation between them.

Plot peak overlaps for multiple features

We can create an UpSet plot to view peak overlaps across multiple genomic features.

library(UpSetR)

res <- genomicElementUpSetR(macs_peak_gr1,
                            TxDb.Hsapiens.UCSC.hg38.knownGene)
upset(res[["plotData"]], 
      nsets = length(colnames(res$plotData)), 
      nintersects = NA)

UpSet plot showing the overlappings of peak distributions among different genomic features

UpSet plot can be considered as a “high dimensional Venn diagram” that allows for the visualization of overlaps for multiple sets. For example, in the above plot, it is evident that the feature set “gene body” (from the “Transcript Level”) and the feature set “intron” (from the “Exon/Intron/Intergenic”) share the highest number of common peaks. Type ?genomicElementUpSetR for the definition of distal_promoter.

Peak-centric method

You can assign the nearest or overlapping features to your peaks by setting the output option to the following values:

• output: criteria for assigning features
  • “nearestLocation”: (default) output the features that are nearest to the peaks based on the absolute value of the difference between the peak location and the feature location, as calculated by abs(PeakLocForDistance - FeatureLocForDistance)
  • “overlapping”: output features that overlap with the peaks within the maximum distance specified by the maxgap parameter
  • “both”: output both the nearest feature and overlapping features that are not the nearest
  • “shortestDistance”: output the features that have the shortest distance to the peaks; the “shortest distance” is determined from either end of the feature to either end of the peak

Other relevant parameters:

• PeakLocForDistance: location of the peak for distance calculation
  • “middle”: (recommended) use the center of the peak
  • “start”: (default) use the peak edge that is smaller in coordinate
  • “end”: use the peak edge that is larger in coordinate
  • “endMinusStart”: use “end” defined above if feature is on the plus strand, and “start” if minus strand
• FeatureLocForDistance: location of the feature for distance calculation
  • “middle”: use the center of the feature
  • “start”: use the feature edge that is smaller in coordinate
  • “end”: use the feature edge that is large in coordinate
  • “TSS”: (default) use “start” defined above if feature is on the plus strand and “end” if on the minus strand
  • “geneEnd”: use “end” defined above if feature is on the plus strand and “start” if on the minus strand
• maxgap: the maximum gap that is allowed between the peak and the feature ranges for them to be considered as overlapping; it is defined as the number of bases that separates the two ranges, the gap between two adjacent ranges is 0

How does peak-centric annotation method work

Feature-centric method

Peaks can also be annotated based on their relative locations to nearby features. For example, by setting output = "upstream", peaks will be annotated to features that they are located upstream of.

Relative peak-to-feature location

You can use the following options to specify the desired relative locations of the peaks to the features.

• output: criteria for assigning features
  • “upstream”: output the feature if the peak resides upstream of it
  • “upstrem&inside”: output the feature if the peak resides upstream of it or contained within it
  • “inside”: output the feature if the peak is completely contained within it
  • “inside&downstream”: output the feature if the peak resides within it or downstream of it
  • “downstream”: output the feature if the peak resides downstream of it
  • “upstreamORdownstream”: output the feature if the peak resides upstream or downstream of it

Other relevant parameter:

• maxgap: allowed gap when determining overlap between two ranges; will be ignored if output = "inside"

The following diagram illustrates how peaks are annotated using the feature-centric method.

How does feature-centric annotation method work

When using the feature-centric method, the annotatePeakInBatch function will output the feature as long as the peak range overlaps with the target region, except when output = "inside". In this case, the entire peak range must be completely contained within the feature range to output it.

More customization with bindingRegion

The bindingRegion parameter, which refers to the regions that the target TF binds to, adds further flexibility to the feature-centric method when defining the target region. For example, if you would like to annotate peaks to features where the peaks are located within specific regions relative to them, such as from 2000bp upstream to 500bp downstream of TSS (where most promoters are found), you can set output = "overlapping", FeatureLocForDistance = "TSS", and bindingRegion = c(-2000, 500). Once bindingRegion is specified, the maxgap will be ignored.

The following diagram demonstrates several examples of how to use bindingRegion to define target regions.

How to define target region with bindingRegion

Bi-directional promoters are genomic regions that are located upstream of the TSS of two adjacent genes that are transcribed in opposite directions (Adachi and Lieber 2002). Those promoters typically regulate the expression of both genes. To annotate peaks located in such regions, you can use output = "overlapping", FeatureLocForDistance = "TSS", and bindingRegion = c(-5000, 3000). By default, select = "all", meaning that all overlapping features will be outputted. If you want to annotate peaks to features with the nearest bi-directional promoters, you can use output = "nearestBiDirectionalPromoters" along with bindingRegion. In this setting, at most one feature will be reported from each direction. When using output = "nearestBiDirectionalPromoters", both maxgap and FeatureLocForDistance will be ignored.

To illustrate, we can use the myPeakList provided by ChIPpeakAnno. As different major genome releases (e.g. hg19 vs hg38) may have variations in feature coordinates, it is highly recommended to use an annotation file that matches the genome version used when generating your peak file. In the case of myPeakList, the peaks were originally called against hg18, so it is necessary to use a matching annotation file created with hg18. Alternatively, you can use the rtracklayer::liftOver() function to convert myPeakList to hg38 coordinates. For step-by-step instructions, refer to “Step3” in Section @ref(custompool).

Example 1: find the nearest features

We can annotate the peaks by assigning the nearby features to them. This can be achieved by setting output = "nearestLocation". When this option is used, the results may include “overlapping” features as long as they are the nearest ones to the peaks.

library(TxDb.Hsapiens.UCSC.hg18.knownGene)

data(myPeakList)
peak_to_anno <- myPeakList[1:100]
anno_data <- transcripts(TxDb.Hsapiens.UCSC.hg18.knownGene)

annotated_peak <- annotatePeakInBatch(peak_to_anno, 
                                      output = "nearestLocation",
                                      PeakLocForDistance = "start",
                                      FeatureLocForDistance = "TSS",
                                      AnnotationData = anno_data)
head(annotated_peak, n = 2)

## GRanges object with 2 ranges and 9 metadata columns:
##                      seqnames        ranges strand |         peak     feature
##                         <Rle>     <IRanges>  <Rle> |  <character> <character>
##   X1_93_556427.ann16     chr1 556660-556760      * | X1_93_556427       ann16
##   X1_41_559455.ann19     chr1 559774-559874      * | X1_41_559455       ann19
##                      start_position end_position feature_strand insideFeature
##                           <integer>    <integer>    <character>   <character>
##   X1_93_556427.ann16         556325       557910              +        inside
##   X1_41_559455.ann19         559396       560212              +        inside
##                      distancetoFeature shortestDistance
##                              <numeric>        <integer>
##   X1_93_556427.ann16               335              335
##   X1_41_559455.ann19               378              338
##                      fromOverlappingOrNearest
##                                   <character>
##   X1_93_556427.ann16          NearestLocation
##   X1_41_559455.ann19          NearestLocation
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

Here is a breakdown of the options:

• output = "nearestLocation": (default)
  • output the nearest features calculated as abs(PeakLocForDistance - FeatureLocForDistance). The output can consist of both “strictly nearest features (non-overlapping)” and “overlapping features” as long as it is the nearest
• PeakLocForDistance = "start" (default):
  • “start” means using the 5’ end (relative to plus strand) of the peak to calculate the distance to features
• FeatureLocForDistance = "TSS" (default):
  • “TSS” means using the 5’ end (relative to plus strand) for features on the plus strand and use the 3’ end (also relative to plus strand) for features on the minus strand to calculate the distance to features

Here is a breakdown of the resulting metadata columns:

• peak: id of the peak
• feature: id of the feature such as Ensembl gene ID
• start_position, end_position, feature_strand: feature coordinates
• insideFeature: relative location of the peak to the feature
  • “upstream”: peak resides upstream of the feature
  • “downstream”: peak resides downstream of the feature
  • “inside”: peak resides inside the feature
  • “overlapStart”: peak overlaps with the start of the feature
  • “overlapEnd”: peak overlaps with the end of the feature
  • “includeFeature”: peak includes the feature entirely
• distancetoFeature: peak-to-feature distance as calculated with abs(PeakLocForDistance - FeatureLocForDistance)
• shortestDistance: the shortest distance from either end of the peak to either end the feature
• fromOverlappingOrNearest: relevant only when output = "both"; “nearestLocation” indicates that the feature is the closest to the peak, can overlap with the peak; “Overlapping” means that the feature overlaps with the peak and it is not the nearest; when output = "nearestLocation", this column should consist of only “NearestLocation”

Example 2: find the nearest and overlapping features

In addition, annotatePeakInBatch can also report both the nearest features and overlapping features by setting output = "both" and the maxgap parameter. For example, the following command outputs the nearest features plus all overlapping features that are within 100bp away.

annotated_peak <- annotatePeakInBatch(peak_to_anno, 
                                      AnnotationData = anno_data,
                                      output = "both",
                                      maxgap = 100)
head(annotated_peak, n = 4)

## GRanges object with 4 ranges and 9 metadata columns:
##                        seqnames        ranges strand |         peak     feature
##                           <Rle>     <IRanges>  <Rle> |  <character> <character>
##     X1_93_556427.ann16     chr1 556660-556760      * | X1_93_556427       ann16
##   X1_93_556427.ann3352     chr1 556660-556760      * | X1_93_556427     ann3352
##     X1_41_559455.ann19     chr1 559774-559874      * | X1_41_559455       ann19
##   X1_41_559455.ann3352     chr1 559774-559874      * | X1_41_559455     ann3352
##                        start_position end_position feature_strand insideFeature
##                             <integer>    <integer>    <character>   <character>
##     X1_93_556427.ann16         556325       557910              +        inside
##   X1_93_556427.ann3352          79158       735413              -        inside
##     X1_41_559455.ann19         559396       560212              +        inside
##   X1_41_559455.ann3352          79158       735413              -        inside
##                        distancetoFeature shortestDistance
##                                <numeric>        <integer>
##     X1_93_556427.ann16               335              335
##   X1_93_556427.ann3352            178753           178653
##     X1_41_559455.ann19               378              338
##   X1_41_559455.ann3352            175639           175539
##                        fromOverlappingOrNearest
##                                     <character>
##     X1_93_556427.ann16          NearestLocation
##   X1_93_556427.ann3352              Overlapping
##     X1_41_559455.ann19          NearestLocation
##   X1_41_559455.ann3352              Overlapping
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

Now, the fromOverlappingOrNearest column consists of both “NearestLocation” and “Overlapping” categories.

Visualize peak-to-feature distances

The relative location distribution of peak-to-feature can be visualized using the information stored in the insideFeature metadata column.

pie1(table(annotated_peak$insideFeature))

Peak distribution around features

Use custom annotation data

You also have the option to create and pass user-defined feature coordinates in GRanges format as annotationData. For example, if you have a list of transcript factor binding sites obtained by literature mining and would like to use it to annotate your peaks. you can convert the TF binding site coordinates into a GRanges object and then pass that object into annotatePeakInBatch.

TF_binding_sites <- GRanges(seqnames = c("1", "2", "3", "4", "5", "6", "1", 
                                         "2", "3", "4", "5", "6", "6", "6", 
                                         "6", "6", "5"),
                            ranges = IRanges(start = c(967659, 2010898, 2496700, 
                                                       3075866, 3123260, 3857500,
                                                       96765, 201089, 249670, 
                                                       307586, 312326, 385750, 
                                                       1549800, 1554400, 1565000, 
                                                       1569400, 167888600),
                                             end = c(967869, 2011108, 2496920, 
                                                     3076166,3123470, 3857780,
                                                     96985, 201299, 249890, 307796, 
                                                     312586, 385960, 1550599, 
                                                     1560799, 1565399, 1571199, 
                                                     167888999),
                                             names = paste("t", 1:17, sep = "")),
                            strand = c("*", "*", "*", "*", "*", "*", "*", "*", "*", 
                                       "*", "*", "*", "*", "*", "*", "*", "*"))

annotated_peak2 <- annotatePeakInBatch(peaks1, AnnotationData = TF_binding_sites)
head(annotated_peak2, n = 2)

## GRanges object with 2 ranges and 9 metadata columns:
##            seqnames          ranges strand |        peak     feature
##               <Rle>       <IRanges>  <Rle> | <character> <character>
##   Site1.t1        1   967654-967754      + |       Site1          t1
##   Site2.t2        2 2010897-2010997      + |       Site2          t2
##            start_position end_position feature_strand insideFeature
##                 <integer>    <integer>    <character>   <character>
##   Site1.t1         967659       967869              *  overlapStart
##   Site2.t2        2010898      2011108              *  overlapStart
##            distancetoFeature shortestDistance fromOverlappingOrNearest
##                    <numeric>        <integer>              <character>
##   Site1.t1                -5                5          NearestLocation
##   Site2.t2                -1                1          NearestLocation
##   -------
##   seqinfo: 6 sequences from an unspecified genome; no seqlengths

Another example of using user-defined AnnotationData is to annotate peaks by promoters. A promoter is typically defined as the DNA sequence located immediately upstream of the TSS of a gene. The specific size of a promoter can vary depending on the gene, its regulatory complexity, and the species being studied. In practice, the promoter region can be defined as 2000bp upstream and 500bp downstream from the TSS. To prepare a custom annotation file containing only promoters, users can leverage the accessor function promoters. Similar results can be obtained using the feature-centric approach mentioned in Section @ref(featurecentric).

promoter_regions <- promoters(TxDb.Hsapiens.UCSC.hg18.knownGene, 
                              upstream = 2000, downstream = 500)
head(promoter_regions, n = 2)

## GRanges object with 2 ranges and 2 metadata columns:
##              seqnames    ranges strand |     tx_id     tx_name
##                 <Rle> <IRanges>  <Rle> | <integer> <character>
##   uc001aaa.2     chr1 -884-1615      + |         1  uc001aaa.2
##   uc009vip.1     chr1 -884-1615      + |         2  uc009vip.1
##   -------
##   seqinfo: 49 sequences (1 circular) from hg18 genome

annotated_peak3 <- annotatePeakInBatch(peak_to_anno, 
                                       AnnotationData = promoter_regions)
head(annotated_peak3, n = 2)

## GRanges object with 2 ranges and 9 metadata columns:
##                           seqnames        ranges strand |         peak
##                              <Rle>     <IRanges>  <Rle> |  <character>
##   X1_93_556427.uc001abb.1     chr1 556660-556760      * | X1_93_556427
##   X1_41_559455.uc001abc.1     chr1 559774-559874      * | X1_41_559455
##                               feature start_position end_position
##                           <character>      <integer>    <integer>
##   X1_93_556427.uc001abb.1  uc001abb.1         556707       559206
##   X1_41_559455.uc001abc.1  uc001abc.1         557396       559895
##                           feature_strand insideFeature distancetoFeature
##                              <character>   <character>         <numeric>
##   X1_93_556427.uc001abb.1              +  overlapStart               -47
##   X1_41_559455.uc001abc.1              +        inside              2378
##                           shortestDistance fromOverlappingOrNearest
##                                  <integer>              <character>
##   X1_93_556427.uc001abb.1               47          NearestLocation
##   X1_41_559455.uc001abc.1               21          NearestLocation
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

Example1: find gene symbols for a vector of entrez IDs

The following example demonstrates how to use the addGeneIDs function to find gene symbols for a vector of entrez IDs using OrgDb. Note that the "org.Hs.eg.db" package name must be quoted.

library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

ucsc.hg38.knownGene <- genes(TxDb.Hsapiens.UCSC.hg38.knownGene)
entrez_ids <- head(ucsc.hg38.knownGene$gene_id, n = 10)
print(entrez_ids)

##  [1] "1"         "10"        "100"       "1000"      "100008586" "100008587"
##  [7] "100009613" "100009667" "100009676" "10001"

res <- addGeneIDs(entrez_ids, 
                  orgAnn = "org.Hs.eg.db", 
                  feature_id_type = "entrez_id",
                  IDs2Add = "symbol")
head(res, n = 3)

##   entrez_id symbol
## 1         1   A1BG
## 2        10   NAT2
## 3       100    ADA

Example2: add gene symbols to annotated peaks

For this example, we annotate the macs_peak_gr2 (obtained in Section @ref(readinpeak)) using transcript information from TxDb, and add gene symbols to them.

txdb.hg38.transcript <- transcripts(TxDb.Hsapiens.UCSC.hg38.knownGene)
head(txdb.hg38.transcript, n = 4)

## GRanges object with 4 ranges and 2 metadata columns:
##       seqnames      ranges strand |     tx_id           tx_name
##          <Rle>   <IRanges>  <Rle> | <integer>       <character>
##   [1]     chr1 11869-14409      + |         1 ENST00000456328.2
##   [2]     chr1 12010-13670      + |         2 ENST00000450305.2
##   [3]     chr1 29554-31097      + |         3 ENST00000473358.1
##   [4]     chr1 30267-31109      + |         4 ENST00000469289.1
##   -------
##   seqinfo: 711 sequences (1 circular) from hg38 genome

head(names(txdb.hg38.transcript), n = 4)

## NULL

It appears that the txdb.hg38.transcript contains a metadata column called tx_name, which holds the Ensembl transcript IDs. However, the annotatePeakInBatch function requires this information to be stored as the names of the txdb.hg38.trancsript object (names(txdb.hg38.transcript)). This is necessary to generate annotated peaks that are compatible with the addGeneIDs function. To achieve this, we need to extract the transcript IDs and assign them as names. Below is how.

tr_id <- txdb.hg38.transcript$tx_name
tr_id <- sub("\\..*$", "", tr_id) # get rid of the trailing version number
names(txdb.hg38.transcript) <- tr_id
head(txdb.hg38.transcript, n = 4)

## GRanges object with 4 ranges and 2 metadata columns:
##                   seqnames      ranges strand |     tx_id           tx_name
##                      <Rle>   <IRanges>  <Rle> | <integer>       <character>
##   ENST00000456328     chr1 11869-14409      + |         1 ENST00000456328.2
##   ENST00000450305     chr1 12010-13670      + |         2 ENST00000450305.2
##   ENST00000473358     chr1 29554-31097      + |         3 ENST00000473358.1
##   ENST00000469289     chr1 30267-31109      + |         4 ENST00000469289.1
##   -------
##   seqinfo: 711 sequences (1 circular) from hg38 genome

Now, we can annotate macs_peak_gr2 with annotatePeakInBatch.

res <- annotatePeakInBatch(macs_peak_gr2, 
                           AnnotationData = txdb.hg38.transcript)
head(res, n = 2)

## GRanges object with 2 ranges and 10 metadata columns:
##                               seqnames      ranges strand |     score
##                                  <Rle>   <IRanges>  <Rle> | <numeric>
##   MACS_peak_1.ENST00000473358     chr1 28341-29610      * |         1
##   MACS_peak_2.ENST00000495576     chr1 90821-91234      * |         1
##                                      peak         feature start_position
##                               <character>     <character>      <integer>
##   MACS_peak_1.ENST00000473358 MACS_peak_1 ENST00000473358          29554
##   MACS_peak_2.ENST00000495576 MACS_peak_2 ENST00000495576          89551
##                               end_position feature_strand insideFeature
##                                  <integer>    <character>   <character>
##   MACS_peak_1.ENST00000473358        31097              +  overlapStart
##   MACS_peak_2.ENST00000495576        91105              -  overlapStart
##                               distancetoFeature shortestDistance
##                                       <numeric>        <integer>
##   MACS_peak_1.ENST00000473358             -1213               56
##   MACS_peak_2.ENST00000495576               284              129
##                               fromOverlappingOrNearest
##                                            <character>
##   MACS_peak_1.ENST00000473358          NearestLocation
##   MACS_peak_2.ENST00000495576          NearestLocation
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

As anticipated, the resulting feature column contains Ensembl transcript IDs. Next, we utilize the mart option to add gene symbols.

library(biomaRt)

mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                dataset = "hsapiens_gene_ensembl")
res <- addGeneIDs(res, 
                  mart = mart,
                  feature_id_type = "ensembl_transcript_id",
                  IDs2Add = "hgnc_symbol")
head(res, n = 2)

## GRanges object with 2 ranges and 11 metadata columns:
##                               seqnames      ranges strand |     score
##                                  <Rle>   <IRanges>  <Rle> | <numeric>
##   MACS_peak_1.ENST00000473358     chr1 28341-29610      * |         1
##   MACS_peak_2.ENST00000495576     chr1 90821-91234      * |         1
##                                      peak         feature start_position
##                               <character>     <character>      <integer>
##   MACS_peak_1.ENST00000473358 MACS_peak_1 ENST00000473358          29554
##   MACS_peak_2.ENST00000495576 MACS_peak_2 ENST00000495576          89551
##                               end_position feature_strand insideFeature
##                                  <integer>    <character>   <character>
##   MACS_peak_1.ENST00000473358        31097              +  overlapStart
##   MACS_peak_2.ENST00000495576        91105              -  overlapStart
##                               distancetoFeature shortestDistance
##                                       <numeric>        <integer>
##   MACS_peak_1.ENST00000473358             -1213               56
##   MACS_peak_2.ENST00000495576               284              129
##                               fromOverlappingOrNearest hgnc_symbol
##                                            <character> <character>
##   MACS_peak_1.ENST00000473358          NearestLocation MIR1302-2HG
##   MACS_peak_2.ENST00000495576          NearestLocation            
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Use listMarts to show available biomart and use listDatasets to show available dataset. Be aware that, unlike using OrgDb option, we must supply hgnc_symbol instead of symbol for the IDs2Add argument.

Obtain sequences surrounding the peaks

Here is an example of how to retrieve the peak sequences, including 20bp upstream and 20bp downstream, for the macs_peak_gr2 peaks obtained in Section @ref(readinpeak).

library(BSgenome.Hsapiens.UCSC.hg38)

sequence_around_peak <- getAllPeakSequence(macs_peak_gr2, 
                                           upstream = 20,
                                           downstream = 20, 
                                           genome = BSgenome.Hsapiens.UCSC.hg38)
head(sequence_around_peak, n = 2)

## GRanges object with 2 ranges and 4 metadata columns:
##               seqnames      ranges strand |     score  upstream downstream
##                  <Rle>   <IRanges>  <Rle> | <numeric> <numeric>  <numeric>
##   MACS_peak_1     chr1 28341-29610      * |         1        20         20
##   MACS_peak_2     chr1 90821-91234      * |         1        20         20
##                             sequence
##                          <character>
##   MACS_peak_1 CTGTAGTTGCTCATCTGAAG..
##   MACS_peak_2 GCTGCTGCTGTCTGTAGCTG..
##   -------
##   seqinfo: 1 sequence from an unspecified genome

The genome argument can accept either a BSgenome object or a Mart object. It is important to ensure that the genome version used matches the one used for creating the peak file. You can find a full list of available BSgenome objects on this site.

The following example demonstrates on how to use the mart option. It is slower compared to using a BSgenome because it queries the BioMart for annotations on the fly if AnnotationData is not set. Refer to Section @ref(usebiomart) for more on Mart.

library(biomaRt)

mart <- useMart(biomart="ENSEMBL_MART_ENSEMBL",
                dataset="hsapiens_gene_ensembl")

sequence_around_peak <- getAllPeakSequence(macs_peak_gr2[1], 
                                           upstream = 20,
                                           downstream = 20, 
                                           genome = mart)

To save the sequences into fasta, use the write2FASTA function.

write2FASTA(sequence_around_peak, file = "macs_peak_gr2.fa")

Discover consensus sequences (motifs) in the peaks

The oligoSummary function utilizes Markov models to ascertain if a motif is enriched in a set of sequences relative to the background. As a prerequisite, we must first calculate the frequencies of all combinations of short oligonucleotides of a specified length in the background (Leung, Marsh, and Speed 1996). This can be accomplished using the oligoFrequency function. In the example below, we aim to identify consensus sequences of length 6.

freqs <- oligoFrequency(BSgenome.Hsapiens.UCSC.hg38$chr1)
motif_summary <- oligoSummary(sequence_around_peak, 
                              oligoLength = 6,
                              MarkovOrder = 3,
                              freqs = freqs,
                              quickMotif = TRUE)

Here is a breakdown of the arguments:

• oligoLength: the length of the motifs to search for
• MarkovOrder: the order of a Markov chain, which determines how much the current state in the chain depends on previous states; simply put, a zero-order Markov chain does not consider any context or dependencies between adjacent elements, while higher-order Markov chains capture longer-range dependencies in sequences, making them more suitable for modeling complex sequence patterns
• freqs: the background motif frequencies used for the Markov model
• quickMotif: whether or not to generate the motif matrices

The resulting motif_summary is a list containing three elements:

• zscore: the Z-score of each oligonucleotide, which quantifies how many standard deviations away the observed count is from the expected count; a high positive score suggests that the motif is significantly over-represented in the peaks
• counts: the count numbe of each oligonucleotide
• motifs: a list of motif matrices when quickMotif = TRUE

We can use histogram (hist) to visualize the resulting Z-scores and labels top hits with the text function. Below, we label the name of the motif that has the highest Z-score.

zscore <- sort(motif_summary$zscore)
h <- hist(zscore, breaks = 100, main = "Histogram of Z-score")
text(x = zscore[length(zscore)], 
     y = h$counts[length(h$counts)] + 1, 
     labels = names(zscore[length(zscore)]), 
     adj = 0, 
     srt = 90)

Histogram showing Z-score distribution

To illustrate how the Z-score can be influenced by the enrichment level of various motifs, the following code utilizes simulated data. We first generate 5000 sequences with lengths ranging between 100 and 2000 nucleotides. Subsequently,we randomly introduce motif 1 into 10% of the sequences and motif 2 into 15% of the sequences.

set.seed(1)

# motifs to simulate
simulate_motif_1 <- "AATTTA"
simulate_motif_2 <- "TGCATG"

# sample 5000 sequences with lengths ranging from 100 to 2000 nucleotides
# randomly insert motif_1 to 10% of the sequences, and motif_2 to 15% of the sequences
simulation_seqs <- sapply(sample(c(1, 2, 0), 
                                 5000,
                                 prob = c(0.1, 0.15, 0.75),
                                 replace = TRUE), 
                           function(x) {
                             seq <- sample(c("A", "T", "C", "G"),
                                           sample.int(1900, 1) + 100, 
                                           replace = TRUE)
                             insert_pos <- sample.int(length(seq) - 6, 1)
                             if (x == 1) {
                               seq[insert_pos:(insert_pos + 5)] <- strsplit(simulate_motif_1, "")[[1]]
                             } else if (x == 2) {
                               seq[insert_pos:(insert_pos + 5)] <- strsplit(simulate_motif_2, "")[[1]]
                             }
                             paste(seq, collapse = "")
                           }
)

motif_summary_simu <- oligoSummary(simulation_seqs, 
                                   oligoLength = 6, 
                                   MarkovOrder = 3, 
                                   quickMotif = TRUE)
zscore_simu <- sort(motif_summary_simu$zscore, 
                    decreasing = TRUE)
h_simu <- hist(x = zscore_simu, 
               breaks = 100, 
               main = "Histogram of Z-score for simulation data")
text(x = zscore_simu[1:2],  
     y = rep(5, 2), 
     labels = names(zscore_simu[1:2]), 
     adj = 0, 
     srt = 90)

Histogram showing Z-score distribution for simulation data

As evident from the simulation results, the higher the probability of the motif, the larger the resulting Z-score. In addition, you can visualize the motif using the motifStack package.

library(motifStack)

pfm <- new("pfm", mat = motif_summary_simu$motifs[[1]],
           name = "sample motif 1")
motifStack(pfm)

To loop through each element in motif_summary$motifs, we can use the mapply function.

pfms <- mapply(function(motif, id) { new("pfm", mat = motif, name = as.character(id)) },
               motif_summary$motifs,
               seq_along(motif_summary$motifs))

motifStack(pfms[[1]])

Plot showing the first motif detected

Scan pre-defined sequence patterns in the peaks

If you have a list of motifs (sequence patterns), the summarizePatternInPeaks function can be utilized to determine whether they are present in the peaks.

example_pattern_file <- system.file("extdata/examplePattern.fa",
                                    package = "ChIPpeakAnno")
readLines(example_pattern_file)

## [1] ">ExamplePattern" "GGNCCK"          ">ExamplePattern" "AACCNM"

pattern_in_peak <- summarizePatternInPeaks(patternFilePath = example_pattern_file,
                                           BSgenomeName = BSgenome.Hsapiens.UCSC.hg38,
                                           peaks = macs_peak_gr2[1:200],
                                           bgdForPerm = "chromosome",
                                           nperm = 100, 
                                           method = "permutation.test")

pattern_in_peak$motif_enrichment

##        patternNum totalNumPatternWithSameLen patternRate cutOffPermutationTest
## GGNCCK       1351                     169140 0.007987466           0.002140105
## AACCNM        586                     169140 0.003464586           0.003781629

head(pattern_in_peak$motif_occurrence, n = 2)

##   motifChr motifStartInChr motifEndInChr      motifName motifPattern
## 1     chr1           28746         28751 ExamplePattern       GGNCCK
## 2     chr1           29049         29054 ExamplePattern       GGNCCK
##   motifStartInPeak motifEndInPeak motifFound motifFoundStrand peakChr peakStart
## 1              406            411     GGACCG                +    chr1     28341
## 2              709            714     GGGCCG                +    chr1     28341
##   peakEnd peakWidth peakStrand
## 1   29610      1270          *
## 2   29610      1270          *

The summarziePatternInPeaks function provides two distinct methods (method = c("binom.test", "permutation.test")) for evaluating the significance of given motif enrichment within target peaks. When utilizing the “binom.test” method, the expected frequencies for each motif must be determined first, which can be achieved via expectFrequencyMethod = c("Markov", "Naive"). When employing the “permutation.test” method, users need to provide parameters such as the number of permutation (“nperm”), the significance level (“alpha”), and methods how background sequences are determined (bgdForPerm = c("shuffle", "chromosome"), chromosome = c("asPeak", "random")).

The outcome from summarizePatternInPeaks comprises two tables. The “motif_enrichment” table succinctly summarizes the enrichment statistics for each motif. Specifically, when method = "binom.test", the last column “pValueBinomTest” represents the P-values resulting from the binomial test. Conversely, when method = "permutation.test", the last column “cutOffPermutationTest” denotes the threshold value derived from the permutation test, determined by the specified signifiance level. The “motif_occurrence” table contains detailed information regarding the occurrences of each motif within every peak.

Alternative tools

Besides using ChIPpeakAnno, users can extract the sequences with the getAllPeakSequence function and employ other tools such as MEME Suite for motif-related analysis.

Use hypergeometric test to determine overlap among peak sets

The hypergeometric test is grounded in the principle of the hypergeometric distribution, which is a probability distribution that describes the likelihood of drawing a specific number of successes (k) from a sample size of n, given a finite population (N) that contains K successes when sampling is done without replacement. The null hypothesis posits that the sample is drawn randomly from the population, implying that the there is no overlap between the two sets of peaks. A small P-value indicates that the null should be rejected, indicating a significant overlap between the two sets of peaks.

In ChIPpeakAnno, the hypergeometric test is incorporated into the makeVennDiagram function, as detailed in Section @ref(venn). The following example illustrates how to compute the hypergeometric test P-values for peak sets from three TF ChIP-seq experiments.

tf1 <- toGRanges(system.file("extdata/TAF.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")
tf2 <- toGRanges(system.file("extdata/Tead4.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")
tf3 <- toGRanges(system.file("extdata/YY1.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")

To effectively apply the hypergeometric test, we first need to estimate the total number of binding sites, i.e. totalTest. The selection of totalTest affects the stringency of the test, with smaller values resulting in a more conservative outcome (larger P-values). For practical guidance on how to choose an appropriate value for totalTest, you can refer to this post.

In our example, we assume that potential binding regions (which include coding regions and promoter regions) constitute 3% (fairly rough estimation) of the entire genome. Since the example data is derived from chromosome 2, we can estimate the number of total binding sites as (length of chr2) * 3% / (mean peak length).

overlapping_peaks <- findOverlapsOfPeaks(tf1, 
                                         tf2, 
                                         tf3, 
                                         connectedPeaks = "keepAll")
mean_peak_width <- mean(width(unlist(GRangesList(overlapping_peaks[["all.peaks"]]))))

total_binding_sites <- length(BSgenome.Hsapiens.UCSC.hg38[["chr2"]]) * 0.03 / mean_peak_width
venn1 <- makeVennDiagram(overlapping_peaks, 
                         totalTest = total_binding_sites, 
                         connectedPeaks = "keepAll", 
                         fill = c("#CC79A7", "#56B4E9", "#F0E442"),
                         col = c("#D55E00", "#0072B2", "#E69F00"),
                         cat.col = c("#D55E00", "#0072B2", "#E69F00"))

Venn diagram1 showing overlapping peaks

For the P-values of each peak pair:

venn1[["p.value"]]

##      tf1 tf2 tf3         pval
## [1,]   0   1   1 4.565072e-52
## [2,]   1   0   1 0.000000e+00
## [3,]   1   1   0 2.744460e-88

For the overlapping peak counts:

venn1[["vennCounts"]]

##      tf1 tf2 tf3   Counts count.tf1 count.tf2 count.tf3
## [1,]   0   0   0 4953.594         0         0         0
## [2,]   0   0   1  621.000         0         0       621
## [3,]   0   1   0 2097.000         0      2097         0
## [4,]   0   1   1  309.000         0       310       319
## [5,]   1   0   0   59.000        59         0         0
## [6,]   1   0   1  166.000       172         0       170
## [7,]   1   1   0    8.000         8         8         0
## [8,]   1   1   1  476.000       545       537       521
## attr(,"class")
## [1] "VennCounts"

The parameter connectedPeaks denotes how to calculate peak counts when multiple peaks from different groups overlap. If connectedPeaks = "keepFirstListConsistent", the counts from the first group will be consistently displayed in the plot. If connectedPeaks = "keepAll", you will see multiple numbers: all original peak counts for each group will be displayed in parentheses, while the count of the minimally involved peaks will be displayed outside the parentheses.

venn2 <- makeVennDiagram(overlapping_peaks, 
                         totalTest = total_binding_sites, 
                         connectedPeaks = "keepFirstListConsistent", 
                         fill = c("#CC79A7", "#56B4E9", "#F0E442"),
                         col = c("#D55E00", "#0072B2", "#E69F00"),
                         cat.col = c("#D55E00", "#0072B2", "#E69F00"))

Venn diagram2 showing overlapping peaks

Given that all the P-values are extremely small, we must reject the null hypothesis. This indicates that each pair of peak sets significantly overlaps with each other. A major drawback of this approach is the necessity to estimate a totalTest, which could dramatically affect the test results. For instance, if we choose “2%” instead of “3%” in the above example, the P-value for tf1 vs. tf2 increases to 0.49, meaning we can no longer reject the null. To circumvent the requirement of totalTest, we have integrated a permutation test in the peakPermTest function.

Use permutation test to determine overlap among peak sets

The permutation test is a non-parametric test, which means it does not require the data to adhere to any specific distribution. The test statistic is determined according to the observed data, and the null distribution of the test statistic is estimated using a permutation (re-sampling) procedure.

In our scenario, the number of overlapping peaks is treated as the test statistic. Its null distribution is estimated by first re-sampling peaks from a random peak list (the peak pool, which represents all potential binding sites) followed by counting the number of overlapping peaks. The random peak list is generated using the distributions found from the input peaks, ensuring that the peak widths and relative binding positions to the features (such as TSS and geneEnd) follow the same distributions as the input peaks. When the null hypothesis is valid, the number of overlapping peaks is not significantly different from what would be expected by chance. Here are some sample codes using peakPermTest.

Example1 demonstrates a permutation test for non-relevant peak sets.

library(TxDb.Hsapiens.UCSC.hg38.knownGene)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
set.seed(123)

# Example1: non-relevant peak sets
utr5 <- unique(unlist(fiveUTRsByTranscript(txdb)))
utr3 <- unique(unlist(threeUTRsByTranscript(txdb)))

utr5 <- utr5[sample.int(length(utr5), 1000)]
utr3 <- utr3[sample.int(length(utr3), 1000)]

pt1 <- peakPermTest(peaks1 = utr3, 
                    peaks2 = utr5,
                    TxDb = txdb, 
                    maxgap = 500,
                    seed = 1)
plot(pt1)

Permutation test for non-relevant peak sets

Example2 demonstrates a permutation test for highly relevant peak sets.

# Example2: highly relevant peak sets
cds <- unique(unlist(cdsBy(txdb)))
ol <- findOverlaps(cds, utr3, maxgap = 1)
peaks2 <- c(cds[sample.int(length(cds), 500)],
            cds[queryHits(ol)][sample.int(length(ol), 500)])
pt2 <- peakPermTest(peaks1 = utr3,
                    peaks2 = peaks2,
                    TxDb = txdb,
                    maxgap = 500,
                    seed = 1)
plot(pt2)

Permutation test for highly relevant peak sets

As observed, for highly relevant peak sets, the P-value is very small.

Use custom peak pool

The peakPermTest function can automatically generate a peak pool given the peak binding type (bindingType = c("TSS", "geneEnd")), annotation type (featureType = c("transcript", "exon")), and annotation data (TxDb). Additionally, users have the option to construct the peak pool from actual binding sites derived from experimental data, with hot spots removed. Hot spots refer to genomic regions that have a high likelihood of being bound by many TFs in ChIP-seq experiments (Yip et al. 2012). We recommend removing hot spots prior to the permutation test to prevent overestimation of the association between two input peak sets. Users are also advised to remove ENCODE blacklist regions. The blacklists were constructed by identifying consistently problematic regions across independent cell lines and types of experiments for each species in the ENCODE and modENCODE datasets (Consortium 2012).

Below is an example of creating a peak pool for the human genome using the TF binding site clusters downloaded from ENCODE. The following steps are involved:

• Step1: download TF binding sites
• Step2: download hot spots
• Step3: liftover hot spots to hg38 (if necessary)
• Step4: select peak sets to test
• Step5: remove hot spots from binding pool
• Step6: perform permutation test

# Step1: download TF binding sites
temp <- tempfile()
download.file(file.path("https://hgdownload.cse.ucsc.edu/",
                        "goldenpath/",
                        "hg38/",
                        "encRegTfbsClustered/",
                        "encRegTfbsClusteredWithCells.hg38.bed.gz"), 
              temp)
df_tfbs <- read.delim(gzfile(temp, "r"), header = FALSE)
unlink(temp)

colnames(df_tfbs)[1:4] <- c("seqnames", "start", "end", "TF")
tfbs_hg38 <- GRanges(as.character(df_tfbs$seqnames),
                     IRanges(df_tfbs$start, df_tfbs$end),
                     TF = df_tfbs$TF)

# Step2: download hot spots
base_url <- "http://metatracks.encodenets.gersteinlab.org/metatracks/"

temp1 <- tempfile()
temp2 <- tempfile()
download.file(file.path(base_url, 
                        "HOT_All_merged.tar.gz"), 
              temp1)
download.file(file.path(base_url, 
                        "HOT_intergenic_All_merged.tar.gz"),
              temp2)
untar(temp1, exdir = dirname(temp1))
untar(temp2, exdir = dirname(temp1))
bedfiles <- dir(dirname(temp1), "bed$")
hot_spots_hg19 <- sapply(file.path(dirname(temp1), bedfiles), toGRanges, format = "BED")
unlink(temp1)
unlink(temp2)

names(hot_spots_hg19) <- gsub("_merged.bed", "", bedfiles)
hot_spots_hg19 <- sapply(hot_spots_hg19, unname)
hot_spots_hg19 <- GRangesList(hot_spots_hg19)

# Step3: liftover hot spots to hg38
library(R.utils)

temp_chain_gz <- tempfile()
temp_chain <- tempfile()
base_url_chain <- "http://hgdownload.cse.ucsc.edu/goldenpath/"
download.file(file.path(base_url_chain, 
                        "hg19/liftOver/",
                        "hg19ToHg38.over.chain.gz"),
              temp_chain_gz)
gunzip(filename = temp_chain_gz, destname = temp_chain)
chain_file <- import.chain(hg19_to_hg38)
unlink(temp_chain_gz)
unlink(temp_chain)

hot_spots_hg38 <- liftOver(hot_spots_hg19, chain_file)

# Step4: select peak sets to test
tfbs_hg38_by_TF <- split(tfbs_hg38, tfbs_hg38$TF)
TAF1 <- tfbs_hg38_by_TF[["TAF1"]]
TEAD4 <- tfbs_hg38_by_TF[["TEAD4"]]

# Step5: remove hot spots from binding pool
tfbs_hg38 <- subsetByOverlaps(tfbs_hg38, hot_spots_hg38, invert=TRUE)
tfbs_hg38 <- reduce(tfbs_hg38)

# Step6: perform permutation test
pool <- new("permPool", 
            grs = GRangesList(tfbs_hg38), 
            N = length(TAF1))
pt3 <- peakPermTest(TAF1, TEAD4, pool = pool, ntimes = 500)
plot(pt3)

Permutation test using custom peak pool

Visualize TSS enrichment signal of multiple experiments

If you have peak files obtained from multiple TF ChIP-seq experiments and would like to compare their enriched signals around TSS using raw signals such as read coverage. ChIPpeakAnno provides two functions, featureAlignedHeatmap and featureAlignedDistribution, to visualize their binding patterns side-by-side.

To illustrate this, we first need to prepare both the peak data and coverage information (bigWig). This involves four steps:

• Step1: read in example peak files
• Step2: find peaks that are shared by all
• Step3: read in example coverage files
• Step4: visualize binding patterns

# Step1: read in example peak files
extdata_path <- system.file("extdata", package = "ChIPpeakAnno")
broadPeaks <- dir(extdata_path, "broadPeak")
gr_TFs <- sapply(file.path(extdata_path, broadPeaks), toGRanges, format = "broadPeak")
names(gr_TFs) <- gsub(".broadPeak", "", broadPeaks)
names(gr_TFs)

## [1] "TAF"   "Tead4" "YY1"

Now, we have imported peak files for three TFs (“TAF”, “Tead”, and “YY1”) into the gr_TFs object. The next step is to identify overlapping peaks to ensure that we are comparing binding patterns over the same genomic regions.

# Step2: find peaks that are shared by all
ol <- findOverlapsOfPeaks(gr_TFs)
gr_TFs_ol <- ol$peaklist$`TAF///Tead4///YY1`

# Step3: read in example coverage files
# here we read in coverage data from -2000bp to -2000bp of each shared peak center
gr_TFs_ol_center <- reCenterPeaks(gr_TFs_ol, width = 4000) # use the center of the peaks and extend 2000bp upstream and downstream to obtain peaks with uniform length of 4000bp

bigWigs <- dir(extdata_path, "bigWig")
coverage_list <- sapply(file.path(extdata_path, bigWigs), 
                        import, # rtracklayer::import
                        format = "BigWig",
                        which = gr_TFs_ol_center,
                        as = "RleList")

names(coverage_list) <- gsub(".bigWig", "", bigWigs)
names(coverage_list)

## [1] "TAF"   "Tead4" "YY1"

Heatmap

# Step4: visualize binding patterns
sig <- featureAlignedSignal(coverage_list, gr_TFs_ol_center)

# since the bigWig files are only a subset of the original files,
# filter to keep peaks that are with coverage data for all peak sets
keep <- rowSums(sig[[1]]) > 0 & rowSums(sig[[2]]) > 0 & rowSums(sig[[3]]) > 0
sig <- sapply(sig, function(x) x[keep, ], simplify = FALSE)
gr_TFs_ol_center <- gr_TFs_ol_center[keep]

featureAlignedHeatmap(sig, gr_TFs_ol_center,
                      upper.extreme=c(3, 0.5, 4))

Heatmap of coverages for selected TFs

By default, the rows in the heatmap are ordered by the total coverage per row from the first sample (“TAF” in this example). We can reorder the rows by tuning the sortBy option. For example, setting sortBy = "YY1" will order the rows by the “YY1” sample in the dataset. You can also sort the rows based on hierarchical clustering results. Here is a demonstration:

# perform hierarchical clustering on rows
sig_rowsums <- sapply(sig, rowSums, na.rm = TRUE)
row_distance <- dist(sig_rowsums)
hc <- hclust(row_distance)

# use hierarchical clustering order to sort
gr_TFs_ol_center$sort_by <- hc$order
featureAlignedHeatmap(sig, gr_TFs_ol_center, 
                      upper.extreme = c(3, 0.5, 4),
                      sortBy = "sort_by")

Heatmap of coverages for selected TFs sorted by hierarchical clustering

Density plot

Additionally, we can create a density plot using the featureAlignedDistribution function.

featureAlignedDistribution(sig, gr_TFs_ol_center,
                           type = "l")

Distribution of coverage densities for selected TFs

Step1: import data

Common peak formats such as BED, GFF, and MACS can be converted into GRanges format using the toGRanges function. After importing the data, concordance across peak replicates will be evaluated with findOverlapsOfPeaks and makeVennDiagram. Be aware that the metadata columns will be dropped for the merged overlapping peaks. To add them back, we can use the addMetadata function. For example, addMetadata(ol, colNames = "score", FUN = mean) will add a “score” column to each merged overlapping peak by taking the mean score of individual peaks involved.

library(ChIPpeakAnno)

# Convert BED/GFF into GRanges
bed1 <- system.file("extdata", "MACS_output_hg38.bed", 
                    package = "ChIPpeakAnno")
gr1 <- toGRanges(bed1, format = "BED", header = FALSE)
gff1 <- system.file("extdata", "GFF_peaks_hg38.gff", 
                    package = "ChIPpeakAnno")
gr2 <- toGRanges(gff1, format = "GFF", header = FALSE)

# Find overlapping peaks
ol <- findOverlapsOfPeaks(gr1, gr2)

# Add "score" metadata column to overlapping peaks
ol <- addMetadata(ol, colNames = "score", FUN = mean) 
head(ol$mergedPeaks, n = 2)

## GRanges object with 2 ranges and 1 metadata column:
##       seqnames        ranges strand |                           peakNames
##          <Rle>     <IRanges>  <Rle> |                     <CharacterList>
##   [1]     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##   [2]     chr1 789471-791811      * |          gr2__003,gr1__MACS_peak_14
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

venn <- makeVennDiagram(ol,
                        fill = c("#009E73", "#F0E442"),
                        col = c("#D55E00", "#0072B2"),
                        cat.col = c("#D55E00", "#0072B2"))

Venn diagram showing the number of overlapping peaks for gr1 and gr2

For the P-value of hypergeometric test:

venn[["p.value"]]

##      gr1 gr2 pval
## [1,]   1   1    0

For overlapping peak counts:

venn[["vennCounts"]]

##      gr1 gr2 Counts count.gr1 count.gr2
## [1,]   0   0      0         0         0
## [2,]   0   1     61         0        61
## [3,]   1   0     63        63         0
## [4,]   1   1    165       167       168
## attr(,"class")
## [1] "VennCounts"

As observed, the extremely small P-value suggests a high relevance between the two peak sets, reflecting good consistency among experimental replicates.

Step2: prepare annotation file

Similar to the peak files, the annotation file must also be converted into a GRanges object. Annotation GRanges can be constructed from not only BED, GFF, and user-defined text files, but also from EnsDb and TxDb objects using the toGRanges function. For EnsDb and TxDb objects, annotation can also be prepared with accessor functions, as detailed in Section @ref(txdbensdb). Note that the version of genome used to create the annotation file must match the genome used for peak calling, as feature positions may vary across different genome releases. As an example, if you are using Mus_musculus.v103 for read mapping, it’s best to use EnsDb.Mmusculus.v103 for annotation.

library(EnsDb.Hsapiens.v86)

ensembl.hs86.gene <- toGRanges(EnsDb.Hsapiens.v86, feature = "gene")
head(ensembl.hs86.gene, n = 2)

## GRanges object with 2 ranges and 1 metadata column:
##                   seqnames      ranges strand |   gene_name
##                      <Rle>   <IRanges>  <Rle> | <character>
##   ENSG00000223972     chr1 11869-14409      + |     DDX11L1
##   ENSG00000227232     chr1 14404-29570      - |      WASH7P
##   -------
##   seqinfo: 357 sequences (1 circular) from 2 genomes (hg38, GRCh38)

Step3: visualize peak distribution

Peak distribution relative to features

Now, given the merged overlapping peaks and annotation data, we can visualize the distribution of the distance from the merged overlapping peaks to the nearest features, such as genes (TSSs), using the binOverFeature function.

binOverFeature(ol$mergedPeaks, 
               nbins = 20,
               annotationData = ensembl.hs86.gene,
               xlab = "peak distance from TSSs (bp)", 
               ylab = "peak count", 
               main = "Distribution of aggregated peak numbers around TSS")

Peak count distribution around transcription start sites

Peak distribution over multiple feature levels

We can use the genomicElementDistribution to summarize the distribution of peaks over different types of genomic features such exon, intron, enhancer, UTR. When inputting a single peak file, a pie graph will be generated.

library(TxDb.Hsapiens.UCSC.hg38.knownGene)

genomicElementDistribution(ol$mergedPeaks, 
                           TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)

Pie graph showing peak distributions over different genomic features

As can be seen, a significant number of peaks originate from promoter regions, which is consistent with the signature of peaks obtained from Tf binding experiments. When inputting a list of peak sets (e.g. replicates), a bar graph will be generated.

macs_peaks <- GRangesList(rep1 = gr1,
                          rep2 = gr2)
genomicElementDistribution(macs_peaks, 
                           TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)

Bar graph showing peak distributions over different genomic features for two replicates

According to the bar graph, the two peak replicates are consistent.

Peak overlappings for multiple features

To visualize peak overlap for multiple feature sets, we can utilize the UpSet plot (for details, see Section @ref(upset)).

library(UpSetR)

res <- genomicElementUpSetR(ol$mergedPeak,
                            TxDb.Hsapiens.UCSC.hg38.knownGene)
upset(res[["plotData"]], 
      nsets = length(colnames(res$plotData)), 
      nintersects = NA)

UpSet plot showing peak overlappings for multiple features

Step4: annotate peaks

Based on the above distribution of aggregated peak numbers around the TSS and the distribution of peaks in different chromosomal regions, most of the peaks locate near the TSS. Therefore, it is reasonable to adopt the feature-centric method to annotate the peaks residing in the promoter regions. Promoters can be specified with the bindingRegion parameter. In the following example, the promoter region is defined as 2000bp upstream and 500bp downstream of the TSS (bindingRegion = c(-2000, 500)).

ol_anno <- annotatePeakInBatch(ol$mergedPeak, 
                               AnnotationData = ensembl.hs86.gene, 
                               output = "nearestBiDirectionalPromoters",
                               bindingRegion = c(-2000, 500))
head(ol_anno, n = 2)

## GRanges object with 2 ranges and 9 metadata columns:
##      seqnames        ranges strand |                           peakNames
##         <Rle>     <IRanges>  <Rle> |                     <CharacterList>
##   X1     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##   X1     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##             peak         feature feature.ranges feature.strand  distance
##      <character>     <character>      <IRanges>          <Rle> <integer>
##   X1          X1 ENSG00000228327  725885-778626              -         0
##   X1          X1 ENSG00000237491  778770-810060              +         0
##      insideFeature distanceToSite     gene_name
##        <character>      <integer>   <character>
##   X1  overlapStart              0 RP11-206L10.2
##   X1  overlapStart              0 RP11-206L10.9
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

You can export the annotation to a CSV file:

ol_anno <- unname(ol_anno) # remove names to avoid duplicate row.names error
ol_anno$peakNames <- NULL # remove peakNames to avoid unimplemented type 'list' error
write.csv(as.data.frame(ol_anno), "ol_anno.csv")

To visualize the distribution of peaks around features:

pie1(table(ol_anno$insideFeature))

Pie chart showing peak distribution around features

peaks_near_BDP <- peaksNearBDP(ol$mergedPeaks, 
                    AnnotationData = ensembl.hs86.gene, 
                    MaxDistance = 5000) 
# MaxDistance will be translated into "bindingRegion = 
# c(-MaxDistance, MaxDistance)" internally

peaks_near_BDP$n.peaks

## [1] 165

peaks_near_BDP$n.peaksWithBDP

## [1] 18

peaks_near_BDP$percentPeaksWithBDP

## [1] 0.1090909

head(peaks_near_BDP$peaksWithBDP, n = 2)

## GRangesList object of length 2:
## $`1`
## GRanges object with 2 ranges and 10 metadata columns:
##      seqnames        ranges strand |                           peakNames
##         <Rle>     <IRanges>  <Rle> |                     <CharacterList>
##   X1     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##   X1     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##        bdp_idx        peak         feature feature.ranges feature.strand
##      <integer> <character>     <character>      <IRanges>          <Rle>
##   X1         1          X1 ENSG00000228327  725885-778626              -
##   X1         1          X1 ENSG00000237491  778770-810060              +
##       distance insideFeature distanceToSite     gene_name
##      <integer>   <character>      <integer>   <character>
##   X1         0  overlapStart              0 RP11-206L10.2
##   X1         0  overlapStart              0 RP11-206L10.9
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
## 
## $`4`
## GRanges object with 2 ranges and 10 metadata columns:
##      seqnames        ranges strand |                  peakNames   bdp_idx
##         <Rle>     <IRanges>  <Rle> |            <CharacterList> <integer>
##   X4     chr1 920981-921619      * | gr1__MACS_peak_17,gr2__005         4
##   X4     chr1 920981-921619      * | gr1__MACS_peak_17,gr2__005         4
##             peak         feature feature.ranges feature.strand  distance
##      <character>     <character>      <IRanges>          <Rle> <integer>
##   X4          X4 ENSG00000223764  916865-921016              -         0
##   X4          X4 ENSG00000187634  923928-944581              +      2308
##      insideFeature distanceToSite   gene_name
##        <character>      <integer> <character>
##   X4  overlapStart              0 RP11-54O7.3
##   X4      upstream           2308      SAMD11
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

library(EnsDb.Hsapiens.v75)

ensembl.hs75.gene <- toGRanges(EnsDb.Hsapiens.v75, feature = "gene")

DNA5C <- system.file("extdata", 
                     "wgEncodeUmassDekker5CGm12878PkV2.bed.gz",
                     package="ChIPpeakAnno")
DNAinteractiveData <- toGRanges(gzfile(DNA5C))
# the example bed.gz file can also be downloaded from:
# https://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/wgEncodeUmassDekker5C/wgEncodeUmassDekker5CGm12878PkV2.bed.gz

peaks_near_enhancer <-  findEnhancers(peaks = ol$mergedPeaks,
                                      annoData = ensembl.hs75.gene, 
                                      DNAinteractiveData = DNAinteractiveData)
head(peaks_near_enhancer, n = 2)

## GRanges object with 1 range and 14 metadata columns:
##      seqnames              ranges strand |                   peakNames
##         <Rle>           <IRanges>  <Rle> |             <CharacterList>
##   X1     chr1 151619224-151619324      * | gr2__229,gr1__MACS_peak_229
##              feature      feature.ranges feature.strand feature.shift.ranges
##          <character>           <IRanges>          <Rle>            <IRanges>
##   X1 ENSG00000143393 151264273-151300191              -  151619414-151655333
##      feature.shift.strand  distance insideFeature distanceToSite   gene_name
##                     <Rle> <integer>   <character>      <integer> <character>
##   X1                    +        89      upstream             89       PI4KB
##      DNAinteractive.idx        peak DNAinteractive.ranges DNAinteractive.blocks
##               <integer> <character>             <IRanges>         <IRangesList>
##   X1                814          X1   151283079-151636526  1-5699,335673-353448
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Step5: perform enrichment analysis

With annotated peaks, we can use the getEnrichedGO and getEnrichedPATH functions to respectively retrieve a list of enriched GO or pathway terms.

library(org.Hs.eg.db)

enriched_go <- getEnrichedGO(annotatedPeak = ol_anno, 
                             orgAnn = "org.Hs.eg.db", 
                             feature_id_type = "ensembl_gene_id",
                             condense = TRUE)
enrichmentPlot(enriched_go)

Histogram showing enriched GO terms

Likewise, the following pertains to pathway enrichment analysis.

library(reactome.db)

enriched_pathway <- getEnrichedPATH(annotatedPeak = ol_anno,
                                    orgAnn = "org.Hs.eg.db", 
                                    pathAnn = "KEGGREST",
                                    maxP = 0.05)
enrichmentPlot(enriched_pathway)

Histogram showing enriched pathways

# To remove the common suffix " - Homo sapiens (human)":
enrichmentPlot(enriched_pathway, label_substring_to_remove = " - Homo sapiens \\(human\\)")

Histogram showing enriched pathways

Step6: conduct motif analysis

To perform motif analysis, we first need to extract sequences surrounding the peaks, then acquire the consensus sequences. We can also visualize the top motifs discovered.

library(BSgenome.Hsapiens.UCSC.hg38)

seq_around_peak <- getAllPeakSequence(ol$mergedPeaks, 
                                      upstream = 20,
                                      downstream = 20, 
                                      genome = BSgenome.Hsapiens.UCSC.hg38)
head(seq_around_peak, n = 2)

## GRanges object with 2 ranges and 4 metadata columns:
##       seqnames        ranges strand |                           peakNames
##          <Rle>     <IRanges>  <Rle> |                     <CharacterList>
##   [1]     chr1 778411-780198      * | gr1__MACS_peak_13,gr2__001,gr2__002
##   [2]     chr1 789471-791811      * |          gr2__003,gr1__MACS_peak_14
##        upstream downstream               sequence
##       <numeric>  <numeric>            <character>
##   [1]        20         20 TGCGCCATGTTGCTAGGCTG..
##   [2]        20         20 AGAATTGAATTGAATGGACT..
##   -------
##   seqinfo: 1 sequence from an unspecified genome

We can use the write2FASTA function to store the result in a FASTA file. The following code snippet generates Z-scores for short oligos of length 6.

freqs <- oligoFrequency(BSgenome.Hsapiens.UCSC.hg38$chr1)
motif_summary <- oligoSummary(seq_around_peak, 
                              oligoLength = 6,
                              MarkovOrder = 3,
                              freqs = freqs,
                              quickMotif = TRUE)
zscore <- sort(motif_summary$zscore)
h <- hist(zscore, 
          breaks = 100, 
          main = "Histogram of Z-score")
text(x = zscore[length(zscore)], 
     y = h$counts[length(h$counts)] + 1, 
     labels = names(zscore[length(zscore)]), 
     adj = 0, 
     srt = 90)

Histogram showing Z-score distribution for oligos of length 6

We can use motifStack to visualize the top discovered motifs.

library(motifStack)

pfm <- new("pfm", mat = motif_summary$motifs[[1]],
           name = "sample motif 1")
motifStack(pfm)

Sequence log of detected motif 1

Step1: import data

tf1 <- toGRanges(system.file("extdata/TAF.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")
tf2 <- toGRanges(system.file("extdata/Tead4.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")
tf3 <- toGRanges(system.file("extdata/YY1.broadPeak", package = "ChIPpeakAnno"),
                 format = "broadPeak")

Step2: determine if there is significant overlap among peak sets

To examine the associations across different peak sets, ChIPpeakAnno implements both hypergeometric test (makeVennDiagram, for details see Section @ref(hypertest)) and permutation test (peakPermTest, for details see Section @ref(permtest)). For demonstration, here we use hypergeometric test.

library(BSgenome.Hsapiens.UCSC.hg38)

overlapping_peaks <- findOverlapsOfPeaks(tf1, 
                                         tf2, 
                                         tf3, 
                                         connectedPeaks = "keepAll")
mean_peak_width <- mean(width(unlist(GRangesList(overlapping_peaks[["all.peaks"]]))))

total_binding_sites <- length(BSgenome.Hsapiens.UCSC.hg38[["chr2"]]) * 0.03 / mean_peak_width
venn1 <- makeVennDiagram(overlapping_peaks, 
                         totalTest = total_binding_sites, 
                         connectedPeaks = "keepAll", 
                         fill = c("#CC79A7", "#56B4E9", "#F0E442"),
                         col = c("#D55E00", "#0072B2", "#E69F00"),
                         cat.col = c("#D55E00", "#0072B2", "#E69F00"))

Venn diagram showing overlapping peaks for three TFs

For the P-values of each peak pair:

venn1[["p.value"]]

##      tf1 tf2 tf3         pval
## [1,]   0   1   1 4.565072e-52
## [2,]   1   0   1 0.000000e+00
## [3,]   1   1   0 2.744460e-88

Given that the P-values are all extremely low, there is a significant overlap between each pair of peak sets.

Step3: visualize and compare the binding patterns

We can leverage heatmap and density plot to effectively visualize and compare the binding patterns of multiple TFs. For detailed instructions, refer to @ref(multiexp).

ol_tfs <- findOverlapsOfPeaks(tf1, tf2, tf3, 
                              connectedPeaks = "keepAll")
gr_ol_tfs <- ol_tfs$peaklist$`tf1///tf2///tf3`

TF_width <- width(gr_ol_tfs)
gr_ol_tfs_center <- reCenterPeaks(gr_ol_tfs, width = 4000)

extdata_path <- system.file("extdata", package = "ChIPpeakAnno")
bigWigs <- dir(extdata_path, "bigWig")
coverage_list <- sapply(file.path(extdata_path, bigWigs), 
                        import, # rtracklayer::import
                        format = "BigWig",
                        which = gr_ol_tfs_center,
                        as = "RleList")

names(coverage_list) <- gsub(".bigWig", "", bigWigs)

sig <- featureAlignedSignal(coverage_list, gr_ol_tfs_center)
# since the bigWig files are only a subset of the original files,
# filter to keep peaks that are with coverage data for all peak sets
keep <- rowSums(sig[[1]]) > 0 & rowSums(sig[[2]]) > 0 & rowSums(sig[[3]]) > 0
sig <- sapply(sig, function(x) x[keep, ], simplify = FALSE)
gr_ol_tfs_center <- gr_ol_tfs_center[keep]

featureAlignedHeatmap(sig, gr_ol_tfs_center,
                      upper.extreme=c(3, 0.5, 4))

Heatmap of coverages for selected TFs

By default, the rows are arranged according to the total coverage per row from the first sample. Below is an example of how to sort by the third sample (“YY1”) in the dataset through tuning the sortBy option.

featureAlignedHeatmap(sig, gr_ol_tfs_center,
                      upper.extreme=c(3, 0.5, 4),
                      sortBy = "YY1")

Heatmap of coverages for selected TFs sorted by YY1 sample

To generate a density plot, use the featureAlignedDistribution function.

featureAlignedDistribution(sig, gr_ol_tfs_center, 
                           type="l")

Distribution of coverage densities for selected TFs

As observed, the binding of “YY1” is significantly stronger, while the binding of “Tead” is considerably weaker.

How to import peaks?

ChIPpeakAnno provides a helper function, toGRanges, specifically designed for importing files. Users also have the option to use other tools, such as rtracklayer::import.

How to properly annotate peaks from ChIP-seq data?

ChIPpeakAnno offers a variety of options for annotating peaks, which can be specified with the output argument in the annotatePeakInBatch function.

• If you are interested in the nearest features surrounding the peaks, set output = "nearestLocation", output = "shortestDistance", or output = "overlapping". This is a peak-centric method, for details see Section @ref(peakcentric).
• If you are interested in the peaks binding to the promoter regions, you can set output = "upstream". For more flexibility, you can adjust the bindingRegion option. This is a feature-centric method, for details see Section @ref(featurecentric)

Should I use annotatePeakInBatch or annoPeaks?

You should use annotatePeakInBatch as it serves as the primary wrapper function that internally calls annoPeaks. Historically, annotatePeakInBatch was primarily composed of peak-centric methods. Over time, feature-centric methods, initially implemented in the annoPeaks function, were integrated into annotatePeakInBatch.

How to prepare annotation data for annotatePeakInBatch?

Ensure that the annotation data is in GRanges format for annotatePeakInBatch to work. You can convert TxDb, EnsDb, or even GFF and BED, etc. into GRanges using the toGRanges function. Additionally, you can leverage the getAnnotation function to fetch annotation data from BioMart. If needed, you can even create a personalized annotation file. For detailed instructions, see Section @ref(prepanno).

Can I output either the nearest or overlapping features?

This question pertains to this post.

When you set output = "both", both the “nearest” and “overlapping” features will be output. If your preference is to assign only one type of feature to each peak (either “nearest” or “overlapping”, with a preference for “overlapping” in cases where the “nearest” feature is not “overlapping”), you can use the following strategy: first, annotate peaks with overlapping features; second, annotate peaks lacking overlapping features with the nearest features; last, concatenate the two sets of results. The example codes are provided below:

library(ensembldb)
library(EnsDb.Hsapiens.v75)
data(myPeakList)
annoData <- annoGR(EnsDb.Hsapiens.v75)

# Step1: annotate peaks to the overlapping features, if "select = 'all'", multiple features can be assigned to a single peak.
anno_overlapping <- annotatePeakInBatch(myPeakList, 
                                        AnnotationData = annoData, 
                                        output = "overlapping", 
                                        select = "first")
anno_overlapping_non_na <- anno_overlapping[!is.na(anno_overlapping$feature)]

# Step2: annotate peaks that are without overlapping features to nearest features
myPeakList_non_overlapping <- myPeakList[!(names(myPeakList) %in% anno_overlapping_non_na$peak)]  
anno_nearest <- annotatePeakInBatch(myPeakList_non_overlapping,
                                    AnnotationData = annoData, 
                                    output = "nearestLocation", 
                                    select = "first")

# Step3: concatenate the two
anno_final <- c(anno_overlapping_non_na, anno_nearest)

The provided code allocates either the “overlapping” or “nearest” feature to each peak. If the “overlapping” feature is not the “nearest”, only the “overlapping” one will be reported.

Why is the sum of peak numbers in the Venn diagram NOT equal to the sum of the peaks in the original peak lists?

This question is a common scenario in calculating the intersection of peaks. Peaks, being a range of continuous points rather than single points, present a challenge in determining the overlapping regions. Consider the intersection of two lists of range objects: list A (1~3, 4~5, 7~9) and list B (2~8), the number of overlapping ranges depends on the reference chosen, if we use list B as the reference, the output would be one range; however, if we use list A as the reference, the output would be three ranges.

How does findOverlapsOfPeaks count the number of overlapping peaks?

This question is in response to this post.

When counting the number of overlapping peaks using findOverlapsOfPeaks, the connectedPeaks option comes into play. If multiple peaks involve in any group of connected or overlapping peaks in any input peak list, setting connectedPeaks = "merge" will increment the overlapping counts by one. On the other hand, setting connectedPeaks = "min" will add the minimal number of involved peaks in each group of connected or overlapped peaks to the overlapping counts. If connectedPeaks = "keepAll", it will add the number of involved peaks for each peak list to the corresponding overlapping counts. In addition, it will output counts as if connectedPeaks was set to “min”.

For examples, if 5 peaks in group1 overlap with 2 peaks in group 2:

• Setting connectedPeaks = "merge" will add 1 to the overlapping counts
• Setting connectedPeaks = "min" will add 2 to the overlapping counts
• Setting connectedPeaks = "keepAll" will add 5 peaks to “count.group1”, 2 to “count.group2”, and 2 to the overlapping counts

Is there a way to show the number of peaks in original peak lists?

You have the option to configure connectedPeaks = "keepAll" in both the findOverlapsOfPeaks and makeVennDiagram functions.

How to extract the original peak IDs of the overlapping peaks?

In the output of findOverlapsOfPeaks, each element in the peak list contains a metadata column names peakNames, which is a CharacterList. This CharacterList is a list of the contributing peak IDs with prefixes, e.g. “peaks1_peakname1”, “peaksi_peaknamej”, where “i” denotes the peak group i and “j” denotes the specific peak j within that group. Users can retrieve the original peak name by splitting these strings. Here are the sample codes:

library(ChIPpeakAnno)
library(reshape2)

# Step1: read in two peak files
bed <- system.file("extdata", 
                   "MACS_output.bed",
                   package = "ChIPpeakAnno")
gff <- system.file("extdata", 
                   "GFF_peaks.gff", 
                   package = "ChIPpeakAnno")
gr1 <- toGRanges(bed, format = "BED", 
                 header = FALSE)
gr2 <- toGRanges(gff, format = "GFF", 
                 header = FALSE, skip = 3)
names(gr2) <- seq(length(gr2)) # add names to gr2 for Step4

# Step2: find overlapping peaks
ol <- findOverlapsOfPeaks(gr1, gr2)
peakNames <- ol$peaklist[['gr1///gr2']]$peakNames

# Step3: extract original peak names
peakNames <- melt(peakNames, 
                  value.name = "merged_peak_id") # reshape df
head(peakNames, n = 2)

##   merged_peak_id NA                NA
## 1              1  1 gr1__MACS_peak_13
## 2              1  1            gr2__1

peakNames <- cbind(peakNames[, 1], 
                   do.call(rbind, 
                           strsplit(as.character(peakNames[, 3]), "__")))

colnames(peakNames) <- c("merged_peak_id", "group", "peakName")
head(peakNames, n = 2)

##      merged_peak_id group peakName      
## [1,] "1"            "gr1" "MACS_peak_13"
## [2,] "1"            "gr2" "1"

# Step4: split by peak group
gr1_subset <- gr1[peakNames[peakNames[, "group"] %in% "gr1", "peakName"]]
gr2_subset <- gr2[peakNames[peakNames[, "group"] %in% "gr2", "peakName"]]
head(gr1_subset, n = 2)

## GRanges object with 2 ranges and 1 metadata column:
##                seqnames        ranges strand |     score
##                   <Rle>     <IRanges>  <Rle> | <numeric>
##   MACS_peak_13     chr1 713791-715332      * |   2550.61
##   MACS_peak_14     chr1 724851-727191      * |     58.34
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

head(gr2_subset, n = 2)

## GRanges object with 2 ranges and 4 metadata columns:
##     seqnames        ranges strand |   source     type     score     phase
##        <Rle>     <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
##   1     chr1 713893-714747      + |  bed2gff region_0         0      <NA>
##   2     chr1 715023-715578      + |  bed2gff region_1         0      <NA>
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Alternatively, the following snippet should also work.

gr1_renamed <- ol$all.peaks$gr1
gr2_renamed <- ol$all.peaks$gr2
head(gr1_renamed, n = 2)

## GRanges object with 2 ranges and 1 metadata column:
##                    seqnames      ranges strand |     score
##                       <Rle>   <IRanges>  <Rle> | <numeric>
##   gr1__MACS_peak_1     chr1 28341-29610      * |    160.81
##   gr1__MACS_peak_2     chr1 90821-91234      * |    133.12
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

head(gr2_renamed, n = 2)

## GRanges object with 2 ranges and 4 metadata columns:
##          seqnames        ranges strand |   source     type     score     phase
##             <Rle>     <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
##   gr2__1     chr1 713893-714747      + |  bed2gff region_0         0      <NA>
##   gr2__2     chr1 715023-715578      + |  bed2gff region_1         0      <NA>
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

peakNames <- melt(ol$peaklist[['gr1///gr2']]$peakNames, 
                  value.name = "merged_peak_id")
gr1_subset <- gr1_renamed[peakNames[grepl("^gr1", peakNames[, 3]), 3]]
gr2_subset <- gr2_renamed[peakNames[grepl("^gr2", peakNames[, 3]), 3]]
head(gr1_subset, n = 2)

## GRanges object with 2 ranges and 1 metadata column:
##                     seqnames        ranges strand |     score
##                        <Rle>     <IRanges>  <Rle> | <numeric>
##   gr1__MACS_peak_13     chr1 713791-715332      * |   2550.61
##   gr1__MACS_peak_14     chr1 724851-727191      * |     58.34
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

head(gr2_subset, n = 2)

## GRanges object with 2 ranges and 4 metadata columns:
##          seqnames        ranges strand |   source     type     score     phase
##             <Rle>     <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
##   gr2__1     chr1 713893-714747      + |  bed2gff region_0         0      <NA>
##   gr2__2     chr1 715023-715578      + |  bed2gff region_1         0      <NA>
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

How to select the proper number for totalTest when using makeVennDiagram?

During the evaluation of the association between two sets of peak data using the hypergeometric distribution, it is essential to determine the total number of potential binding sites, the totalTest parameter in the makeVennDiagram. This parameter should have a value larger than the largest number of peaks in the peak list. The choice of totalTest influences teh stringency of the test, with smaller values yielding more stringent tests and larger P-values. The computation time for calculating P-values remains independent of the totalTest value. For practical guidance on selecting an appropriate totalTest value, please refer to this post.