Cooltools: Enabling high-resolution Hi-C analysis in Python

1This is a PLOS Computational Biology Software paper.

2Chromosome conformation capture technologies are powerful tools for experimentally characterizing higher-order genome organization. Multiple labs and consortia, including the 4D Nucleome project (4DN), the International Nucleome Consortium (INC), and ENCODE, are generating increasingly larger datasets to probe genome organization at ever higher resolutions. Larger datasets increase the challenges at each step of computational analysis, from storage to memory, to researchers’ time—at the highest resolutions, the data for even a single chromosome does not fit easily in memory.

3Computational analyses of genome folding typically employ multiple approaches to quantify the variety of patterns in contact maps. These include: (i) enrichment of contacts within chromosomes, referred to as chromosome territories, (ii) decreasing contact frequency versus distance within chromosomes, or P(s); (ii) plaid patterns of preferential contacts, referred to as compartmentalization, (iii) squares of enriched contact frequency, referred to as topologically associated domains (TADs), (iv) focal peaks of pairwise contact frequencies, referred to as loops or dots. To quantitatively describe these features, and assay how they change across conditions, maps are summarized in many ways. Summaries include 2D pairwise averages, 1D profiles, and lists of genomic positions. The variety of Hi-C features and summaries benefits from a set of modular tools, rather than end-to-end pipelines, for data analysis.

4The features in genome folding maps vary substantially through cell states, across organisms, and between protocols. As C-data is obtained from a growing variety of contexts, computational analyses must co-evolve and expand as well. Because of this, computational researchers require APIs that can not only perform standard analyses but also enable the design of new analyses. This includes flexibility over parameters and thresholds, but also encompasses more dramatic extensions.

5Python provides one of the most popular environments for flexible data analysis and software development. This is in part due to the rich ecosystem of numpy, scipy, scikit-learn, and pandas. For Hi-C analysis in Python, the widely-adopted cooler format and library present an ideal foundation. Cooler handles storage and access for high-resolution datasets via a sparse data model. To fully leverage the advantages of cooler, however, requires a corresponding analysis framework in Python.

6While there are existing computational suites for extracting features of genome folding via the CLI or an R API, only FAN-C and TADbit have Python APIs (Table 1). However, FAN-C and TADbit neither take full advantage of the Python ecosystem nor provide consistent low-level access to the data, limiting their flexibility. We note the TADbit framework was developed before cooler was available, and hence without sparse storage considerations, and had an orthogonal primary focus on building 3D models.

