Pairtools: From sequencing data to chromosome contacts

1Chromosome conformation capture technologies (3C+), particularly Hi-C, revolutionized the study of genome folding by using high-throughput sequencing to measure spatial proximity. In all 3C+ approaches, spatial proximity is captured via DNA ligation, resulting in libraries of chimeric DNA molecules. Contact frequency between pairs of genomic loci is then read out by sequencing these libraries. All 3C+ protocols involve four key steps: (i) chemical cross-linking of chromatin, (ii) partial digestion of DNA followed by (iii) DNA ligation, and (iv) sequencing. Ligation is the pivotal step that records the spatial proximity of DNA loci in chimeric DNA sequences. In a typical 3C-based experiment, the DNA ligation products are subsequently shredded into smaller pieces, purified, and amplified. The resulting libraries are typically sequenced in short-read paired-end mode (around 50–300 bp on each side), producing millions to tens of billions of sequencing reads.

2Processing 3C+ sequencing data requires highly-performant and specialized computational tools. First, the quantity of 3C+ data is increasing rapidly. A growing number of labs and consortia (4DN, ENCODE, DANIO-CODE DCC ) use Hi-C to produce large quantities of proximity ligation data. At the same time, new protocols, such as Micro-C and Hi-C 3.0, improve resolution and sensitivity and generate even larger datasets. This requires software to be fast, parallelizable, and storage-efficient. Emerging 3C+ technologies introduce additional challenges, including measuring contacts within individual homologs and sister chromatids, single cells, and multi-way contacts (Pore-C and MC-3C ). This large variation in the growing family of 3C+ methods thus requires software to be versatile and flexible.

33C+ data are typically computationally processed in three stages (Figure 1a). First, sequencing reads are aligned against the reference genome. Next, pairs of genomic locations are extracted from the alignments. These pairs may be interpreted as genomic contact events. For various statistical and technical reasons, pairs are normally aggregated or binned to form contact matrices at various lower genomic resolutions. Binned contact maps can be stored, manipulated, and analyzed using downstream tools and software packages, such as cooler and cooltools from Open2C. However, a great deal of pertinent experimental information in 3C+ data accessible at the level of pairs is lost during the process of binning. This includes classification and quality control information, such as exact mapping positions, strand orientation, mapping quality, and pair type. As a result, pairs serve as a resource for generating custom filtered contact matrices. For these reasons, it is important to be able to flexibly generate, interpret, store, and manipulate pairs-level data.

4Here we introduce Pairtools, a suite of flexible and performant tools for converting sequenced 3C+ libraries into captured chromosome contacts. Pairtools provides modules to (i) extract and classify pairs from sequences of chimeric DNA molecules, (ii) deduplicate, filter, and manipulate resulting contact lists, and (iii) generate summary statistics for quality control (QC). Pairtools enables the processing of standard Hi-C as well as many Hi-C-derived protocols, including homolog- and sister-sensitive, multi-contact, single-cell Hi-C protocols. Pairtools uses a well-defined data format and may be arranged in highly efficient data pipes (Figure 1a, b). The benchmark against several popular 3C+ data mappers shows advantages of pairtools for high-performance and flexible 3C+ analysis. Pairtools is implemented in Python, powered by common data analysis libraries such as NumPy and pandas, offers a CLI, and is available as open-source software at: https://github.com/open2c/pairtools/.

5Pairtools provides tools for each step of data processing between the sequence alignment and the contacts binning (Figure 1a, b): extraction, processing, filtering, and quality control of contact pairs.

6Pairtools adheres to the following design principles, aligned with Unix style principles: Split functionality into tools that can be used independently or combined into pipelines.Focus on modularity, flexibility, and clarity first and performance second. Data processing should be as fast as alignment, but not necessarily faster.Outsource functionality when possible. Rely on existing software for alignment and workflow managers for data pipelining.Leverage the rich ecosystem of Python and data analysis libraries, including NumPy, pandas, scipy and scikit-learn.Use a standardized tabular format for pairs.Accommodate existing Hi-C protocol modifications by generalizing existing tools. When not possible, introduce protocol-specific tools.Take advantage of multi-processing and streaming to improve performance.

7Split functionality into tools that can be used independently or combined into pipelines.

