Requirements¶

Install using conda¶

We highly recommend using the conda package manager to install pre-compiled pairtools together with all its dependencies. To get it, you can either install the full Anaconda Python distribution or just the standalone conda package manager.

With conda, you can install pre-compiled pairtools and all of its dependencies from the bioconda channel:

$ conda install -c conda-forge -c bioconda pairtools

Install using pip¶

Alternatively, compile and install pairtools and its Python dependencies from PyPI using pip:

$ pip install pairtools

Install the development version¶

Finally, you can install the latest development version of pairtools from github. First, make a local clone of the github repository:

$ git clone https://github.com/mirnylab/pairtools

Then, you can compile and install pairtools in the development mode, which installs the package without moving it to a system folder and thus allows immediate live-testing any changes in the python code. Please, make sure that you have cython installed!

$ cd pairtools
$ pip install -e ./

Overview¶

Hi-C experiments aim to measure the frequencies of contacts between all pairs of loci in the genome. In these experiments, the spacial structure of chromosomes if first fixed with formaldehyde crosslinks, after which DNA is partially digested with restriction enzymes and then re-ligated back. Then, DNA is shredded into smaller pieces, released from nucleus, sequenced and aligned to the reference genome. The resulting sequence alignments reveal if DNA molecules were formed through ligations between DNA from different locations in the genome. These ligation events imply that ligated loci were close to each other when the ligation enzyme was active, i.e. they formed "a contact".

pairtools parse detects ligation events in the aligned sequences of DNA molecules formed in Hi-C experiments and reports them in the .pairs/.pairsam format.

Terminology¶

Throughout this document we will be using the same visual language to describe how DNA sequences (in the .fastq format) are transformed into sequence alignments (.sam/.bam) and into ligation events (.pairs).

Unmapped/multimapped reads¶

Sometimes, one side or both sides of a read pair may not align to the reference genome:

Multiple ligations (walks)¶

Finally, a read pair may contain more than two alignments:

Interpreting gaps between alignments¶

