HiCool
The 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.1
bowtie2 2.4.5
samtools 1.7
hicstuff 3.1.5
cooler 0.8.11
chromosight 1.6.3
The 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