8Focus on modularity, flexibility, and clarity first and performance second. Data processing should be as fast as alignment, but not necessarily faster.

9Outsource functionality when possible. Rely on existing software for alignment and workflow managers for data pipelining.

10Leverage the rich ecosystem of Python and data analysis libraries, including NumPy, pandas, scipy and scikit-learn.

11Use a standardized tabular format for pairs.

12Accommodate existing Hi-C protocol modifications by generalizing existing tools. When not possible, introduce protocol-specific tools.

13Take advantage of multi-processing and streaming to improve performance.

14Pairtools processes 3C+ data in three essential steps (Fig 1b). First, the genomic alignments of 3C+ products are parsed into individual contact events, or pairs. Second, the resulting pairs are sorted to facilitate data access and analysis. Third, pairs are deduplicated, resulting in the final list of 3C+ contacts.

15The minimal pairtools-based pipeline is expressed concisely as:

16Below, we describe these three steps and the corresponding pairtools functionality in detail.

17The DNA molecules in a 3C+ library are chimeric by design: spatial proximity between different genomic segments is captured as DNA ligation events, which are then read out via DNA sequencing. Pairtools makes use of existing software for alignment, taking.sam/.bam files as input. Sequence alignments in.sam/.bam comprehensively describe the structure of chimeric DNA molecules in the 3C+ library. Each entry in these files stores an alignment of a continuous read segment to the reference genome. Entries include mapping position, as well as flags and tags describing mapping properties (such as uniqueness of mapping, nucleotide variations, error probability, and more). The properties can be read with tools like pysam. However, alignments in the.sam/.bam files are reported sequentially and are not structured as contact pairs. Extraction of proximity ligation events from alignments requires additional processing, which pairtools achieves with parse.

18Pairtools parse is developed and optimized for Hi-C libraries with chimeric DNA molecules formed via single ligation events. Pairtools relies on local sequence aligners (e.g., bwa mem) that can split chimeric reads into segments, each aligned to a different location in the genome. Pairtools parse extracts pairs of adjacent alignments and reports them as individual contacts, categorizing them based on uniqueness (Fig 2a, Supplementary Figure 1a–c) and exposing the properties of alignments relevant to Hi-C. Parse also detects cases when one of the DNA fragments in a pair was sequenced on both sides of the read and thus is not a contact but an artifact of paired-end sequencing. Parse detects such cases, and merges the alignments, “rescuing” the true contact pair (Supplementary Figure 1d). The output of parse adheres to the standard format.pairs, discussed below.

19The engine of pairtools parse uses pysam to extract tags and flags from.sam/.bam files. Pairtools parse can run in combination with a variety of popular local sequence aligners, such as bwa mem, bwa mem2, minimap2, and others, as long as their output complies with the.sam/.bam format.

20Early Hi-C protocols aimed to produce libraries of chimeric molecules containing single ligation sites. However, newer experiments routinely produce molecules containing multiple ligation events, also known as multi-way contacts (Figure 2b, Supplementary Figure 1d–g), for the following reasons: The latest 3C+ protocols use more frequently cutting restriction enzymes (and their combinations ) or non-specific endonucleases, resulting in smaller ligated DNA fragments and more reads from multiple ligation events.The technologies for paired-end “short read” Illumina sequencing can now generate reads of several hundred nucleotides per side, which often exceeds the average fragment length.More dramatically, new 3C+ protocols use long-read sequencing technology to purposefully generate molecules with multiple ligation sites to detect co-occurring proximity events.

21The latest 3C+ protocols use more frequently cutting restriction enzymes (and their combinations ) or non-specific endonucleases, resulting in smaller ligated DNA fragments and more reads from multiple ligation events.

22The technologies for paired-end “short read” Illumina sequencing can now generate reads of several hundred nucleotides per side, which often exceeds the average fragment length.

23More dramatically, new 3C+ protocols use long-read sequencing technology to purposefully generate molecules with multiple ligation sites to detect co-occurring proximity events.

24The increasing prevalence of datasets with multi-way contacts calls for a tool dedicated to their extraction. This process differs substantially from the extraction of single proximity ligation events.

