bio-hi-c-analysis-contact-pairs

What this skill does

## Version Compatibility

Reference examples tested with: pairtools 1.1+, bwa 0.7.17+ (or bwa-mem2 2.2+), chromap 0.2+, samtools 1.19+, cooler 0.10+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

pairtools defaults have shifted across releases (e.g. `parse --max-molecule-size` is 750 bp in 1.1.x; `dedup --backend` defaults to scipy). `parse` and `parse2` report DIFFERENT positions by default - confirm `--report-position` before mixing outputs. If a command errors, introspect with `pairtools <subcommand> --help` and adapt rather than retrying.

# Hi-C Contact Pairs

**"Turn my Hi-C reads into clean contacts and tell me if the library worked"** -> Align mates independently through ligation junctions, classify and deduplicate pairs to a 5'-canonical .pairs file, then read the cis/trans and orientation statistics to decide library quality before any matrix is built.
- CLI: `bwa mem -SP5M ref.fa R1.fq R2.fq | pairtools parse -c chrom.sizes | pairtools sort | pairtools dedup | pairtools stats`

## The Single Most Important Modern Insight -- The Read Count Is a Lie Until Pairs Are Classified; Library Quality Is Decided in pairtools, Not in the Aligner

A Hi-C library's *usable signal* is not "reads sequenced" - it is the **uniquely-mapped, deduplicated, long-range cis** contacts. Everything between FASTQ and the matrix exists to strip a *specific* artifact of proximity-ligation chemistry, and the diagnostic ratios from `pairtools stats` reveal whether the experiment succeeded **before** any compute is spent binning it. Three load-bearing consequences:

1. **% long-range cis is the one-number quality metric; trans is the noise floor.** True crosslink-ligation contacts are overwhelmingly cis and distance-decaying. Random ligation between two unrelated molecules in solution is as likely to be trans as cis-far, so **trans% is a direct readout of the spurious-ligation floor**. A good in-situ human library runs cis>=1kb ~50-65%, inter-chromosomal <10%. But these numbers are **genome-size dependent** - a yeast or bacterial genome legitimately has higher expected trans (more inter-chromosomal volume per cis distance). Never apply a human trans threshold to a microbe.

2. **The ligation junction lives INSIDE the read, so a local/split aligner aligning mates independently is required.** A single read often sequences *through* a ligation junction (locus A | locus B within one read). An end-to-end aligner soft-clips or mis-maps it and the contact is lost. `bwa mem -SP5M` aligns R1 and R2 as **independent single-end reads** (`-SP` skips mate rescue and pairing - proper-pair logic would destroy every long-range and trans contact) and marks the **5'-most chimeric segment primary** (`-5`, the anchor for pairtools' 5' convention). The chimera fraction rises with read length, so on 150bp PE and on Micro-C long reads this is a large, real chunk of contacts.

3. **Keep UU AND rescued UC; selecting only UU silently discards every rescued ligation.** A naive `select pair_type=="UU"` throws away the chimeric reads pairtools successfully reconstructed (UC = combined-unique) - a meaningful fraction on long reads. The 4DN standard keeps **UU and UC**.

## Aligner Taxonomy

| Aligner | Role | Hi-C invocation | When |
|---------|------|-----------------|------|
| bwa mem -SP5M | reference standard; local/split alignment reconstructs in-read junctions | `bwa mem -SP5M -t N ref.fa R1 R2` | default; best inter-contig accuracy in benchmarks |
| bwa-mem2 | drop-in faster reimplementation, identical output, same flags | `bwa-mem2 mem -SP5M -t N ref.fa R1 R2` | when speed matters and the larger index fits RAM |
| chromap --preset hic | ultrafast integrated align + dedup + pairs (4DN .pairs out) | `chromap --preset hic -x idx -r ref.fa -1 R1 -2 R2 -o out.pairs` | ~10x faster scans; trades fine walk-policy control for speed |

