bio-long-read-sequencing-long-read-alignment

What this skill does

## Version Compatibility

Reference examples tested with: minimap2 2.28+, samtools 1.19+, winnowmap 2.03+, pbmm2 1.13+.

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

Version-driven behavior to record:
- `lr:hq` and `map-iclr` were added in minimap2 2.27; `lr:hqae` in 2.28. Use >=2.28.
- `--MD` was broken by the 2.27 `--ds` addition and fixed in 2.28; use >=2.28 for any MD-dependent caller.
- A prebuilt `.mmi` index bakes in k/w/H/I - it must be built with the same preset used for alignment.

If code throws an error, introspect the installed tool (`minimap2 --help`, man page) and adapt the example to the actual API rather than retrying.

# Long-Read Alignment with minimap2

**"Align my long reads to the reference"** -> Map with the preset that matches the reads' ERROR RATE (not just platform), keeping the supplementary alignments and tags that downstream callers need.
- CLI: `minimap2 -ax lr:hq --MD -Y ref.fa reads.fq | samtools sort -o aln.bam` (accurate ONT/R10), `minimap2 -ax map-ont` (noisy R9 ONT), `minimap2 -ax map-hifi` (PacBio HiFi)

## The Single Most Important Modern Insight -- The Preset Rewrites the Scoring Model, So the Wrong One Fabricates or Erases Variants

`-x <preset>` is not a label. The man page defines each preset as a literal bundle that rewrites k-mer/window AND the entire scoring model (match `-A`, mismatch `-B`, gap-open `-O`, gap-extend `-E`), Z-drop `-z`, and chaining bandwidth `-r`. So the wrong preset does not merely "align worse" - it changes which gaps the chainer will span, and thereby fabricates or erases the exact insertions, deletions, introns, and SV breakpoints the downstream caller is built to find. Three corollaries an expert holds:

1. **Preset = read error rate, not platform.** "ONT" is no longer one regime: noisy R9/fast/hac = `map-ont`; accurate Q20+/duplex/R10-sup = `lr:hq` (2.27+, ~4x fewer CPU-hours, equal/better accuracy). `map-hifi` is literally `lr:hq` + HiFi scoring.
2. **Supplementary alignments ARE the SV signal.** SV callers read the split-read pattern (primary + supplementary chimeric pieces), not the tidy primary. Feeding them secondaries, or hard-clipping supplementaries, silently degrades SV sensitivity.
3. **A tag absent at alignment time is unrecoverable.** MM/ML, MD, cs - if minimap2 did not write them, no downstream tool can reconstruct them; the pipeline succeeds and produces empty/wrong results.

## Preset Taxonomy

| Preset | Read type / when correct | Notes |
|--------|--------------------------|-------|
| `map-ont` | ONT noisy genomic (R9, fast/hac) | the historic default; ~10% error scoring |
| `lr:hq` | accurate long reads <1% err (ONT Q20+/duplex/R10 sup) | 2.27+; the modern accurate-ONT default |
| `map-hifi` | PacBio HiFi/CCS genomic | = `lr:hq` + HiFi scoring (2.27+) |
| `map-pb` | PacBio CLR (legacy, ~15% err) | homopolymer-compressed minimizers; NEVER for HiFi |
| `splice` | noisy long RNA (ONT cDNA/direct RNA) | add `-uf` for stranded direct RNA |
| `splice:hq` | accurate long RNA (PacBio Iso-Seq, R10 cDNA) | |
| `asm5` / `asm10` / `asm20` | assembly-to-ref at ~0.1% / ~1% / ~5% divergence | PAF output; `--cs` for paftools call |
| `ava-ont` / `ava-pb` | all-vs-all read overlap (miniasm) | overlaps only, no base alignment |
| `lr:hqae` | accurate reads back to THEIR OWN assembly | 2.28+; fixes centromere self-mapping mismaps |

## Aligner Decision Tree

| Situation | Aligner | Why |
|-----------|---------|-----|
| Standard ONT/HiFi to a normal reference (SNV/SV/general) | minimap2 | the de-facto standard; default for Sniffles2, cuteSV, Clair3 |
| Accurate ONT (Q20+/duplex/R10 sup) | minimap2 `-x lr:hq` | ~4x faster than map-ont, equal/better |
| Centromeres / satellite arrays / segmental dups / T2T reference | Winnowmap2 | minimap2 minimizer-masking mismaps long tandem repeats; Winnowmap down-weights via meryl repetitive k-mers |
| Complex/nested SVs, inversions, tandem dups | VACmap (or lra) | variant-aware nonlinear chaining resolves CSVs minimap2 splits |
| Accurate reads -> a diploid assembly built from them | minimap2 `-x lr:hqae` (2.28+) | avoids self-assembly centromere mismaps |
| PacBio-native (.bam/.xml, want sorted+indexed in one call) | pbmm2 | minimap2 + PacBio plumbing; presets SUBREAD/CCS/HIFI/ISOSEQ |
| Legacy Sniffles1 reproduction | NGMLR | the 2018 standard, now superseded by minimap2+Sniffles2 |