25First, the structure of multi-way contact is much more complex than a single proximity event (Figure 2b). Each library DNA molecule is expected to contain multiple fragments connected by pairwise direct ligation junctions. Multi-way contacts are ordered by the position within the chimeric DNA molecule. Since the order of the ligation events along a 3C+ chimeric molecule can contain meaningful information, multi-way contacts with a defined order are referred to as walks (Figure 2b). For example, this information is used to quantify either chromosomal entanglements or traversals between A and B compartments with walks.

26Further complexity arises when the contacts are categorized by the number of ligation events. When two fragments are sequenced one after another on the same read, they form direct ligation because the ligation junction between them was observed (Figure 2b). They are also called single hops because there is a single disruption of continuous DNA fragment. Some types of 3C+ analyses augment the data to include indirect ligations (mediated by multiple ligation events, two-, three- and many-hops) by combinatorial combinations of fragments in the same walk.

27Long chimeric DNA molecules can be sequenced by either traditional paired-end or single-end long-read sequencing. Paired-end sequencing introduces additional complications and uncertainty into analysis.

28In paired-end sequencing, if the 3C+ library DNA molecules are much longer than two read lengths, then the central side of the DNA molecule can remain unsequenced (Figure 2b, Supplementary Figure 1e–f). This insert can potentially contain other DNA fragments, resulting in indirect ligation. Thus, contacts formed by the ends of paired reads should be interpreted as unobserved ligations.

29The opposite scenario in paired-end sequencing is when the length of the DNA molecule is smaller than two read lengths, resulting in a read-through. The same subregion of a molecule is sequenced twice, producing internal duplicates (Supplementary Figure 1g).

30Pairtools parse2 extracts multi-way contacts from single 3C+ chimeric molecules resolving these issues (Figure 2a) and provides rich and flexible reporting (Supplementary Figure 1h). Parse2 works with both paired “short” reads and long read technologies and has the following features: Parse2 preserves and reports the order of ligation events (a walk) within a molecule (Figure 2a).Parse2 removes internal duplicates ensuring pairs are not reported twice if they were sequenced by the overlapping reads on both sides of the molecule (Supplementary Figure 1g).Parse2 distinguishes direct and unconfirmed contacts and reports the type, e.g., originating from the first side of paired-end read (R1), the second side (R2), present on both sides (R1&R2), or formed by end alignments on two sides (R1-R2) (Figure 2a, Supplementary Figure 1e–g). Parse2 reports this information and, therefore, allows filtration of only direct ligation events.Parse2 can combinatorially expand contacts of a walk by generating all pairs of DNA fragments of the same molecule, including indirect ones.Parse2 provides flexibility for reporting orientation and mapping positions of alignments in pairs. Pairs can either be reported either traditionally (as 5’-ends of alignments) or as ligation junctions (precise nucleotide-resolved positions of contacting alignment ends) (Supplementary Figure 1b).

31Parse2 preserves and reports the order of ligation events (a walk) within a molecule (Figure 2a).

32Parse2 removes internal duplicates ensuring pairs are not reported twice if they were sequenced by the overlapping reads on both sides of the molecule (Supplementary Figure 1g).

33Parse2 distinguishes direct and unconfirmed contacts and reports the type, e.g., originating from the first side of paired-end read (R1), the second side (R2), present on both sides (R1&R2), or formed by end alignments on two sides (R1-R2) (Figure 2a, Supplementary Figure 1e–g). Parse2 reports this information and, therefore, allows filtration of only direct ligation events.

34Parse2 can combinatorially expand contacts of a walk by generating all pairs of DNA fragments of the same molecule, including indirect ones.

35Parse2 provides flexibility for reporting orientation and mapping positions of alignments in pairs. Pairs can either be reported either traditionally (as 5’-ends of alignments) or as ligation junctions (precise nucleotide-resolved positions of contacting alignment ends) (Supplementary Figure 1b).

36While parse now supports some of these features for backward compatibility reasons, we strongly recommend using parse2 to analyze multi-contact data, such as PORE-C and MC-3C (Figure 2a). When applied to regular Hi-C or Micro-C, parse2 captures more contacts than regular parse (Figure 2c) because it captures the contacts originating from walks. We observe a substantial gain (5–25%) of contacts for both 150 bps 3C+ reads and 250 bps long reads. A common approach is to include all pairwise contacts within each walk via combinatorial expansion. However, we do not recommend this method because we find evidence of less desirable properties of indirect contacts generated this way (see section on quality control below).

