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

name: bio-long-read-sequencing-long-read-alignment description: Aligns Oxford Nanopore and PacBio long reads (and assemblies) to a reference with minimap2 using the error-rate-matched preset (map-ont, lr:hq, map-hifi, map-pb, splice/splice:hq, asm5/10/20, ava), producing a sorted/indexed BAM for variant, SV, methylation, or isoform analysis. Covers why the preset rewrites the scoring/chaining model, why SV calling rides on supplementary not secondary alignments, carrying MM/ML methylation tags through with -y, the multi-part-index MAPQ trap, and when to swap in Winnowmap/VACmap/lra/pbmm2. Use when mapping ONT or PacBio reads, choosing a minimap2 preset by platform/chemistry, preparing input for Clair3/medaka/Sniffles/modkit, aligning into repeats/centromeres, or spliced-aligning cDNA/Iso-Seq. tool_type: cli primary_tool: minimap2

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

Aligner Decision Tree

Tag Requirements by Downstream Tool

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

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

# 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 use dorado aligner.

Multi-part index destroys MAPQ

Trigger: reference larger than -I (default 8G) - large plant/polyploid or concatenated refs. Mechanism: minimap2 builds a multi-part index and scores batches independently, so cross-batch best hits are invisible and MAPQ is wrong. Symptom: "no @SQ lines ... use --split-prefix"; spurious MAPQ. Fix: -I <bigger-than-ref> or --split-prefix.

Wrong preset on accurate reads

Trigger: map-ont on Q20/R10/duplex, or map-pb on HiFi. Mechanism: noisy-read scoring on accurate reads (or CLR scoring on HiFi). Symptom: ~4x slower for no gain (map-ont case), or spurious clips/indels (map-pb-on-HiFi). Fix: lr:hq for accurate ONT, map-hifi for HiFi.

Direct-RNA junctions on the wrong strand

Trigger: -ax splice on direct RNA without -uf. Mechanism: splice defaults to -ub (GT-AG on both strands), but dRNA is stranded. Symptom: invented/misplaced introns. Fix: add -uf (and usually -k14).

Centromere/SD mismapping looks fine in flagstat

Trigger: plain minimap2 into long tandem repeats. Mechanism: minimizer masking collapses minimizer density, so reads map to the wrong paralog/copy. Symptom: reads still "map" (flagstat clean) but produce false SVs/heterozygosity in repeats. Fix: Winnowmap2 with a meryl repetitive-k-mer set.

Quantitative Thresholds

Common Errors

References

• Li H. 2018. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34(18):3094-3100.
• Li H. 2021. New strategies to improve minimap2 alignment accuracy. Bioinformatics 37(23):4572-4574.
• Jain C, Rhie A, Hansen NF, et al. 2022. Long-read mapping to repetitive reference sequences using Winnowmap2. Nat Methods 19:705-710.
• Ren J, Chaisson MJP. 2021. lra: a long read aligner for sequences and contigs. PLoS Comput Biol 17(6):e1009078.
• Ding H, et al. 2026. VACmap: an accurate long-read aligner for unraveling complex genomic rearrangements. Nat Commun 16:11198.
• Sedlazeck FJ, et al. 2018. Accurate detection of complex structural variations using single-molecule sequencing (NGMLR/Sniffles). Nat Methods 15:461-468.

Related Skills

• basecalling - The basecaller chemistry/error rate that picks the preset; carries MM/ML to pass with -y
• long-read-qc - Read length/quality before mapping; % identity from the aligned BAM
• structural-variants - Consumes the supplementary (split-read) signal this preserves with -Y
• clair3-variants - Small-variant calling on this BAM (needs the matched basecaller model)
• nanopore-methylation - Pileup of the MM/ML tags carried through with -y
• isoseq-analysis - Spliced alignment of full-length cDNA/Iso-Seq
• alignment-files/sam-bam-basics - Sort/index/inspect the BAM this produces
• alignment-files/alignment-filtering - Filter by MAPQ and secondary/supplementary flags
• genome-assembly/long-read-assembly - Assemble the reads instead of reference-mapping