`-SP5M`, letter by letter: `-S` skip mate rescue; `-P` skip pairing (no proper-pair rescue); together `-SP` align mates as independent single-end reads. `-5` mark the 5'-most split segment primary (anchors the 5' convention). `-M` is legacy compatibility only (secondary flag 256 vs supplementary 2048); pairtools handles either - never agonize over `-M`, never drop `-SP5`.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Standard in-situ/Omni-C, FASTQ -> matrix | bwa mem -SP5M -> parse (5unique) -> sort -> dedup -> stats | the 4DN/distiller default; restriction-agnostic |
| Fastest scan, control not critical | chromap --preset hic | integrated align+dedup+pairs ~10x faster |
| Multi-way contacts (Pore-C, MC-3C, Micro-C walks) | `parse --walks-policy all` or `parse2 --expand` | 5unique COLLAPSES concatemers to pairwise silently |
| Micro-C / DNase Hi-C | NO fragment map; do NOT apply a 1kb min-distance cut | sub-1kb (nucleosome ladder) is signal, not artifact |
| Arima / Hi-C 3.0 dual-enzyme, fragment-level | fragment file must encode BOTH motifs (`restrict -f`) | a single-enzyme digest file is silently wrong |
| Repeat-heavy genome / stringent loop anchors | raise `parse --min-mapq 30` | reclassifies borderline reads M, drops repeat false anchors |
| Allele-specific / diploid folding | diploid ref + `-XA` suboptimal hits -> `pairtools phase` -> two coolers | needs the suboptimal-score gap to resolve haplotypes |
| Decide whether to sequence deeper | complexity curve from dup model / preseq lc_extrap | a bare dup% without depth is meaningless |
| Build the matrix from clean pairs | -> hic-data-io (`cooler cload pairs`) | binning happens after classification/dedup |
| Annotate boundary/anchor coordinates | -> genome-intervals/bed-file-basics | pairs are 1-based, half-open conventions differ |

## Align: Mates as Independent Single-End Reads

```bash
bwa index ref.fa                                            # or bwa-mem2 index ref.fa
bwa mem -SP5M -t 16 ref.fa R1.fq.gz R2.fq.gz | \            # -SP: independent SE; -5: 5' segment primary
    samtools view -b -@ 8 - > aligned.bam
# chromap fast path (integrated align + dedup + 4DN pairs, no separate pairtools needed):
# chromap -i -r ref.fa -o idx && \
# chromap --preset hic -x idx -r ref.fa -1 R1.fq.gz -2 R2.fq.gz -o sample.pairs
```

## Parse, Sort, Dedup, Select: the pairtools Core

```bash
# Parse alignments into a 5'-canonical .pairs. min-mapq 1 (default) is permissive: only MAPQ-0 is "multi".
pairtools parse -c chrom.sizes --walks-policy 5unique --min-mapq 1 \
    --add-columns mapq --drop-sam aligned.bam | \
pairtools sort --nproc 8 | \                                # block-sort; flips to upper-triangular (5'-canonical)
pairtools dedup --max-mismatch 3 --mark-dups \              # within 3bp on both sides = duplicate; tag DD
    --output-stats sample.dedup.stats | \
pairtools select '(pair_type=="UU") or (pair_type=="UC")' \ # keep both-unique AND rescued chimeric
    -o sample.valid.pairs.gz
```

Dedup MUST run on a flipped, 5'-canonical file (`sort` does the flip); on a non-canonical file dedup under-collapses and dup% reads falsely low. `--max-molecule-size` (750 bp in 1.1.x) governs single-ligation chimera rescue; `--max-inter-align-gap` (20 bp) sets when a coverage gap becomes a null alignment.

## Library QC: the Decision This Skill Owns

`pairtools stats` is the canonical readout. Read it as a funnel, not a single number.

```bash
pairtools stats --bytile-dups -o sample.stats.tsv sample.valid.pairs.gz
# Key fields: frac_dups; frac_cis; cis_1kb+/cis_20kb+; trans; pair_types; dist_freq orientation FF/FR/RF/RR.
```

- **% long-range cis (cis>=1kb, often cis>=20kb)** = signal quality. **trans = noise floor** (genome-size dependent).
- **Orientation vs distance = fragment-map-free dangling-end/self-circle QC.** Above ~1kb the four orientations FF/FR/RF/RR each converge to ~25% (random, the positive QC signal). A short-range **FR (inward) spike** = dangling ends /