37Both parse and parse2 output contact tables in a text tab-separated format called.pairs, designed by the NIH 4DN consortium. As a text format,.pairs has several advantages over custom binary formats: (i) text tables can be processed in all programming languages, (ii) are easily piped between individual CLI tools, and (iii) have a set of highly efficient utilities for sorting (Unix sort), compression (bgzip/lz4), and random access (tabix /pairix ).

38Each tab-separated row in a.pairs file describes a single observed contact. The required columns contain the id of the read and the genomic locations of the two sites that formed the contact. Pairtools parse augments these data with the pair type (Supplemental Figure 1a–g) and optional columns with details of the genomic alignments supporting the contact.

39Headers of.pairs files can store metadata, which by default includes the names of columns and the description of chromosomes of the reference genome. To ensure data provenance, pairtools extends this standard header with (i) the header of the.sam files that stored the original alignments and (ii) the complete history of data processing, with a separate entry for each CLI command of pairtools that was applied to the file. Pairtools provides a set of CLI and Python functions to parse and manipulate the header. Pairtools header can generate a new header, validate an existing one, transfer it between.pairs files, as well as set new column names (Figure 3a).

40Sorting the pairs in contact tables facilitates data processing, and analyses as it (i) enables fast random access via indexing and (ii) simplifies the detection of duplicated pairs, as they end up in adjacent rows of the table. Pairtools sorts pairs in two steps. First, individual pairs are “flipped,” i.e., arranged so that the alignment with the lower coordinate in the genome comes first in the pair (Figure 3b). Flipping is necessary as the two sides of a DNA molecule are sequenced in random order. Flipping is performed by default during parsing and can be done manually by pairtools flip. Second, pairtools sort performs Unix-based sorting of pairs in contact tables according to their genomic positions (on chrom1, chrom2, pos1, pos2). This sorting scheme has multiple advantages: it arranges pairs in blocks corresponding to contacts between a given pair of chromosomes, separates within (cis) from between (trans) chromosome contacts, and facilitates access to unmapped and multi-mapped pairs.

41A key issue for sequencing-based protocols, including Hi-C and other 3C+, is that the same DNA product can be duplicated by PCR and then sequenced and reported more than once, thus introducing an error into their quantitative measurements. Pairtools provides a computationally efficient tool for duplicate removal called pairtools dedup (Figure 1b). It detects clusters of pairs with matching genomic coordinates and strands and removes all but one pair. Pairtools dedup additionally permits minor mismatches between coordinates to account for potential losses of a few nucleotides from the ends of a molecule. To enable such imperfect coordinate matching, pairtools dedup relies on a KD-tree-based fixed radius nearest neighbor search algorithm, implemented in scipy and scikit-learn. To reduce its memory footprint, pairtools dedup processes input data in overlapping chunks. Finally, pairtools dedup can take additional columns from a.pair file and require additional properties of pairs to match, such as type of contact (direct/indirect), presence of mutations, or phased haplotype ).

42Tracking duplicates with pairtools enables an estimate of library complexity, i.e. the total number of unique DNA molecules prior to PCR amplification, an important QC for 3C+. Library complexity can guide the choice of sequencing depth of the library and provide an estimate of library quality. To estimate library complexity, pairtools assumes that each sequencing read is randomly chosen with replacement from a finite pool of fragments in DNA library.

43In addition to supporting the parse-sort-dedup steps (Figure 1b) that are sufficient to build a minimalistic 3C+ processing pipeline, pairtools also provides tools to build feature-rich pipelines. This includes tools for automated QC reporting, filtering of high-quality 3C+ data, as well as merging of replicates and conditions into meta contact maps, required for a complete and convenient end-to-end data processing pipeline.

44Pairtools provides a collection of tools for the manipulation of tabular pairs data.

45pairtools select splits and subsets datasets according to arbitrary filter expressions (Figure 3c). These expressions are provided as Python functions, enabling expressive and powerful filters. Filters can include wildcard and regex matching on string columns, custom functions from 3rd-party libraries (for examples of advanced usage, see Sister-C, scsHi-C ), as well as filtering pairs to a given subset of chromosomes.

