ChIPseeker is an R package for annotating ChIP-seq data analysis. It supports annotating ChIP peaks and provides functions to visualize ChIP peaks coverage over chromosomes and profiles of peaks binding to TSS regions. Comparison of ChIP peak profiles and annotation are also supported. Moreover, it supports evaluating significant overlap among ChIP-seq datasets. Currently, ChIPseeker contains 17,000 bed file information from GEO database. These datasets can be downloaded and compare with user’s own data to explore significant overlap datasets for inferring co-regulation or transcription factor complex for further investigation.
If you use ChIPseeker(Yu, Wang, and He 2015) in published research, please cite:
- Q Wang#, M Li#, T Wu, L Zhan, L Li, M Chen, W Xie, Z Xie, E Hu, S Xu, G Yu*. Exploring epigenomic datasets by ChIPseeker. Current Protocols, 2022, 2(10): e585.
- G Yu*, LG Wang, QY He*. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparision and visualization. Bioinformatics. 2015, 31(14):2382-2383.
Chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) has become standard technologies for genome wide identification of DNA-binding protein target sites. After read mappings and peak callings, the peak should be annotated to answer the biological questions. Annotation also create the possibility of integrating expression profile data to predict gene expression regulation. ChIPseeker(Yu, Wang, and He 2015) was developed for annotating nearest genes and genomic features to peaks.
ChIP peak data set comparison is also very important. We can use it as an index to estimate how well biological replications are. Even more important is applying to infer cooperative regulation. If two ChIP seq data, obtained by two different binding proteins, overlap significantly, these two proteins may form a complex or have interaction in regulation chromosome remodelling or gene expression. ChIPseeker(Yu, Wang, and He 2015) support statistical testing of significant overlap among ChIP seq data sets, and incorporate open access database GEO for users to compare their own dataset to those deposited in database. Protein interaction hypothesis can be generated by mining data deposited in database. Converting genome coordinations from one genome version to another is also supported, making this comparison available for different genome version and different species.
Several visualization functions are implemented to visualize the coverage of the ChIP seq data, peak annotation, average profile and heatmap of peaks binding to TSS region.
The datasets CBX6 and CBX7 in this vignettes were downloaded from GEO (GSE40740)(Pemberton et al. 2014) while ARmo_0M, ARmo_1nM and ARmo_100nM were downloaded from GEO (GSE48308)(Urbanucci et al. 2012) . ChIPseeker provides
readPeakFile to load the peak and store in
## $ARmo_0M ##  "/tmp/Rtmpg88Cs2/Rinstac2c760da5ec5/ChIPseeker/extdata/GEO_sample_data/GSM1174480_ARmo_0M_peaks.bed.gz" ## ## $ARmo_1nM ##  "/tmp/Rtmpg88Cs2/Rinstac2c760da5ec5/ChIPseeker/extdata/GEO_sample_data/GSM1174481_ARmo_1nM_peaks.bed.gz" ## ## $ARmo_100nM ##  "/tmp/Rtmpg88Cs2/Rinstac2c760da5ec5/ChIPseeker/extdata/GEO_sample_data/GSM1174482_ARmo_100nM_peaks.bed.gz" ## ## $CBX6_BF ##  "/tmp/Rtmpg88Cs2/Rinstac2c760da5ec5/ChIPseeker/extdata/GEO_sample_data/GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz" ## ## $CBX7_BF ##  "/tmp/Rtmpg88Cs2/Rinstac2c760da5ec5/ChIPseeker/extdata/GEO_sample_data/GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz"
## GRanges object with 1331 ranges and 2 metadata columns: ## seqnames ranges strand | V4 V5 ## <Rle> <IRanges> <Rle> | <character> <numeric> ##  chr1 815093-817883 * | MACS_peak_1 295.76 ##  chr1 1243288-1244338 * | MACS_peak_2 63.19 ##  chr1 2979977-2981228 * | MACS_peak_3 100.16 ##  chr1 3566182-3567876 * | MACS_peak_4 558.89 ##  chr1 3816546-3818111 * | MACS_peak_5 57.57 ## ... ... ... ... . ... ... ##  chrX 135244783-135245821 * | MACS_peak_1327 55.54 ##  chrX 139171964-139173506 * | MACS_peak_1328 270.19 ##  chrX 139583954-139586126 * | MACS_peak_1329 918.73 ##  chrX 139592002-139593238 * | MACS_peak_1330 210.88 ##  chrY 13845134-13845777 * | MACS_peak_1331 58.39 ## ------- ## seqinfo: 24 sequences from an unspecified genome; no seqlengths
ChIP peaks coverage plot
After peak calling, we would like to know the peak locations over the whole genome,
covplot function calculates the coverage of peak regions over chromosomes and generate a figure to visualize. GRangesList is also supported and can be used to compare coverage of multiple bed files.
Profile of ChIP peaks binding to TSS regions
First of all, for calculating the profile of ChIP peaks binding to TSS regions, we should prepare the TSS regions, which are defined as the flanking sequence of the TSS sites. Then align the peaks that are mapping to these regions, and generate the tagMatrix.
In the above code, you should notice that tagMatrix is not restricted to TSS regions. The regions can be other types that defined by the user. ChIPseeker expanded the scope of region. Users can input the
by parameters to get the regions they want.
Heatmap of ChIP binding to TSS regions
ChIPseeker provide a one step function to generate this figure from bed file. The following function will generate the same figure as above.
Users can use
nbin parameter to speed up.
Users can also use ggplot method to change the details of the figures.
Users can also profile genebody regions with
Sometimes there will be a need to explore the comparison of the peak heatmap over two regions, for example, the following picture is the peak over two gene sets. One possible scenery of using this method is to compare the peak heatmap over up-regulating genes and down-regulating genes. Here
txdb2 is the simulated gene sets obtain from
list object containing different regions information. The length of each part is correlated to the amount of regions.
txdb1 <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb2 <- unlist(fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene))[1:10000,] region_list <- list(geneX = txdb1, geneY = txdb2) peakHeatmap_multiple_Sets(peak = files[], upstream = 1000,downstream = 1000, by = c("geneX","geneY"), type = "start_site", TxDb = region_list,nbin = 800)
We also meet the need of ploting heatmap and peak profiling together.
Exploring several regions with heatmap and peak profiling is also supported.
txdb1 <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb2 <- unlist(fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene))[1:10000,] region_list <- list(geneX = txdb1, geneY = txdb2) peak_Profile_Heatmap(peak = files[], upstream = 1000, downstream = 1000, by = c("geneX","geneY"), type = "start_site", TxDb = region_list,nbin = 800)
Average Profile of ChIP peaks binding to TSS region
## >> plotting figure... 2023-04-25 16:23:55
plotAvgProf2 provide a one step from bed file to average profile plot. The following command will generate the same figure as shown above.
Confidence interval estimated by bootstrap method is also supported for characterizing ChIP binding profiles.
Profile of ChIP peaks binding to different regions
Referring to the issue #16 , we developed and improved several functions support start site region, end site region and body region of Gene/Transcript/Exon/Intron/3UTR/5UTR.
getBioRegion can prepare the different regions for ChIP peaks to bind.
getTagMatrix can accept
downstream parameters to get tagmatrix according to different needs.
plotPeakProf2 supports the plotting of profiles of peaks binding to different regions.Users can also create heatmap or average profile of ChIP peaks binding to these regions.
In order to plot body regions, a new methond
binning,was introduced to
getTagMatrix. The idea of
binning was derived from deeptools(Ramı́rez et al. 2016).
binning scaled the regions having different lengths to the equal length by deviding the regions into the same amounts of boxes. Because the amount of boxes is equal, the regions can be thought of scaling to equal length.
binning method can speed up the
getTagMatrix by changing the precision from bp to box(several bps).
There are three ways to plot these regions. First, users can use
getBioRegion to prepare the regions. Then align the peaks that are mapping to these regions, and generate the tagMatrix by
getTagMatrix. At Last, plot the figures by
plotPeakProf. Second, users can input
by parameters to
getTagMatrix to get the tagMatrix and plot the figures. Third, users can use
plotPeakProf2 to do everything in one step.
Binning method for profile of ChIP peaks binding to TSS regions
Here uses the method of inputting
type = "start_site" means the start site region.
by = "gene" means that it is the start site region of gene(TSS regions). If users want to use binning method, the
nbin method must be set.
Profile of ChIP peaks binding to body regions
We improved and developed several functions to plot body region of Gene/Transcript/Exon/Intron/3UTR/5UTR. If users want to get more information from the body region, we added
downstream parameters to functions in order to get flank extension of body regions.
downstream can be NULL(default), rel object and actual numbers. NULL(default) reflects body regions with no flank extension. Rel object reflects the percentage of total length of body regions. Actual numbers reflects the actual length of flank extension.
## Here uses `plotPeakProf2` to do all things in one step. ## Gene body regions having lengths smaller than nbin will be filtered ## A message will be given to warning users about that. ## >> 9 peaks(0.872093%), having lengths smaller than 800bp, are filtered... ## the ignore_strand is FALSE in default. We put here to emphasize that. ## We will not show it again in the below example plotPeakProf2(peak = peak, upstream = rel(0.2), downstream = rel(0.2), conf = 0.95, by = "gene", type = "body", nbin = 800, TxDb = txdb, weightCol = "V5",ignore_strand = F)
Users can also get the profile ChIP peaks binding to gene body regions with no flank extension or flank extension decided by actual length.
## The first method using getBioRegion(), getTagMatrix() and plotPeakProf() to plot in three steps. genebody <- getBioRegion(TxDb = txdb, by = "gene", type = "body") matrix_no_flankextension <- getTagMatrix(peak,windows = genebody, nbin = 800) plotPeakProf(matrix_no_flankextension,conf = 0.95) ## The second method of using getTagMatrix() and plotPeakProf() to plot in two steps matrix_actual_extension <- getTagMatrix(peak,windows = genebody, nbin = 800, upstream = 1000,downstream = 1000) plotPeakProf(matrix_actual_extension,conf = 0.95)
Users can also get the body region of 5UTR/3UTR.
## >> loading peak file... 2023-04-25 16:23:55 ## >> preparing features information... 2023-04-25 16:23:55 ## >> identifying nearest features... 2023-04-25 16:23:55 ## >> calculating distance from peak to TSS... 2023-04-25 16:23:56 ## >> assigning genomic annotation... 2023-04-25 16:23:56 ## >> adding gene annotation... 2023-04-25 16:24:09 ## >> assigning chromosome lengths 2023-04-25 16:24:09 ## >> done... 2023-04-25 16:24:09
Note that it would also be possible to use Ensembl-based
EnsDb annotation databases created by the ensembldb package for the peak annotations by providing it with the
TxDb parameter. Since UCSC-style chromosome names are used we have to change the style of the chromosome names from Ensembl to UCSC in the example below.
Peak Annotation is performed by
annotatePeak. User can define TSS (transcription start site) region, by default TSS is defined from -3kb to +3kb. The output of
csAnno instance. ChIPseeker provides
as.GRanges to convert
GRanges instance, and
as.data.frame to convert
data.frame which can be exported to file by
TxDb object contained transcript-related features of a particular genome. Bioconductor provides several package that containing
TxDb object of model organisms with multiple commonly used genome version, for instance TxDb.Hsapiens.UCSC.hg38.knownGene, TxDb.Hsapiens.UCSC.hg19.knownGene for human genome hg38 and hg19, TxDb.Mmusculus.UCSC.mm10.knownGene and TxDb.Mmusculus.UCSC.mm9.knownGene for mouse genome mm10 and mm9, etc. User can also prepare their own
TxDb object by retrieving information from UCSC Genome Bioinformatics and BioMart data resources by R function
TxDb object should be passed for peak annotation.
All the peak information contained in peakfile will be retained in the output of
annotatePeak. The position and strand information of nearest genes are reported. The distance from peak to the TSS of its nearest gene is also reported. The genomic region of the peak is reported in annotation column. Since some annotation may overlap, ChIPseeker adopted the following priority in genomic annotation.
- 5’ UTR
- 3’ UTR
Downstream is defined as the downstream of gene end.
ChIPseeker also provides parameter genomicAnnotationPriority for user to prioritize this hierachy.
annotatePeak report detail information when the annotation is Exon or Intron, for instance “Exon (uc002sbe.3/9736, exon 69 of 80)”, means that the peak is overlap with an Exon of transcript uc002sbe.3, and the corresponding Entrez gene ID is 9736 (Transcripts that belong to the same gene ID may differ in splice events), and this overlaped exon is the 69th exon of the 80 exons that this transcript uc002sbe.3 prossess.
Parameter annoDb is optional, if provided, extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added. The geneId column in annotation output will be consistent with the geneID in TxDb. If it is ENTREZID, ENSEMBL will be added if annoDb is provided, while if it is ENSEMBL ID, ENTREZID will be added.
Visualize Genomic Annotation
To annotate the location of a given peak in terms of genomic features,
annotatePeak assigns peaks to genomic annotation in “annotation” column of the output, which includes whether a peak is in the TSS, Exon, 5’ UTR, 3’ UTR, Intronic or Intergenic. Many researchers are very interesting in these annotations. TSS region can be defined by user and
annotatePeak output in details of which exon/intron of which genes as illustrated in previous section.
Pie and Bar plot are supported to visualize the genomic annotation.
Since some annotation overlap, user may interested to view the full annotation with their overlap, which can be partially resolved by
We extend UpSetR to view full annotation overlap. User can user
We can combine
upsetplot by setting vennpie = TRUE.
Visualize distribution of TF-binding loci relative to TSS
The distance from the peak (binding site) to the TSS of the nearest gene is calculated by
annotatePeak and reported in the output. We provide
plotDistToTSS to calculate the percentage of binding sites upstream and downstream from the TSS of the nearest genes, and visualize the distribution.
Functional enrichment analysis
Once we have obtained the annotated nearest genes, we can perform functional enrichment analysis to identify predominant biological themes among these genes by incorporating biological knowledge provided by biological ontologies. For instance, Gene Ontology (GO)(Ashburner et al. 2000) annotates genes to biological processes, molecular functions, and cellular components in a directed acyclic graph structure, Kyoto Encyclopedia of Genes and Genomes (KEGG)(Kanehisa et al. 2004) annotates genes to pathways, Disease Ontology (DO)(Schriml et al. 2011) annotates genes with human disease association, and Reactome(Croft et al. 2013) annotates gene to pathways and reactions.
ChIPseeker also provides a function, seq2gene, for linking genomc regions to genes in a many-to-many mapping. It consider host gene (exon/intron), promoter region and flanking gene from intergenic region that may under control via cis-regulation. This function is designed to link both coding and non-coding genomic regions to coding genes and facilitate functional analysis.
Enrichment analysis is a widely used approach to identify biological themes. I have developed several Bioconductor packages for investigating whether the number of selected genes associated with a particular biological term is larger than expected, including DOSE(Yu et al. 2015) for Disease Ontology, ReactomePA for reactome pathway, clusterProfiler(Yu et al. 2012) for Gene Ontology and KEGG enrichment analysis.
## ID Description ## R-HSA-186712 R-HSA-186712 Regulation of beta-cell development ## R-HSA-452723 R-HSA-452723 Transcriptional regulation of pluripotent stem cells ## GeneRatio BgRatio pvalue p.adjust qvalue ## R-HSA-186712 12/476 41/10955 9.857442e-08 9.936301e-05 9.936301e-05 ## R-HSA-452723 8/476 24/10955 4.760174e-06 2.399128e-03 2.399128e-03 ## geneID ## R-HSA-186712 2494/5080/3651/3175/6928/390874/3642/4821/4825/2255/222546/389692 ## R-HSA-452723 27022/2034/57167/27086/6657/63978/2516/2649 ## Count ## R-HSA-186712 12 ## R-HSA-452723 8
## ID Description GeneRatio ## R-HSA-186712 R-HSA-186712 Regulation of beta-cell development 10/386 ## R-HSA-383280 R-HSA-383280 Nuclear Receptor transcription pathway 10/386 ## BgRatio pvalue p.adjust qvalue ## R-HSA-186712 41/10955 1.107266e-06 0.001038615 0.001012857 ## R-HSA-383280 52/10955 1.105103e-05 0.004110375 0.004008434 ## geneID Count ## R-HSA-186712 2494/3651/4821/4825/2255/222546/389692/5080/6928/3642 10 ## R-HSA-383280 2494/5241/2100/9611/4306/7025/7101/2516/2649/9971 10
More information can be found in the vignettes of Bioconductor packages DOSE(Yu et al. 2015), ReactomePA, clusterProfiler(Yu et al. 2012), which also provide several methods to visualize enrichment results. The clusterProfiler(Yu et al. 2012) is designed for comparing and visualizing functional profiles among gene clusters, and can directly applied to compare biological themes at GO, DO, KEGG, Reactome perspective.
ChIP peak data set comparison
Profile of several ChIP peak data binding to TSS region
plotPeakProf can accept a list of
tagMatrix and visualize profile or heatmap among several ChIP experiments, while
plotPeakProf2 can accept a list of bed files and perform the same task in one step.
## >> plotting figure... 2023-04-25 16:24:50
## normal method plotPeakProf2(files, upstream = 3000, downstream = 3000, conf = 0.95, by = "gene", type = "start_site", TxDb = txdb, facet = "row") ## binning method plotPeakProf2(files, upstream = 3000, downstream = 3000, conf = 0.95, by = "gene", type = "start_site", TxDb = txdb, facet = "row", nbin = 800)
Profile of several ChIP peak data binding to body region
plotPeakProf2 also support to plot profile of several ChIP peak data binding to body region.