Reads that are only partially aligned to the genome can be interpreted in two different ways. One possibility is to assume that this molecule was formed via at least two ligations (i.e. it's a walk) but the non-aligned part (a gap) was missing from the reference genome for one reason or another. Another possibility is to simply ignore this gap (for example, because it could be an insertion or a technical artifact), thus assuming that our molecule was formed via a single ligation and has to be reported:

Rescuing single ligations¶

Importantly, some of DNA molecules containing only one ligation junction may still end up with three alignments:

Sorting¶

pairtools sort arrange pairs in the order of (chrom1, chrom2, pos1, pos2). This order is also known as block sorting, because all pairs between any given pair of chromosomes become grouped into one continuous block. Additionally, pairtools sort also sorts pairs with identical positions by pair_type. This does not really do much for mapped reads, but it nicely splits unmapped reads into blocks of null-mapped and multi-mapped reads.

We note that there is an alternative to block sorting, called row sorting, where pairs are sorted by (chrom1, pos1, chrom2, pos2). In pairtools sort, we prefer block-sorting since it cleanly separates cis interactions from trans ones and thus is a more optimal solution for typical use cases.

Flipping¶

In a typical paired-end experiment, side1 and side2 of a DNA molecule are defined by the order in which they got sequenced. Since this order is essentially random, any given Hi-C pair, e.g. (chr1, 1.1Mb; chr2, 2.1Mb), may appear in a reversed orientation, i.e. (chr2, 2.1Mb; chr1, 1.1Mb). If we were to preserve this order of sides, interactions between same loci would appear in two different locations of the sorted pair list, which would complicate finding PCR/optical duplicates.

To ensure that Hi-C pairs with similar coordinates end up in the same location of the sorted list, we flip pairs, i.e. we choose side1 as the side with the lowest genomic coordinate. Thus, after flipping, for trans pairs (chrom1!=chrom2), order(chrom1)<order(chrom2); and for cis pairs (chrom1==chrom2), pos1<=pos2. In a matrix representation, flipping is equal to reflecting the lower triangle of the Hi-C matrix onto its upper triangle, such that the resulting matrix is upper-triangular.

In pairtools, flipping is done during parsing - that's why pairtools parse requires a .chromsizes file that specifies the order of chromosomes for flipping. Importantly, pairtools parse also flips one-sided pairs such that side1 is always unmapped; and unmapped pairs such that side1 always has a "poorer" mapping type (i.e. null-mapping<multi-mapping).

Chromosomal order for sorting and flipping¶

Importantly, the order of chromosomes for sorting and flipping can be different. Specifically, pairtools sort uses the lexicographic order for chromosomes (chr1, chr10, chr11, ..., chr2, chr21,...) instead of the "natural" order (chr1, chr2, chr3, ...); at the same time, flipping is done in pairsamtools parse using the chromosomal order specified by the user.

pairtools sort uses the lexicographic order for sorting chromosomes. This order is used universally to sorting strings in all languages and tools [1], which makes it easy to design tools for indexing and searching in sorted pair lists.

At the same time, pairtools parse uses a custom user-provided chromosomal order to flip pairs. For performance considerations, for flipping, we recommend ordering chromosomes in a way that will be used in plotting contact maps.

[1]

Unfortunately, many existing genomes use rather unconventional choices in chromosomal naming schemes. For example, in sacCer3, chromosomes are enumerated with Roman numerals; in dm3, big autosomes are split 4 different contigs each. Thus, it is impossible to design a universal algorithm that would order chromosomes in a "meaningful" way for all existing genomes.

.pairs¶

.pairs is a simple tabular format for storing DNA contacts detected in a Hi-C experiment. The detailed .pairs specification is defined by the 4DN Consortium.

The body of a .pairs contains a table with a variable number of fields separated by a "\t" character (a horizontal tab). The .pairs specification fixes the content and the order of the first seven columns:

A .pairs file starts with a header, an arbitrary number of lines starting with a "#" character. By convention, the header lines have a format of "#field_name: field_value". The .pairs specification mandates a few standard header lines (e.g., column names, chromosome order, sorting order, etc), all of which are automatically filled in by pairtools.

The entries of a .pairs file can be flipped and sorted. "Flipping" means that the sides 1 and 2 do not correspond to side1 and side2 in sequencing data. Instead, side1 is defined as the side with the alignment with a lower sorting index (using the lexographic order for chromosome names, followed by the numeric order for positions and the lexicographic order for pair types). This particular order of "flipping" is defined as "upper-triangular flipping", or "triu-flipping". Finally, pairs are typically block-sorted: i.e. first lexicographically by chrom1 and chrom2, then numerically by pos1 and pos2.

Pairtools' flavor of .pairs¶

.pairs files produced by pairtools extend .pairs format in a few ways.

Pair types¶

pairtools use a simple two-character notation to define all possible pair types, according to the quality of alignment of the two sides. The type of a pair can be defined unambiguously using the table below. To use this table, identify which side has an alignment of a "poorer" quality (unmapped < multimapped < unique alignment) and which side has a "better" alignment and find the corresponding row in the table.

[1]

"walks", or, C-walks are Hi-C molecules formed via multiple ligation events which cannot be reported as a single pair.

[2]

"corrupt" pairs are those with technical issues - e.g. missing a FASTQ sequence/SAM entry from one side of the molecule.

[2]

"rescued" pairs have two non-overlapping alignments on one of the sides (referred below as the chimeric side/read), but the inner (3'-) one extends the only alignment on the other side (referred as the non-chimeric side/read). Such pairs form when one of the two ligated DNA fragments is shorter than the read length. In this case, one of the reads contains this short fragment entirely, together with the ligation junction and a chunk of the other DNA fragment (thus, this read ends up having two non-overlapping alignments). Following the procedure introduced in HiC-Pro and Juicer, pairtools parse rescues such Hi-C molecules, reports the position of the 5' alignment on the chimeric side, and tags them as "NU", "MU", "UR" or "RU" pair type, depending on the type of the 5' alignment on the chimeric side. Such molecules can and should be used in downstream analysis. Read more on the rescue procedure in the section on parsing.

[3]

pairtools dedup detects molecules that could be formed via PCR duplication and tags them as "DD" pair type. These pairs should be excluded from downstream analyses.

.pairsam¶

pairtools also define .pairsam, a valid extension of the .pairs format. On top of the pairtools' flavor of .pairs, .pairsam format adds two extra columns containing the alignments from which the Hi-C pair was extracted:

Note that, normally, the fields of a sam alignment are separated by a horizontal tab character (\t), which we already use to separate .pairs columns. To avoid confusion, we replace the tab character in sam entries stored in sam1 and sam2 columns with a UNIT SEPARATOR character (\031).

Finally, sam1 and sam2 can store multiple .sam alignments, separated by a string '\031NEXT_SAM\031'

.pairs format¶

The motivation behind some of the technical decisions in the pairtools' flavor of .pairs/.pairsam:

CLI¶