46pairtools sample can generate random subsets of pairs, e.g., to equalize coverage between libraries or assess the statistical robustness of analyses (Figure 3d).

47pairtools merge combines multiple input datasets into one; for pre-sorted inputs, it efficiently produces sorted outputs.

48Together, pairtools select and merge enable the split-apply-combine pattern for distributed data processing.

493C+ experiments have multiple failure modes and thus require tight quality controls. Many experimental issues can be inferred from the statistics of the resulting 3C+ data (for a detailed discussion, see ).

50A particularly rich source of information about 3C+ experiments is the decay of contact frequency with the genomic distance referred to as the P(s) curve or scaling (borrowing the physics terminology for power-law relationships). Scalings are used both to characterize mechanisms of genome folding and reveal issues with QC. For example, early flattening of the scaling in Micro-C revealed the importance of long cross-linkers. Scalings can also be used to determine that combinatorial expansion of walks produces undesirable contacts because indirect contacts result in flatter scaling (Supplementary Figure 2e).

51Pairtools scaling calculates strand-oriented scalings that can be used for by-product quality control and filtration (Figure 3e). After the ligation step, some fragments can form a valid pair or produce unwanted 3C+ by-products, such as self-circles, dangling ends (unligated DNA) (Supplementary Figure 2c), and mirror reads (potential PCR artifacts). A short-range peak in divergent orientation is a sign of self-circled DNA, while a short-range peak in convergent orientation is a sign of dangling ends (Figure 3e, Supplementary Figure 2d). For example, early Hi-C variants with a low concentration of cross-linker caused the prevalence of self-circles. At larger genomic separations, pairs are formed in all four orientations with equal probabilities, and strand-oriented scalings converge. The orientation convergence distance indicates the minimum distance where pairs can simply be interpreted as contacts for a 3C+ dataset. Removing contacts below the orientation convergence distance removes nearly all by-products marked by restriction fragment annotation (see below, Supplementary Figure 2b). For DpnII Hi-C datasets, orientation convergence usually occurs by ~2kb. We note that for analyses downstream of QC, scaling can also be calculated from corrected binned data, e.g., using cooltools.

52For convenience and workflow reproducibility, pairtools stats automatically reports contact scalings for individual chromosomes. It also generates additional summary statistics, including the total number of pairs based on their type, and the number of trans contacts. This has been used to understand the impact of various protocol decisions. For example, information about the frequency of trans and different ranges of cis-contacts demonstrated that extra cross-linking yields more intra-chromosomal contacts. The frequency of contacts between the nuclear and mitochondrial genomes reflects the noise introduced by various digestion strategies. pairtools stats produces a human-readable nested dictionary of statistics stored in a YAML file or a tab-separated text table (used in, ) (Figure 3f). These outputs can be visualized with MultiQC, an interactive web-based tool that aggregates a wide set of sequencing QC statistics and provides an overview of whole collections of samples.

53Chromosome capture is a growing field, with novel protocol modifications emerging regularly. Thanks to its flexible architecture, pairtools can process data from many such experiments. For example, data from chemical modification-based protocols, such as scsHi-C, sn-m3C-seq, or Methyl-HiC can be processed by (i) extracting sequence mismatches into separate columns of.pairs by pairtools parse and (ii) filtering and subsetting pairs based on these columns with pairtools select (Figure 3c). For other popular and actively developing protocol variants, such as Micro-C, haplotype-resolved and single-cell Hi-C, pairtools provides specialized utilities.

54Many 3C+ protocols, particularly original Hi-C, rely on cutting DNA by restriction enzymes and theoretically should generate ligations only between restriction sites. Thus, early 3C+ analysis pipelines included filters that detected and removed (i) unligated single restriction fragments and (ii) ligations between pieces of DNA located far from any restriction sites. Pairtools restrict enables such filters by assigning the nearest restriction sites to each alignment in a pair (Figure 3g).

55However, we find restriction-based filters unnecessary for more recently-published Hi-C and do not include them in the standard pairtools pipeline. First, in our tests of recently published Hi-C datasets, the statistical properties of pairs located far from and close to restriction sites proved nearly the same (Supplementary Figure 2a). Second, we found that unligated pieces of DNA can be removed by a simpler filter against short-distance pairs, which can be calibrated using strand-oriented scalings (Supplementary Figure 2b).