7Here we introduce cooltools, a suite of computational methods that enable the analysis of high-resolution genome folding data accessible using cooler (Fig 1). Cooltools features a Python API paired with command line tools. The methods in cooltools are modular and provide access to both high-level and low-level APIs. The high-level API enables users to reproducibly perform routine analysis. The low-level API allows for troubleshooting challenges posed by new organisms or cell states and developing new advanced analyses. cooltools is part of a research community and broader ecosystem of open-source chromosome data analysis tools, Open2C (https://open2c.github.io/). cooltools provides detailed introductions to key concepts in Hi-C-data analysis with interactive notebook documentation. Collectively, cooltools enables computational analysis and methods development to keep pace with the rapidly-increasing quantity of experimental data and rapidly-evolving understanding of genome organization.

8The design of cooltools stems from its goal to enable flexible analyses of large quantities of high-resolution Hi-C data.

9Provide a modular set of tools with a balance between comprehensiveness and flexibility. Modular tools can be used separately or combined in a custom order. The set of core tools extracts and quantifies the majority of common features in the Hi-C literature.

10Separate processing from data analysis. Cooltools works with processed Hi-C data and is not concerned with mapping workflows.

11Separate data analysis from plotting, with a focus on the former. For Hi-C data, plotting often has a greater range of user-specific needs than feature extraction and quantification.

12Provide a paired Python API and CLI. CLI tools enable their use in computational pipelines and have a one-to-one correspondence with user-facing Python API functions. The user-facing API and CLI both have a standard format for arguments and take a cooler as their first input. For other types of genomic data besides Hi-C, the cooltools API accepts pandas dataframes, while the CLI reuses existing formats when possible: BED-like files for interval data and bedgraphs or bigwigs for genomic tracks.

13Split the Python API into two layers: a high-level API for user-facing functions and a flexible low-level APIs. A modular, consistent, and well-documented internal Python API allows experienced users to customize existing analysis routines or create entirely new ones.

14Enable analyses for custom genomes and non-model organisms. This includes user-defined subsets of regions from a reference genome, e.g. chromosomal arms. Leveraging its use of standard Python data structures, cooltools enables these analyses with bioframe. A view specifies a unique 1D coordinate axis with an ordered set of non-overlapping genomic intervals, also referred to as regions.

15Leverage indexed sparse storage. Cooltools is built directly on top of the cooler storage format and library, which allows it to operate on sparse matrices and/or out-of-core, either on raw counts or normalized contact matrices.

16Allow parallel processing to take advantage of modern multi-core CPU architectures. Most operations are performed via iteration over genomic regions or chunks of non-zero pixels, enabling straightforward multiprocessing.

17Leverage the Python data science stack to promote interoperability. Cooltools delegates to numpy, scipy, and pandas when possible and aims to follow the style of these APIs. To accelerate specific low-level computations, cooltools uses the just-in-time compiler numba.

18Return widely-used Python data structures, pandas dataframes, and numpy arrays, for user-facing outputs, and avoid the use of custom classes. This allows cooltools to operate on low-level and easily modifiable Python objects with an intuitive, consistent, and flexible interface.

19Cooltools is open-source and structured to encourage future contributions. For experimental code, cooltools provides a sandbox to enable easy sharing, safe testing, and eventual adoption of novel analysis routines.

20The main cooltools modules were initially developed with mammalian interphase Hi-C maps in mind, but the flexible implementation of cooltools has made it useful for analysis of data from different organisms (e.g., Drosophila, yeast, nematodes, chicken) and cellular states (mitosis, meiosis). We explain cooltools in tripartite chapters for each module below, describing the relevant features, the implementation, and analyses where the flexible API has proven useful. Conceptually, cooltools modules can be grouped in two main ways. First, modules can be grouped by whether they perform feature extraction (compartments, insulation, dots) versus feature quantification (expected, saddles, pileups). Alternatively, modules can be grouped by whether they make use of chromosome-wide information (expected, compartments, saddles) or more local information (insulation, dots, pileups). The design choice to separate analysis from plotting facilitated the creation of interactive notebook documentation (https://cooltools.readthedocs.io/en/latest/notebooks). These notebooks also provide flexible starting points for custom visualization and analysis. To highlight the practical application of cooltools, the figures throughout this manuscript are all generated using the same HFF micro-C dataset from Krietenstein et al., stored as an.mcool file, which provides a convenient access to sparse matrices binned at multiple resolutions, and the code to generate all figures is additionally available for download as interactive notebooks (https://github.com/open2c/open2c_vignettes/tree/main/cooltools_manuscript).

21Contact frequency decreases rapidly with genomic separation within chromosomes, which arises because chromosomes are long polymers This feature is interchangeably referred to as the: (i) expected because locus-specific features increase or decrease contact frequency relative to this major trend, (ii) scaling which is borrowed from the polymer physics literature, (iii) P(s), or contact probability, P, as a function of genomic separation, s. The shape of P(s) reflects folding patterns at different length scales, from those of adjacent nucleosomes to full chromosomes, and can be quantified with its derivative. P(s) varies through the cell cycle, across cell types, and is modified by structural maintenance of chromosomes complexes (SMCs) in interphase, mitosis, and meiosis (review).

22To quantify P(s) at sub-kilobase resolution (Fig 2), cooltools avoids conversion from sparse-to-dense matrices. Instead, the calculation is done by iterating over chunks of pixels stored in coolers, assigning pixels to regions, and aggregating at each distance. Avoiding conversion to dense also requires keeping track of filtered-out bins and a non-trivial accounting of the number of intersecting filtered bin pairs per-region per-distance. By allowing users to specify a correction weight, cooltools also enables P(s) calculation on raw data. To reduce variation at large genomic separations due to increasing sparsity of the data, cooltools returns P(s) curves smoothed in log-space using a numba-accelerated approach. This also allows cooltools to more robustly compute the derivative of P(s) and observed-over-expected maps. To quantify P(s) for arbitrary sets of regions, users can specify a view. The resulting curves for each region are then optionally averaged together. To enable interoperability with downstream analyses that require accounting for P(s), cooltools defines an expected format with specifications and checks. This includes both pairs of regions within chromosomes (cis expected) or pairs of regions on different chromosomes (trans expected).

23The flexible approach for computing P(s) used in cooltools has proven useful in multiple situations. For example, the ability to specify regions with a view has been used to calculate both within-arm P(s) and within-compartment P(s). Within-arm P(s) was particularly useful in mitosis, as P(s) is the main point of comparison with models of large-scale chromosome organization and is altered greatly across centromeric regions. Within-compartment P(s) revealed that extrusion occurs in compartment B as well as in A.

24A second prominent signal in mammalian interphase contact maps is plaid patterns that appear both within and between chromosomes, termed chromosome compartments. These plaid patterns reflect tendencies of chromosome regions to make more frequent contacts with regions of the same type: active regions have increased contact frequency with other active regions, and inactive regions tend to contact other inactive regions more frequently. The strength of compartmentalization varies through the cell cycle, across cell types, and after the degradation of components of the cohesin complex (for review).

25Cooltools calculates compartmentalization profiles using matrix eigenvector decomposition (Fig 3). For cis profiles, cooltools first adjusts the contact matrix for distance dependence before decomposition. For trans profiles, cooltools implements the ICE approach which masks cis portions of the map with randomly sampled trans data. Since eigenvectors are determined up to a sign and compartmentalization patterns are sometimes better captured by second or third eigenvectors, calculated vectors often need to be oriented and sorted according to their correlation with known markers of activity (e.g., GC content, density of genes). Cooltools allows users to request an arbitrary number of eigenvectors and can sort and orient eigenvectors according to a provided “phasing track”, most commonly the GC content per genomic bin.

26The flexible implementation makes it useful for downstream analysis. For example, modifies cooltools to extract 50+ eigenvectors, which are then used to find multiple classes of compartments. In Drosophila, GC content did not always reliably orient eigenvectors, but the flexible interface allowed gene density to be used instead.

27Contact frequency preferences between pairs of genomic loci can depend on multiple properties, including chromatin state and position along a chromosomal arm. Pairwise summaries, obtained by averaging over pairs of regions assigned to a set of classes, are a straightforward way to quantify these preferences. Pairwise summaries are commonly visualized as a 2D heatmap. For genome folding, these are most frequently used to quantify the strength of compartmentalization after assigning classes from the extracted compartment profile. The pattern of preferences revealed by this averaging (where contacts are disfavored between active and inactive chromatin) led to these summaries being referred to as “saddle plots”. Pairwise summaries have also been used to quantify GC content and fragment length biases.

28To quantify these preferences with cooltools saddle (Fig 4), the approach is to: (i) group bins with similar properties into classes, (ii) average over all map pixels specified by the bins in each pair of classes (i.e. for each pair of classes p, q, summarize the collection of pixels {Aij | i in p, j in q}). Bins can be assigned to classes either by discrete properties, or by discretizing a quantitative signal. cooltools operates on observed-over-expected contact frequencies, either based on raw counts or balanced data, and can be restricted to either cis or trans portions of the genome wide map.

29The flexibility of cooltools saddle enables 2D summaries to be used in many ways. This includes considering the average pairwise preferences over a reduced set of four histone modification states. Saddle plots have also been used to quantify contact preferences relative to the nuclear lamina and speckles, as well as differences between experimental protocols. Finally, since extraction of compartment profiles is decoupled from quantification, cooltools can be used to quantify how compartmental preferences shift with respect to a reference map, either through the cell cycle or in single cells.

30At small genomic distances (~2 Mb and shorter), mammalian interphase contact maps display square-like features along the diagonal. These 2D features, also referred to as domains, or TADs, are often summarized with 1D profiles of insulation. Insulation profiles quantify the frequency of contacts across each genomic bin averaged over a sliding diamond-shaped window. Sites with lower scores, reflecting reduced contact frequencies between upstream and downstream loci, are often referred to as insulating boundaries and coincide with the edges of domains. Insulation scores are popular for analysis due to their simplicity and relative robustness to Hi-C sequencing depth. Obtaining the positions of boundaries along the genome has in turn helped uncover relationships between genome folding and key genome-organizing proteins like CTCF and cohesin.

31To enable calculation of insulation profiles for high resolution data (Fig 5), cooltools iterates and aggregates over the sparse cooler pixel table. The size of the sliding window is chosen by the user. cooltools detects insulating boundaries as local minima of the insulation profile and quantifies their strength as their topographic prominence. Random variations in insulation along the genome lead to many local minima with very small prominences that are not desired in downstream analysis. To select strong boundaries from this set of minima, cooltools relies on Li and Otsu automated thresholding criteria borrowed from the image processing field (from scikit-image).

32Cooltools can calculate insulation profiles and boundaries at extremely high resolution (100 and 200 bp resolution). The flexible window size has also been used to analyze Micro-C boundaries at multiple levels. Finally, the flexibility of cooltools enabled its use to calculate insulation even with assays modified to detect sister chromatids.

33Punctate pairwise peaks of contact frequency are another prevalent feature of mammalian interphase contact maps. These features are also referred to as ‘dots’ or ‘loops’ in the literature, and can appear either in isolation or as parts of grids and at the corners of domains. In mammalian interphase, peaks seen between loci at distances <2Mb are often enriched for binding of CTCF and cohesin. In the context of loop extrusion, such peaks can emerge by the reproducible positioning of extruded cohesin loops at CTCF barriers. Peaks in mammalian interphase maps can also emerge due to other complexes, like polycomb. Peaks have also been observed and quantified in other contexts, including yeast meiosis and mitosis.

34Cooltools implements a HiCCUPs-style approach for peak detection (Fig 6): briefly, (i) maps are queried with a dense set of tiles that fit easily into memory, (ii) convolutional kernels are swept across map tiles to score all pixels for local enrichment (iii) significant pixels are determined by thresholding relative to their local background expected values, using a modified Benjamini-Hochberg FDR procedure (iv) adjacent significant pixels are optionally clustered and further filtered by enrichment. Cooltools enables control over dot calling at multiple levels, including: the area of the map queried by tiling, allowing user-specified kernels, choosing the FDR, and whether to perform the final clustering and filtering.

35The cooltools dots module has been used in multiple settings, including: calling peaks in MicroC data, calling peaks across conditions, and in yeast meiosis. In yeast meiosis, dot-calling benefited from the flexibility in cooltools to choose a smaller kernel and relaxed filters in the clustering step.

36Different classes of genomic elements display distinct local contact patterns depending on their propensity to form features like boundaries or dots. These patterns are typically quantified by extracting local 2D contact maps, or snippets, and then averaging snippets to create pileups or average Hi-C maps. Averaging can reveal genome-wide trends that might otherwise be masked by locus-specific structures or noise in individual snippets. Due to the symmetry of Hi-C maps, pileups come in two varieties: (i) on-diagonal pileups around a set of anchor positions, and (ii) off-diagonal pileups around a set of paired-anchors. On-diagonal pileups are useful for quantifying the strength of boundaries, which define a set of anchors, and off-diagonal pileups are useful for quantifying the strength of dots, which define a set of paired-anchors. Both types of pileup can be useful for determining the relationship between other genomic data (e.g., CTCF binding) with genome folding.

37To calculate pileups (Fig 7), cooltools first retrieves snippets for the specified set of anchors, a procedure termed snipping. Pileups can be built either directly from observed contact map snippets or from snippets normalized by a specified expected table. This is especially useful when the pileups are created for regions with different P(s), and specific expected tables are used for different regions. Users can also specify their choice of balancing weights or whether to calculate pileups directly on raw data. Multiprocessing can be specified to speed up snipping for large sets of anchors.

38The modular cooltools API facilities interactive analysis of snippets. The flexibility of using raw data for cooltools snipping has enabled its use with single-cell Hi-C. The flexibility offered by cooltools enables more complex analyses, including considering aggregate-peak-analysis (APA) as a function of genomic distance between anchors. In, minimal modifications to the cooltools code enabled nucleotide-resolution pileups, revealing patterns of nucleosomes around dots, CTCF anchors and gene starts. Leveraging the flexibility of cooltools, coolpuppy enables distance-stratified analyses as well as additional functionality for pileups. The separation of feature extraction and quantification has also enabled the quantification of features like TAD strength and peak strength in more sparsely sequenced datasets relative to a more deeply sequenced reference (e.g., for CTCF mutants).

39In addition to the main modules described above, cooltools provides additional tools for computing coverage, adaptive smoothing, and resampling.

40An early observation from Hi-C maps was that contacts are more frequent within chromosomes than between them, even at large genomic separations. This corresponds to the established concept of chromosome territoriality (for review). In cooltools, cis and total coverage profiles are calculated by iterating over chunks of pixels followed by summation. One usage of this feature was to determine that dot anchors in yeast meiotic Hi-C maps had higher proportions of cis contacts.

41Since Hi-C maps are often sparse, especially at high resolutions, understanding and mitigating the impact of sampling noise is often a key component of computational analysis. To aid such analyses, cooltools provides a module for generating randomly downsampled coolers. This allowed controlled comparisons across stages of Drosophila germline cell differentiation. To mitigate sampling noise, cooltools provides an adaptive smoothing method that iteratively splits observed counts across multiple bins until reaching a user-defined threshold, similar to the approach in Serpentine. Adaptive smoothing has proven useful for applications ranging from visualization to preprocessing data for neural network training.

42In summary, cooltools provides a suite of tools for genome folding analysis that are seamlessly integrated with cooler. Available on GitHub since 2017, cooltools modernize, extend, and replace the analysis routines in hiclib, the first open-source package for processing and analyzing Hi-C data. Currently, cooltools and FAN-C are the only available Hi-C analysis suites that provide a Python API and operate on the.cool data format (Table 1). Benchmarked on Hi-C data, cooltools shows better performance, especially at higher resolutions (S1 Fig).

43The methods in cooltools are modular, enabling the development of new analyses. The Python API in cooltools yields interpretable and structured outputs (Pandas dataframes or NumPy arrays), facilitating both data exploration as well as future tool development. For example, since the output of cooltools insulation is a pandas DataFrame of genomic intervals with additional score columns, the set of boundaries set can be interactively intersected with other genomic data (e.g., CTCF occupancy or motifs) using bioframe.

44The flexibility and breadth of the library make cooltools useful across a broad range of contexts. Cooltools have been used to analyze Hi-C data across cell states, including mitosis and meiosis, and for organisms ranging from bacteria to yeast to invertebrates to mammals. Analysis of data for non-standard genomes is facilitated by the use of genomic views, which support flexible ordering and naming of chromosomes, as implemented in bioframe. Indeed, the flexibility of cooltools made it the tool-of-choice for the largest evaluation of experimental Hi-C protocols to date.

45Cooltools provides a foundation for multiple future directions. First, integration with dask could be used to further improve the performance of snippers and other tools. Second, cooltools could be extended to support emerging types of contact data, including asymmetric contacts (e.g., RNA-vs-DNA), and single-cell data. Finally, cooltools can be used to develop reproducible analyses for Hi-C datasets. This includes reports, either on single datasets or for cross-dataset analysis, workflow routines for reproducible analysis (e.g., https://github.com/open2c/quaich), and integration with visualization platforms (e.g., HiGlass).