EpiCompare Report
Sep-29-2022
knitr::opts_chunk$set(echo = FALSE)
### Create Output Directory ###
# if save_output = TRUE
if(params$save_output){
outfile_dir <- file.path(params$output_dir,"EpiCompare_file")
if(!dir.exists(outfile_dir)){
dir.create(outfile_dir, showWarnings = FALSE, recursive = TRUE)
}
}
#### ------ Prepare genome builds ------ ####
# e.g. genome_build <- list(reference="hg19",peakfiles="hg38",blacklist="hg19")
# or... genome_build <- "hg19"
builds <- prepare_genome_builds(genome_build = params$genome_build)
## Standardise all data to hg19 build
output_build <- prepare_output_build(params$genome_build_output)
#### ------ Prepare peaklist(s) ------ ####
# and check that the list is named, if not default filenames are used
peaklist <- prepare_peaklist(peaklist = params$peakfile)
peaklist <- liftover_grlist(grlist = peaklist,
input_build = builds$peaklist,
output_build = output_build)
#### ------ Prepare reference(s) ------ ####
reference <- prepare_reference(reference = params$reference)
reference <- liftover_grlist(grlist = reference,
input_build = builds$reference,
output_build = output_build)#### ------ Prepare blacklist ------ ####
blacklist <- liftover_grlist(grlist = params$blacklist,
input_build = builds$blacklist,
output_build = output_build)
### Standardise peaklist(s) ###
peaklist_tidy <- tidy_peakfile(peaklist = peaklist,
blacklist = blacklist)
### Standardise reference(s) ###
# and include in peaklist
reference_tidy <- reference
if (!is.null(reference)){
reference_tidy <- tidy_peakfile(peaklist = reference,
blacklist = blacklist)
peaklist_tidy <- c(peaklist_tidy, reference_tidy)
}
### Obtain Genome Annotation ###
txdb <- check_genome_build(genome_build = output_build)
### Dynamic Figure Height ###
fig_height <- fig_length(default_size = 7,
number_of_items = length(peaklist_tidy),
max_items = 10)EpiCompare compares epigenomic datasets for quality control and benchmarking purposes. The report consists of three sections:
- General Metrics: Metrics on peaks - percentage of blacklisted and non-standard peaks and peak widths - and fragments - duplication rate - of samples.
- Peak Overlap: Frequency, percentage, statistical significance of overlapping and non-overlapping peaks. This also includes Upset, precision-recall and correlation plots.
- Functional Annotation: Functional annotation (ChromHMM, ChIPseeker and enrichment analysis) of peaks. This also includes peak enrichment around TSS.
Input Datasets
- Reference peakfile: reference1
- Total of 2 peak files:
## [1] "File1: sample1"
## [1] "File2: reference1"Code
The function call used to create the report.
EpiCompare(peakfiles = list(sample1, reference1),
genome_build = hg19, hg38, hg19,
genome_build_output = hg19,
blacklist = blacklist,
picard_files = list(),
reference = reference1,
upset_plot = TRUE,
stat_plot = TRUE,
chromHMM_plot = TRUE,
chromHMM_annotation = "K562",
chipseeker_plot = TRUE,
enrichment_plot = TRUE,
tss_plot = TRUE,
precision_recall_plot = TRUE,
corr_plot=TRUE,
interact = TRUE,
save_output = TRUE,
output_dir = "/shared/bms20/projects/CUT_n_TAG/EpiCompare_pooled_peaks")1. General Metrics
Peak Information
Column Description:
- PeakN before tidy: Total number of peaks including those blacklisted and those in non-standard chromosomes.
- Blacklisted peaks removed (%): Percentage of blacklisted peaks in samples. E.g. ENCODE hg19 blacklist includes regions in the hg19 genome that have anomalous and/or unstructured signals independent of the cell-line or experiment.
- Non-standard peaks removed (%): Percentage of peaks in non-standard
and/or mitochondrial chromosomes. Identified using
BRGenomics::tidyChromosomes().
- PeakN after tidy: Total number of peaks after removing those in blacklisted regions and non-standard chromosomes.
NB: EpiCompare uses filtered peakfiles (i.e. datasets after removing peaks in blacklisted regions and non-standard chromosomes)
| PeakN Before Tidy | Blacklisted Peaks Removed (%) | Non-standard Peaks Removed (%) | PeakN After Tidy | |
|---|---|---|---|---|
| sample1 | 266274 | 0.463 | 0.342 | 264132 |
Fragment Information
Metrics on fragments is shown only if Picard summary is provided. See manual for help.
Column Description:
- Mapped_Fragments: Number of mapped read pairs in the file.
- Duplication_Rate: Percentage of mapped sequence that is marked as duplicate.
- Unique_Fragments: Number of mapped sequence that is not marked as duplicate.
Peak widths
Distribution of peak widths in samples
2. Peak Overlap
Percentage Overlap
Percentage of overlapping peaks between samples. Hover over heatmap for percentage values.
N.B. How to interpret heatmap: [Samples in x-axis of heatmap] peaks in [Samples in y-axis of heatmap] peaks
Statistical Significance
Depending on the format of the reference file, EpiCompare outputs different plots:
- Reference dataset has BED6+4 format (peak calling performed with MACS2): EpiCompare generates paired boxplot per sample showing the distribution of -log10(q-value) of reference peaks that are overlapping and non-overlapping with the sample dataset.
- Reference dataset does not have BED6+4 format: EpiCompare generates a barplot of percentage of overlapping sample peaks with the reference, coloured by statistical significance (adjusted p-value) of the overlap.
Reference peakfile: reference1
Keys:
- Overlap: Sample peaks in Reference peaks
- Unique: Sample peaks not in Reference peaks
Precision-Recall Curves
The first plot shows the balance between precision and recall across multiple peak calling stringency thresholds.
- Precision: Number of sample peaks in reference.
- Recall: Number of reference peaks in sample.
The second plot shows F1 score (a score that combines precision and recall) across the different peak calling stringency thresholds.
- F1:
2*(precision*recall) / (precision+recall)
Correlation Plot
The correlation plot shows the correlation between the quantiles when the genome is binned at a set size. These quantiles are based on the intensity of the peak, dependent on the peak caller used (q-value for MACS2):
3. Functional Annotation
3.1 ChromHMM
ChromHMM annotates and characterises peaks into different chromatin states. ChromHMM annotations used in EpiCompare were obtained from here.
- Cell-type annotation file used in this analysis: K562
All samples
ChromHMM annotation of individual samples.
Overlap: Sample peaks in Reference peaks
Percentage of Sample peaks found in Reference peaks (Reference peakfile: reference1)
| query | subject | overlap | total | Percentage | type | peaklist1 | peaklist2 |
|---|---|---|---|---|---|---|---|
| sample1 | reference1 | 25066 | 264132 | 9.49 | precision | sample1 | reference1 |
| reference1 | reference1 | 52544 | 52544 | 100.00 | precision | reference1 | reference1 |
| reference1 | sample1 | 51152 | 52544 | 97.40 | recall | sample1 | reference1 |
| reference1 | reference1 | 52544 | 52544 | 100.00 | recall | reference1 | reference1 |
ChromHMM annotation of sample peaks found in reference peaks.
Overlap: Reference peaks in Sample peaks
Percentage of Reference peaks found in Sample peaks (Reference peakfile: reference1)
| query | subject | overlap | total | Percentage | type | peaklist1 | peaklist2 |
|---|---|---|---|---|---|---|---|
| reference1 | sample1 | 51152 | 52544 | 97.40 | precision | reference1 | sample1 |
| reference1 | reference1 | 52544 | 52544 | 100.00 | precision | reference1 | reference1 |
| sample1 | reference1 | 25066 | 264132 | 9.49 | recall | reference1 | sample1 |
| reference1 | reference1 | 52544 | 52544 | 100.00 | recall | reference1 | reference1 |
ChromHMM annotation of reference peaks found in sample peaks.
Unique: Sample peaks not in Reference peaks
Percentage of sample peaks not found in reference peaks (Reference peakfile: reference1)
| query | subject | overlap | total | Percentage | type | peaklist1 | peaklist2 |
|---|---|---|---|---|---|---|---|
| sample1 | reference1 | 239066 | 264132 | 90.50 | precision | sample1 | reference1 |
| reference1 | reference1 | 0 | 52544 | 0.00 | precision | reference1 | reference1 |
| reference1 | sample1 | 1392 | 52544 | 2.65 | recall | sample1 | reference1 |
| reference1 | reference1 | 0 | 52544 | 0.00 | recall | reference1 | reference1 |
ChromHMM annotation of sample peaks not found in reference peaks.
Unique: Reference peaks not in Sample peaks
Percentage of reference peaks not found in sample peaks (Reference peakfile: reference1)
| query | subject | overlap | total | Percentage | type | peaklist1 | peaklist2 |
|---|---|---|---|---|---|---|---|
| reference1 | sample1 | 1392 | 52544 | 2.65 | precision | reference1 | sample1 |
| reference1 | reference1 | 0 | 52544 | 0.00 | precision | reference1 | reference1 |
| sample1 | reference1 | 239066 | 264132 | 90.50 | recall | reference1 | sample1 |
| reference1 | reference1 | 0 | 52544 | 0.00 | recall | reference1 | reference1 |
ChromHMM annotation of reference peaks not found in sample peaks.
3.2 ChIPseeker
EpiCompare uses ChIPseeker::annotatePeak() to annotate
peaks with the nearest gene and genomic regions where the peak is
located. The peaks are annotated with genes taken from human genome
annotations (hg19 or hg38) provided by Bioconductor.
3.3 Functional Enrichment Analysis
EpiCompare performs KEGG pathway and GO enrichment analysis using
clusterProfiler. ChIPseeker::annotatePeak() is
first used to assign peaks to nearest genes. Biological themes amongst
the genes are identified using ontologies (KEGG and GO). The peaks are
annotated with genes taken from annotations of human genome (hg19 or
hg38) provided by Bioconductor.
KEGG
GO
3.4 Peak Frequency around TSS
This plots peaks that are mapping to transcriptional start sites (TSS). TSS regions are defined as the flanking sequence of the TSS sites. The frequency of peaks in downstream (-3000bp) and upstream (+3000bp) of TSS is plotted. Faint color line around the main frequency line represents the 95% confidence interval estimated by bootstrap method.
## >> Running bootstrapping for tag matrix... 2022-09-29 17:37:55
## >> Running bootstrapping for tag matrix... 2022-09-29 17:38:28
4. Session Info
## R Under development (unstable) (2022-02-25 r81808)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rtracklayer_1.57.0 GenomicRanges_1.49.1 GenomeInfoDb_1.33.5 IRanges_2.31.2
## [5] S4Vectors_0.35.3 BiocGenerics_0.43.1 EpiCompare_1.1.1 PeakyFinders_0.99.2
## [9] stringr_1.4.1 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 R.methodsS3_1.8.2
## [3] ragg_1.2.2 tidyr_1.2.1
## [5] ggplot2_3.3.6 bit64_4.0.5
## [7] knitr_1.40 DelayedArray_0.23.1
## [9] R.utils_2.12.0 data.table_1.14.2
## [11] KEGGREST_1.37.3 RCurl_1.98-1.8
## [13] GEOquery_2.65.2 generics_0.1.3
## [15] GenomicFeatures_1.49.6 RSQLite_2.2.17
## [17] shadowtext_0.1.2 BSgenome.Hsapiens.UCSC.hg19_1.4.3
## [19] bit_4.0.4 tzdb_0.3.0
## [21] enrichplot_1.17.1 xml2_1.3.3
## [23] httpuv_1.6.6 SummarizedExperiment_1.27.2
## [25] assertthat_0.2.1 gargle_1.2.1
## [27] viridis_0.6.2 xfun_0.33
## [29] hms_1.1.2 jquerylib_0.1.4
## [31] evaluate_0.16 promises_1.2.0.1
## [33] fansi_1.0.3 restfulr_0.0.15
## [35] progress_1.2.2 caTools_1.18.2
## [37] dbplyr_2.2.1 readxl_1.4.1
## [39] igraph_1.3.5 DBI_1.1.3
## [41] geneplotter_1.75.0 htmlwidgets_1.5.4
## [43] googledrive_2.0.0 purrr_0.3.4
## [45] ellipsis_0.3.2 crosstalk_1.2.0
## [47] dplyr_1.0.10 annotate_1.75.0
## [49] gridBase_0.4-7 biomaRt_2.53.2
## [51] MatrixGenerics_1.9.1 vctrs_0.4.1
## [53] Biobase_2.57.1 remotes_2.4.2
## [55] cachem_1.0.6 withr_2.5.0
## [57] ggforce_0.3.4 BSgenome_1.65.2
## [59] genomation_1.29.0 vroom_1.5.7
## [61] GenomicAlignments_1.33.1 treeio_1.21.2
## [63] prettyunits_1.1.1 DOSE_3.23.2
## [65] ExperimentHub_2.5.0 ape_5.6-2
## [67] dir.expiry_1.5.0 lazyeval_0.2.2
## [69] crayon_1.5.1 basilisk.utils_1.9.2
## [71] genefilter_1.79.0 pkgconfig_2.0.3
## [73] labeling_0.4.2 tweenr_2.0.2
## [75] nlme_3.1-159 rlang_1.0.6
## [77] lifecycle_1.0.2 downloader_0.4
## [79] filelock_1.0.2 BiocFileCache_2.5.0
## [81] seqPattern_1.29.0 AnnotationHub_3.5.0
## [83] cellranger_1.1.0 rprojroot_2.0.3
## [85] polyclip_1.10-0 matrixStats_0.62.0
## [87] Matrix_1.4-1 aplot_0.1.7
## [89] boot_1.3-28 png_0.1-7
## [91] viridisLite_0.4.1 rjson_0.2.21
## [93] bitops_1.0-7 gson_0.0.8
## [95] R.oo_1.25.0 KernSmooth_2.23-20
## [97] Biostrings_2.65.3 blob_1.2.3
## [99] MACSr_0.99.15 qvalue_2.29.0
## [101] regioneR_1.29.1 readr_2.1.2
## [103] gridGraphics_0.5-1 echodata_0.99.14
## [105] scales_1.2.1 memoise_2.0.1
## [107] magrittr_2.0.3 plyr_1.8.7
## [109] gplots_3.1.3 zlibbioc_1.43.0
## [111] compiler_4.2.0 scatterpie_0.1.8
## [113] echoconda_0.99.7 BiocIO_1.7.1
## [115] RColorBrewer_1.1-3 plotrix_3.8-2
## [117] DESeq2_1.37.6 Rsamtools_2.13.4
## [119] cli_3.4.1 XVector_0.37.1
## [121] patchwork_1.1.2 MASS_7.3-58.1
## [123] tidyselect_1.1.2 stringi_1.7.8
## [125] textshaping_0.3.6 highr_0.9
## [127] yaml_2.3.5 GOSemSim_2.23.0
## [129] locfit_1.5-9.6 ggrepel_0.9.1
## [131] grid_4.2.0 sass_0.4.2
## [133] fastmatch_1.1-3 tools_4.2.0
## [135] parallel_4.2.0 rstudioapi_0.14
## [137] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 piggyback_0.1.4
## [139] gridExtra_2.3 BRGenomics_1.9.0
## [141] farver_2.1.1 ggraph_2.0.6
## [143] digest_0.6.29 BiocManager_1.30.18
## [145] shiny_1.7.2 Rcpp_1.0.9
## [147] BiocVersion_3.16.0 later_1.3.0
## [149] org.Hs.eg.db_3.15.0 httr_1.4.4
## [151] AnnotationDbi_1.59.1 colorspace_2.0-3
## [153] fs_1.5.2 XML_3.99-0.10
## [155] reticulate_1.26 splines_4.2.0
## [157] yulab.utils_0.0.5 tidytree_0.4.1
## [159] graphlayouts_0.8.1 basilisk_1.9.2
## [161] ggplotify_0.1.0 plotly_4.10.0
## [163] systemfonts_1.0.4 xtable_1.8-4
## [165] jsonlite_1.8.0 ggtree_3.5.3
## [167] tidygraph_1.2.2 UpSetR_1.4.0
## [169] ggfun_0.0.7 R6_2.5.1
## [171] pillar_1.8.1 htmltools_0.5.3
## [173] mime_0.12 glue_1.6.2
## [175] fastmap_1.1.0 clusterProfiler_4.5.2
## [177] DT_0.25 BiocParallel_1.31.12
## [179] interactiveDisplayBase_1.35.0 codetools_0.2-18
## [181] ChIPseeker_1.33.1 fgsea_1.23.0
## [183] utf8_1.2.2 lattice_0.20-45
## [185] bslib_0.4.0 tibble_3.1.8
## [187] curl_4.3.2 gtools_3.9.3
## [189] zip_2.2.1 GO.db_3.15.0
## [191] openxlsx_4.2.5 survival_3.4-0
## [193] limma_3.53.6 rmarkdown_2.16
## [195] munsell_0.5.0 DO.db_2.9
## [197] GenomeInfoDbData_1.2.8 impute_1.71.0
## [199] reshape2_1.4.4 gtable_0.3.1