56For downstream analyses in cooler, such by-products are removed by dropping pairs of bins with separations below a cutoff distance, which corresponds to removing a few central diagonals of a binned contact matrixf. Finally, the annotation of restriction sites becomes less accurate and unproductive for libraries generated with frequent and/or flexible cutters (e.g., DpnII, MboI, and DdeI), cocktails thereof, and impossible in restriction-free protocols, such as Micro-C and DNase Hi-C.

57Haplotype-resolved Hi-C leverages sequence variation between homologous chromosomes to resolve their contacts in cis and trans. In particular, single nucleotide variants (SNVs) can be used to differentiate contacts on the same chromosome (cis-homologous) from contacts between two homologs (trans-homologous).

58Pairtools phase is designed to resolve reads mapped to maternal and/or paternal genomes (Figure 3h). To enable analyses with pairtools phase, reads must be mapped to a reference genome that contains both haplotypes, reporting suboptimal alignments; these suboptimal alignments will be extracted into separate.pairs columns by pairtools parse. By checking the scores of the two best suboptimal alignments, pairtools phase distinguishes true multi-mappers from unresolved pairs (i.e., cases without a distinguishing SNV/indel) and reports the phasing status of each alignment in a pair: non-resolved, resolved as the first haplotype or the second haplotype, or multi-mapper. Further downstream, pairtools select and pairtools stats enable splitting and analyzing phased pairs. Pairtools phase has been successfully used to study allele-specific folding of chromosomes in D. melanogaster and mouse.

59Single-cell 3C+ experimental approaches shed light on variation and regularities in chromatin patterns among individual cells. Single-cell 3C+ data can be processed with pairtools almost the same way as bulk Hi-C, with the addition of one filter. In single cells, the number of contacts produced by each genomic site is limited by the chromosome copy number. Thus, regions with irregularly high coverage indicate amplification artifacts or copy number variations and must be excluded from downstream analyses. Pairtools filterbycov calculates genome-wide coverage per cell and excludes regions with coverage exceeding the specified threshold. This procedure helped to remove regions with anomalous coverage in single-cell Hi-C studies in Drosophila.

60Contact extraction from raw sequencing data is the first and typically most time-consuming step of the 3C+ data processing. Pairtools is one of the fastest methods (lagging behind only Chromap ) without consuming more memory (in combination with bwa mem), making it the best candidate for scalable 3C+ data processing (Figure 4). Notably, Pairtools is the only tool that combines high performance with the flexibility to enable adaptations to a broad range of 3C+protocols (Table 1).

61Pairtools provides a set of interoperable and high-performance utilities for processing 3C+ contact data, particularly for extraction of contacts from sequence alignments, manipulating pair files, deduplication, data filtering, generation of summary statistics and contact scalings, quality control, and treatment of data generated with 3C+ protocol modifications.

62pairtools is easy to install via pip and conda. We provide extensive documentation of pairtools, including example scripts of minimal pairtools-based pipelines and Jupyter tutorials with explanations of the working of pairtools, restriction, and phasing in pairtools GitHub repository.

63The modular structure of pairtools and its usage of the.pairs format already make it useful in many pipelines. pairtools is used in the 4DN pipeline (standard Hi-C), the PORE-C pipeline (multi-way Hi-C), HI-CAR nf-core pipeline (open-chromatin-associated contacts), and iMARGI pipelines (RNA-DNA contacts). pairtools also serve as the foundation of distiller, a feature-rich fastq-to-cooler pipeline, based on the nextflow workflow framework and maintained by Open2C. Distiller takes advantage of piping pairtools command outputs and can parallelize 3C+ data processing within a single machine, on a cluster, or in the cloud.

64In the future, a binary storage format for pairs could substantially speed up 3C+ contact extraction. Currently, zarr is the best candidate as it allows variable length strings (not supported by hdf5 ) and allows appending columns and storing multiple tables in a single file (not supported by parquet ).

65To summarize, Pairtools provides an adaptable framework for future development to enable the expanding universe of 3C+ protocols.