## Tag Requirements by Downstream Tool

A missing tag is a silent failure. Add the tag at alignment time.

| Tag / flag | What it does | Needed for |
|------------|--------------|-----------|
| `--MD` | mismatch positions vs ref | many small-variant callers, IGV mismatch coloring (use minimap2 >=2.28) |
| `-Y` | soft-clip supplementary (default hard-clips) | SV callers: keeps breakpoint/insertion SEQ on the split read |
| `-y` | copy MM/ML (and other) tags from the input | methylation: carries Dorado MM/ML through alignment |
| `--cs` | minimap2 difference string | `paftools.js call` (assembly/long-read variant calling) requires it |
| `--eqx` | `=`/`X` CIGAR instead of `M` | tools that read match/mismatch from CIGAR |
| `-L` | move >65535-op CIGAR to CG:B tag | ultra-long ONT reads (else unrepresentable in BAM) |

Supplementary (flag 0x800) = split piece of one read across loci = the SV substrate, controlled by chaining + `-Y`. Secondary (flag 0x100) = multi-mapping alternative, controlled by `--secondary`/`-N`/`-p`. SV work keeps primary+supplementary and is fine with `--secondary=no`.

## Core Commands

```bash
# Accurate ONT (Q20+/R10 sup) -> genome, SV+variant ready, sorted+indexed
minimap2 -ax lr:hq -t 16 --MD -Y -R '@RG\tID:s1\tSM:s1' ref.fa reads.fq.gz \
  | samtools sort -@4 -o aln.bam && samtools index aln.bam

# Noisy ONT (R9 / fast / hac)
minimap2 -ax map-ont -t 16 --MD -Y ref.fa r9.fq.gz | samtools sort -o ont.bam

# PacBio HiFi (minimap2, or pbmm2 in one sorted+indexed call)
minimap2 -ax map-hifi -t 16 --MD -Y ref.fa hifi.fq.gz | samtools sort -o hifi.bam
pbmm2 align --preset HIFI --sort -j 16 ref.fa hifi.bam hifi.aligned.bam

# Methylation passthrough: carry Dorado MM/ML through alignment (the -y trap). -Y soft-clips
# supplementary records so hard-clipping does not break the MM per-base skip counting.
samtools fastq -T MM,ML dorado.mod.bam \
  | minimap2 -ax lr:hq -y -Y --MD ref.fa - \
  | samtools sort -o meth.bam        # then modkit pileup meth.bam ...

# Direct RNA (ONT): stranded forward-only, small k for terminal-exon sensitivity
minimap2 -ax splice -uf -k14 -G500k ref.fa dRNA.fq.gz | samtools sort -o drna.bam
#   -G500k raises max-intron above the 200k default only for genes with long introns

# Assembly-to-reference: PAF is correct here; --cs enables paftools variant calling
minimap2 -cx asm5 --cs ref.fa asm.fa > asm.paf
paftools.js call asm.paf > asm.var.vcf

# Repeats / centromeres / T2T: Winnowmap (precompute repetitive k-mers)
meryl count k=15 output merylDB ref.fa
meryl print greater-than distinct=0.9998 merylDB > repetitive_k15.txt
winnowmap -W repetitive_k15.txt -ax map-ont ref.fa reads.fq.gz | samtools sort -o wm.bam

# Prebuild index - bake the SAME preset's k/w in (else the preset's k/w is ignored)
minimap2 -x lr:hq -d ref.lrhq.mmi ref.fa
```

## Per-Method Failure Modes

### Hard-clipped supplementaries break SV insertion calls
**Trigger:** mapping for SV calling without `-Y`. **Mechanism:** minimap2 hard-clips supplementary records, discarding the breakpoint-spanning bases. **Symptom:** imprecise/missing insertions and translocations. **Fix:** add `-Y` (soft-clip) so split reads keep full SEQ.

### Methylation tags silently dropped
**Trigger:** aligning a Dorado mod BAM without preserving tags. **Mechanism:** `samtools fastq` strips MM/ML unless `-T MM,ML`; minimap2 ignores them unless `-y`. **Symptom:** aligned BAM has no MM/ML; modkit produces empty bedMethyl, no error. **Fix:** `samtools fastq -T MM,ML | minimap2 -y -Y`, or us