HiCoolThe HiCool R/Bioconductor package provides an end-to-end interface to
process and normalize Hi-C paired-end fastq reads into .(m)cool files.
hicstuff python library
(https://github.com/koszullab/hicstuff).hicstuff.cooler (https://github.com/open2c/cooler)
library is used to parse pairs into a multi-resolution, balanced .mcool file.
.(m)cool is a compact, indexed HDF5 file format specifically tailored
for efficiently storing HiC-based data. The .(m)cool file format was
developed by Abdennur and Mirny and
published in 2019.basilisk environment.The main processing function offered in this package is HiCool().
To process .fastq reads into .pairs & .mcool files, one needs to provide:
r1 and r2);.fasta sequence file, a path to a pre-computed bowtie2 index
or a supported ID character (hg38, mm10, dm6, R64-1-1, WBcel235, GRCz10,
Galgal4);x <- HiCool(
r1 = '<PATH-TO-R1.fq.gz>',
r2 = '<PATH-TO-R2.fq.gz>',
restriction = '<RE1(,RE2)>',
resolutions = "<resolutions of interest>",
genome = '<GENOME_ID>'
)
Here is a concrete example of Hi-C data processing.
HiContactsData package..mcool file will have three levels of resolutions, from 1000bp to 8000bp.R64-1-1, the yeast genome reference.output/ directory.library(HiCool)
hcf <- HiCool(
r1 = HiContactsData::HiContactsData(sample = 'yeast_wt', format = 'fastq_R1'),
r2 = HiContactsData::HiContactsData(sample = 'yeast_wt', format = 'fastq_R2'),
restriction = 'DpnII,HinfI',
resolutions = c(4000, 8000, 16000),
genome = 'R64-1-1',
output = './HiCool/'
)
#> see ?HiContactsData and browseVignettes('HiContactsData') for documentation
#> loading from cache
#> see ?HiContactsData and browseVignettes('HiContactsData') for documentation
#> loading from cache
#> HiCool :: Recovering bowtie2 genome index from AWS iGenomes...
#> + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'create' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/HiCool/1.2.0/env' 'python=3.7.12' '--quiet' '-c' 'conda-forge' '-c' 'bioconda'
#> + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/HiCool/1.2.0/env' 'python=3.7.12'
#> + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/HiCool/1.2.0/env' '-c' 'conda-forge' '-c' 'bioconda' 'python=3.7.12' 'python=3.7.12' 'bowtie2=2.5.0' 'samtools=1.16.1' 'hicstuff=3.1.5' 'chromosight=1.6.3' 'cooler=0.9.1'
#> HiCool :: Initiating processing of fastq files [tmp folder: /tmp/Rtmp9rhi34/SWFYOC]...
#> HiCool :: Mapping fastq files...
#> HiCool :: Removing unwanted chromosomes...
#> HiCool :: Parsing pairs into .cool file...
#> HiCool :: Generating multi-resolution .mcool file...
#> HiCool :: Balancing .mcool file...
#> HiCool :: Tidying up everything for you...
#> HiCool :: .fastq to .mcool processing done!
#> HiCool :: Check ./HiCool/folder to find the generated files
#> HiCool :: Generating HiCool report. This might take a while.
#> HiCool :: Report generated and available @ /tmp/RtmpNTzYfK/Rbuild298c941c4975a5/HiCool/vignettes/HiCool/27bee7466dcb3a_7833^mapped-R64-1-1^SWFYOC.html
#> HiCool :: All processing successfully achieved. Congrats!
hcf
#> CoolFile object
#> .mcool file: ./HiCool//matrices/27bee7466dcb3a_7833^mapped-R64-1-1^SWFYOC.mcool
#> resolution: 4000
#> pairs file: ./HiCool//pairs/27bee7466dcb3a_7833^mapped-R64-1-1^SWFYOC.pairs
#> metadata(3): log args stats
S4Vectors::metadata(hcf)
#> $log
#> [1] "./HiCool//logs/27bee7466dcb3a_7833^mapped-R64-1-1^SWFYOC.log"
#>
#> $args
#> $args$r1
#> [1] "/home/biocbuild/.cache/R/ExperimentHub/27bee7466dcb3a_7833"
#>
#> $args$r2
#> [1] "/home/biocbuild/.cache/R/ExperimentHub/27bee77a85b142_7834"
#>
#> $args$genome
#> [1] "/tmp/Rtmp9rhi34/R64-1-1"
#>
#> $args$resolutions
#> [1] "4000"
#>
#> $args$resolutions
#> [1] "8000"
#>
#> $args$resolutions
#> [1] "16000"
#>
#> $args$restriction
#> [1] "DpnII,HinfI"
#>
#> $args$iterative
#> [1] TRUE
#>
#> $args$balancing_args
#> [1] " --min-nnz 10 --mad-max 5 "
#>
#> $args$threads
#> [1] 1
#>
#> $args$output
#> [1] "./HiCool/"
#>
#> $args$exclude_chr
#> [1] "Mito|chrM|MT"
#>
#> $args$keep_bam
#> [1] FALSE
#>
#> $args$scratch
#> [1] "/tmp/Rtmp9rhi34"
#>
#> $args$wd
#> [1] "/tmp/RtmpNTzYfK/Rbuild298c941c4975a5/HiCool/vignettes"
#>
#>
#> $stats
#> $stats$nFragments
#> [1] 1e+05
#>
#> $stats$nPairs
#> [1] 73993
#>
#> $stats$nDangling
#> [1] 10027
#>
#> $stats$nSelf
#> [1] 2205
#>
#> $stats$nDumped
#> [1] 83
#>
#> $stats$nFiltered
#> [1] 61678
#>
#> $stats$nDups
#> [1] 719
#>
#> $stats$nUnique
#> [1] 60959
#>
#> $stats$threshold_uncut
#> [1] 7
#>
#> $stats$threshold_self
#> [1] 7
Extra optional arguments can be passed to the hicstuff workhorse library:
iterative TRUE): By default, hicstuff first truncates your set of reads to 20bp and attempts to align the truncated reads, then moves on to aligning 40bp-truncated reads for those which could not be mapped, etc. This procedure is longer than a traditional mapping but allows for more pairs to be rescued. Set to FALSE if you want to perform standard alignment of fastq files without iterative alignment;balancing_args " --min-nnz 10 --mad-max 5 "): Specify here any balancing argument to be used by cooler when normalizing the binned contact matrices. Full list of options available at cooler documentation website;threads 1L): Number of CPUs to use to process data;exclude_chr 'Mito|chrM|MT'): List here any chromosome you wish to remove from the final contact matrix file;keep_bam FALSE): Set to TRUE if you wish to keep the pair of .bam files;scratch tempdir()): Points to a temporary directory to be used for processing.The important files generated by HiCool are the following:
<output_folder>/logs/<prefix>^mapped-<genome>^<hash>.log<output_folder>/matrices/<prefix>^mapped-<genome>^<hash>.mcool.pairs file: <output_folder>/pairs/<prefix>^mapped-<genome>^<hash>.pairs<output_folder>/plots/<prefix>^mapped-<genome>^<hash>_*.pdf.The diagnosis plots illustrate how pairs were filtered during the processing,
using a strategy described in Cournac et al., BMC Genomics 2012. The event_distance
chart represents the frequency of ++, +-, -+ and -- pairs in the library, as a function
of the number of restriction sites between each end of the pairs, and shows the inferred filtering threshold.
The event_distribution chart indicates the proportion of each type of pairs (e.g. dangling, uncut, abnormal, …)
and the total number of pairs retained (3D intra + 3D inter).
Notes:
.pairs file format is defined by the 4DN consortium;.(m)cool file format is defined by cooler authors in the supporting publication.Processing Hi-C sequencing libraries into .pairs and .mcool files requires
several dependencies, to (1) align reads to a reference genome, (2) manage
alignment files (SAM), (3) filter pairs, (4) bin them to a specific resolution
and (5)
All system dependencies are internally managed by basilisk. HiCool maintains
a basilisk environment containing:
python 3.9.1bowtie2 2.4.5samtools 1.7hicstuff 3.1.5cooler 0.8.11chromosight 1.6.3The first time HiCool() is executed, a fresh basilisk environment will
be created and required dependencies automatically installed. This ensures
compatibility between the different system dependencies needed to process
Hi-C fastq files.
sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] GenomicRanges_1.54.0 S4Vectors_0.40.0 IRanges_2.36.0
#> [4] HiContactsData_1.3.0 ExperimentHub_2.10.0 AnnotationHub_3.10.0
#> [7] BiocFileCache_2.10.0 dbplyr_2.3.4 BiocGenerics_0.48.0
#> [10] HiCool_1.2.0 HiCExperiment_1.2.0 BiocStyle_2.30.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.1.3 bitops_1.0-7
#> [3] rlang_1.1.1 magrittr_2.0.3
#> [5] matrixStats_1.0.0 compiler_4.3.1
#> [7] RSQLite_2.3.1 dir.expiry_1.10.0
#> [9] png_0.1-8 vctrs_0.6.4
#> [11] stringr_1.5.0 pkgconfig_2.0.3
#> [13] crayon_1.5.2 fastmap_1.1.1
#> [15] ellipsis_0.3.2 XVector_0.42.0
#> [17] rmdformats_1.0.4 utf8_1.2.4
#> [19] promises_1.2.1 rmarkdown_2.25
#> [21] sessioninfo_1.2.2 tzdb_0.4.0
#> [23] strawr_0.0.91 purrr_1.0.2
#> [25] bit_4.0.5 xfun_0.40
#> [27] zlibbioc_1.48.0 cachem_1.0.8
#> [29] GenomeInfoDb_1.38.0 jsonlite_1.8.7
#> [31] blob_1.2.4 later_1.3.1
#> [33] rhdf5filters_1.14.0 DelayedArray_0.28.0
#> [35] Rhdf5lib_1.24.0 BiocParallel_1.36.0
#> [37] interactiveDisplayBase_1.40.0 parallel_4.3.1
#> [39] R6_2.5.1 bslib_0.5.1
#> [41] stringi_1.7.12 reticulate_1.34.0
#> [43] jquerylib_0.1.4 Rcpp_1.0.11
#> [45] bookdown_0.36 SummarizedExperiment_1.32.0
#> [47] knitr_1.44 httpuv_1.6.12
#> [49] Matrix_1.6-1.1 tidyselect_1.2.0
#> [51] abind_1.4-5 yaml_2.3.7
#> [53] codetools_0.2-19 curl_5.1.0
#> [55] lattice_0.22-5 tibble_3.2.1
#> [57] withr_2.5.1 KEGGREST_1.42.0
#> [59] shiny_1.7.5.1 InteractionSet_1.30.0
#> [61] Biobase_2.62.0 basilisk.utils_1.14.0
#> [63] evaluate_0.22 Biostrings_2.70.0
#> [65] pillar_1.9.0 BiocManager_1.30.22
#> [67] filelock_1.0.2 MatrixGenerics_1.14.0
#> [69] stats4_4.3.1 plotly_4.10.3
#> [71] generics_0.1.3 vroom_1.6.4
#> [73] RCurl_1.98-1.12 BiocVersion_3.18.0
#> [75] ggplot2_3.4.4 munsell_0.5.0
#> [77] scales_1.2.1 xtable_1.8-4
#> [79] glue_1.6.2 lazyeval_0.2.2
#> [81] tools_4.3.1 BiocIO_1.12.0
#> [83] data.table_1.14.8 rhdf5_2.46.0
#> [85] grid_4.3.1 tidyr_1.3.0
#> [87] crosstalk_1.2.0 AnnotationDbi_1.64.0
#> [89] colorspace_2.1-0 GenomeInfoDbData_1.2.11
#> [91] basilisk_1.14.0 cli_3.6.1
#> [93] rappdirs_0.3.3 fansi_1.0.5
#> [95] S4Arrays_1.2.0 viridisLite_0.4.2
#> [97] dplyr_1.1.3 gtable_0.3.4
#> [99] sass_0.4.7 digest_0.6.33
#> [101] SparseArray_1.2.0 htmlwidgets_1.6.2
#> [103] memoise_2.0.1 htmltools_0.5.6.1
#> [105] lifecycle_1.0.3 httr_1.4.7
#> [107] mime_0.12 bit64_4.0.5