From de60e8c9f23e403699d0a78091300beb4984b04b Mon Sep 17 00:00:00 2001 From: Jay Hesselberth Date: Mon, 5 Oct 2026 06:14:40 -0600 Subject: [PATCH] feat: pair orientation, MQ, discordant count, fixmate-exact TLEN Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 8 +- CLAUDE.md | 4 +- crates/eschalign-cli/src/align.rs | 103 ++++- crates/eschalign-cli/tests/align_paired.rs | 386 ++++++++++++++++-- crates/eschalign-cli/tests/align_unpaired.rs | 100 ----- .../tests/fixtures/pairs/gen_pair_fixture.py | 2 + .../tests/fixtures/pairs/pairs_R1.fq | 4 + .../tests/fixtures/pairs/pairs_R2.fq | 4 + crates/eschalign/src/lib.rs | 3 +- crates/eschalign/src/mapper.rs | 2 +- crates/eschalign/src/mapper/pair.rs | 136 +++++- docs/align.md | 6 +- 12 files changed, 576 insertions(+), 182 deletions(-) delete mode 100644 crates/eschalign-cli/tests/align_unpaired.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 6434479..09658ce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,10 @@ per-mate `NM`/`MD`/`AS`/`XS`/`XA`, MAPQ 60 only for a unique best reference, mates adjacent in input order, `--sort coordinate` and `--secondary` supported. See [docs/align.md](docs/align.md#paired-end). +- `--pair-orientation fr|rf|ff` (default `fr`): how the mates must face each + other in a proper pair. `MQ` (the mate's MAPQ) is written beside `MC`, and the + run summary counts discordant pairs (mates whose best references differ, + each reported at its own) (#4). - 3' adapter read-through is detected before alignment (TruSeq Read 1/2 prefixes; `--adapter-r1`, `--adapter-r2`, `--min-adapter-overlap`, `--no-adapter-trim`), the adapter is kept as soft clip with an `ad:i` tag, and @@ -22,7 +26,7 @@ - Library: `eschalign::adapter`, `Aligner::map_pairs`/`map_pairs_scored`, `PairOptions`, `PairMapping`, `PairHit`. - Tests: `align_paired.rs` over a synthesized pair fixture in - `tests/fixtures/pairs/`; `align_unpaired.rs` pins the unpaired output to a + `tests/fixtures/pairs/`; `align_paired::single_end_output_unchanged` pins the unpaired output to a golden record-stream digest. ### Changed @@ -30,7 +34,7 @@ - A BAM whose first record is flagged paired (0x1 with 0x40/0x80) is now a paired run; it used to be read as unpaired reads with the flags discarded. Every other input, and every existing flag, behaves as before (output - byte-identical, `unpaired_output_unchanged`). + byte-identical, `single_end_output_unchanged`). ## 0.1.0 diff --git a/CLAUDE.md b/CLAUDE.md index 5743af3..72ce19b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -183,7 +183,7 @@ would treat it as on. `clippy.toml`'s `disallowed-methods` enforces this (and - `crates/eschalign/` — alphabet / `Scoring` / `Mode` / scalar oracle / `simd` (two kernel families) / `sam` (CIGAR, `MD`/`NM`) / `panel` / `mapper` (`Aligner`) / `cuda` (`gpu`) / `fasta` / `error`. See its `CLAUDE.md` for the map, tests and traps. - `crates/eschalign-cli/src/` — `main.rs` (clap, tracing, rayon pool, allocator), `align.rs` (the pipeline), `device.rs` (`--device`), `env.rs`, `threads.rs`, `progress.rs`, `style.rs`, `uuid_util.rs`. -- `crates/eschalign-cli/tests/` — `align_e2e.rs`, `align_sort.rs`, `align_unpaired.rs` over `tests/fixtures/` (60 real tRNA reads as a uBAM-shaped BAM plus the 47-reference panel and an ambiguity-code variant; **no pod5** — these tests never read signal), and `align_paired.rs` over `tests/fixtures/pairs/` (a synthesized 4-reference panel with two paralogs and one hand-built read pair per behaviour; `gen_pair_fixture.py` regenerates it). `align_unpaired.rs` pins the unpaired output to a golden digest of the uncompressed record stream (`fixtures/unpaired_golden.tsv`, `ESCHALIGN_BLESS=1` rewrites it) — a change to it is a change to unpaired behaviour and must be intended. +- `crates/eschalign-cli/tests/` — `align_e2e.rs`, `align_sort.rs` over `tests/fixtures/` (60 real tRNA reads as a uBAM-shaped BAM plus the 47-reference panel and an ambiguity-code variant; **no pod5** — these tests never read signal), and `align_paired.rs` over `tests/fixtures/pairs/` (a synthesized 4-reference panel with two paralogs and one hand-built read pair per behaviour; `gen_pair_fixture.py` regenerates it). `align_paired::single_end_output_unchanged` pins the single-end output to a golden digest of the uncompressed record stream (`fixtures/unpaired_golden.tsv`, `ESCHALIGN_BLESS=1` rewrites it) — a change to it is a change to unpaired behaviour and must be intended. - `crates/eschalign/tests/fixtures/` — the parasail and calmd goldens and the scripts that regenerate them (`gen_align_golden.py` needs `parasail`; `gen_calmd_golden.py` needs samtools). ### The command @@ -192,7 +192,7 @@ would treat it as on. `clippy.toml`'s `disallowed-methods` enforces this (and ### Paired-end (#1) -`-1/-2`, `--interleaved`, or a BAM whose first primary record has 0x1 with 0x40/0x80 (`bam_is_paired` opens it once to peek) selects `Input::Paired` and a `PairCfg`; **nothing else changes for unpaired input** (`align_unpaired::unpaired_output_unchanged`). Where the pieces live: +`-1/-2`, `--interleaved`, or a BAM whose first primary record has 0x1 with 0x40/0x80 (`bam_is_paired` opens it once to peek) selects `Input::Paired` and a `PairCfg`; **nothing else changes for unpaired input** (`align_paired::single_end_output_unchanged`). Where the pieces live: - **Library** (`eschalign`, still std only): `adapter.rs` (`find`, `detect_pair`: leftmost prefix match, `overlap/10` mismatches, the two mates' insert ends must agree or the weaker detection is dropped) and `mapper/pair.rs` (`Aligner::map_pairs_scored`, `PairOptions`, `PairMapping::{Proper, Separate}`). A *candidate* is `(reference, template strand)`; both mates must clear `min_score`; candidates are traced **tier by tier in descending summed score, every open pair's next tier through one `align_pairs` call** (so the traceback lanes still fill), and the first tier with a candidate that passes the constraints on the *traced* coordinates (no outie, insert ≤ `max_insert`, both mates align a base) is the tie set. Constraints need coordinates, so a score alone can never decide a pair. No qualifying candidate → each mate goes through `map_reads_scored` with `both_strands: true`. - **CLI** (`align.rs`): `PairSource` (`Bam` / `Interleaved` / `Files`) → `read_pair_batches` (whole pairs per batch, mate 1 then mate 2 adjacent, names stripped of `/1` `/2`, orphan/mismatch = error naming the read) → `cut_chunks(batch, paired)` (**a cut never falls between mates — the pair-in-one-chunk invariant; do not add a barrier either**) → `process_pair_chunk` → `build_pair_records` (flags, `RNEXT`/`PNEXT`/`TLEN`, `MC`, `ad`, per-mate `XA` entry *k* = tied pair *k*, secondaries per pair). The adapter bases are *not* aligned: the aligned mate is `seq[..len - ad]`, and `clipped_span` shifts a reverse hit's coordinates by `ad` so the soft clip lands at the right end. diff --git a/crates/eschalign-cli/src/align.rs b/crates/eschalign-cli/src/align.rs index 8e0f21b..3b4de0f 100644 --- a/crates/eschalign-cli/src/align.rs +++ b/crates/eschalign-cli/src/align.rs @@ -120,8 +120,8 @@ use sam::header::record::value::map::{self, Program, ReferenceSequence, program: use std::sync::Arc; use eschalign::{ - Aligner, Backend, CigarOp, Hit, MapOptions, Mode, PairMapping, PairOptions, Panel, ReadMapping, - ScoreMatrix, Scoring, sam as align_sam, + Aligner, Backend, CigarOp, Hit, MapOptions, Mode, PairMapping, PairOptions, PairOrientation, + Panel, ReadMapping, ScoreMatrix, Scoring, sam as align_sam, }; use rayon::prelude::*; @@ -160,6 +160,13 @@ pub struct AlignArgs { #[arg(long, value_name = "N")] pub max_insert: Option, + /// Paired input only: how the mates face each other in a proper pair. + /// `fr`: opposite strands, facing each other (mate 2 is scored on the + /// strand opposite mate 1); `rf`: opposite strands, facing away; `ff`: + /// the same strand, mate 1 upstream [default: fr] + #[arg(long, value_enum, value_name = "ORIENTATION")] + pub pair_orientation: Option, + /// Paired input only: mate 1's 3' adapter (its prefix is enough) /// [default: TruSeq Read 1, AGATCGGAAGAGCACACGTCTGAACTCCAGTCA] #[arg(long, value_name = "SEQ")] @@ -278,6 +285,23 @@ pub enum StrandArg { Both, } +#[derive(Clone, Copy, Debug, PartialEq, Eq, clap::ValueEnum)] +pub enum OrientationArg { + Fr, + Rf, + Ff, +} + +impl From for PairOrientation { + fn from(o: OrientationArg) -> Self { + match o { + OrientationArg::Fr => PairOrientation::Fr, + OrientationArg::Rf => PairOrientation::Rf, + OrientationArg::Ff => PairOrientation::Ff, + } + } +} + #[derive(Clone, Copy, Debug, PartialEq, Eq, clap::ValueEnum)] pub enum SortOrder { Unsorted, @@ -767,7 +791,7 @@ impl PairSource { if mate_stem(&name_a) != mate_stem(&name_b) { bail!( "read {} has no mate: the next record is {} (mates must be adjacent, \ - name-grouped)", + name-grouped: run `samtools collate` or `samtools sort -n` first)", show(&name_a), show(&name_b) ); @@ -885,6 +909,7 @@ fn read_pair_batches( #[derive(Clone)] struct PairCfg { max_insert: usize, + orientation: PairOrientation, /// `None` under `--no-adapter-trim`. adapters: Option, } @@ -959,6 +984,9 @@ fn pair_cfg_from(args: &AlignArgs) -> anyhow::Result { }; Ok(PairCfg { max_insert: args.max_insert.unwrap_or(500), + orientation: args + .pair_orientation + .map_or(PairOrientation::Fr, Into::into), adapters, }) } @@ -1166,6 +1194,8 @@ fn build_records( /// adapter length. const MC_TAG: Tag = Tag::new(b'M', b'C'); const AD_TAG: Tag = Tag::new(b'a', b'd'); +/// The mate's mapping quality, set beside `MC` when the mate is placed. +const MQ_TAG: Tag = Tag::new(b'M', b'Q'); /// How one mate was placed, however its pair resolved. struct MatePlacement<'a> { @@ -1184,15 +1214,22 @@ struct MatePlacement<'a> { unique: bool, } -/// `TLEN` of `hit` against its partner: the span from the leftmost aligned -/// base of either to the rightmost, positive for the leftmost mate (mate 1 -/// when both start together) and negative for the other. -fn template_length(hit: &Hit, partner: &Hit, first: bool) -> i32 { - let (a, b) = (&hit.alignment, &partner.alignment); - let len = a.ref_end.max(b.ref_end) - a.ref_start.min(b.ref_start); - let left = a.ref_start < b.ref_start || (a.ref_start == b.ref_start && first); - let len = i32::try_from(len).unwrap_or(i32::MAX); - if left { len } else { -len } +/// A placed mate's MAPQ: 60 when its reference is unique, else 0. +fn mate_mapq(unique: bool) -> i64 { + if unique { 60 } else { 0 } +} + +/// `TLEN` of `hit` against its partner, as `samtools fixmate` computes it: the +/// partner's 5' end minus this mate's, where a mate's 5' end is the start of a +/// forward alignment and the (exclusive) end of a reverse one. For a proper +/// pair that is the insert, positive on the forward-strand mate and negative on +/// the other; it is what keeps a `fixmate` round trip a no-op. +fn template_length(hit: &Hit, partner: &Hit) -> i32 { + let five = |h: &Hit| { + let a = &h.alignment; + i64::try_from(if h.reverse { a.ref_end } else { a.ref_start }).unwrap_or(i64::MAX) + }; + i32::try_from(five(partner) - five(hit)).unwrap_or(i32::MAX) } /// Rewrite a read pair's two input records into its output records: the two @@ -1217,7 +1254,7 @@ fn build_pair_records( for m in 0..2 { let rec = &mut recs[m]; kept[m] = rec.flags() & (Flags::QC_FAIL | Flags::DUPLICATE); - for t in OWNED_TAGS.iter().chain([&MC_TAG, &AD_TAG]) { + for t in OWNED_TAGS.iter().chain([&MC_TAG, &MQ_TAG, &AD_TAG]) { rec.data_mut().remove(t); } *rec.mate_reference_sequence_id_mut() = None; @@ -1293,14 +1330,18 @@ fn build_pair_records( clipped_span(oh, adn[o]).1, ); rec.data_mut().insert(MC_TAG, Value::String(mc.into())); + rec.data_mut() + .insert(MQ_TAG, int_value(mate_mapq(mates[o].unique))); if h.reference == oh.reference { - *rec.template_length_mut() = template_length(h, oh, m == 0); + *rec.template_length_mut() = template_length(h, oh); } } (Some(h), None) => { flags[m] |= Flags::MATE_UNMAPPED; *rec.mate_reference_sequence_id_mut() = Some(h.reference); *rec.mate_alignment_start_mut() = Some(position(h.alignment.ref_start)); + // The mate has no CIGAR; `samtools fixmate` writes `MC:Z:*`. + rec.data_mut().insert(MC_TAG, Value::String("*".into())); } (None, Some(oh)) => { // An unmapped mate sits at its partner's position, as samtools @@ -1312,6 +1353,15 @@ fn build_pair_records( *rec.alignment_start_mut() = Some(position(oh.alignment.ref_start)); *rec.mate_reference_sequence_id_mut() = Some(oh.reference); *rec.mate_alignment_start_mut() = Some(position(oh.alignment.ref_start)); + let mc = align_sam::cigar_string( + &oh.alignment.ops, + full[o], + clipped_span(oh, adn[o]).0, + clipped_span(oh, adn[o]).1, + ); + rec.data_mut().insert(MC_TAG, Value::String(mc.into())); + rec.data_mut() + .insert(MQ_TAG, int_value(mate_mapq(mates[o].unique))); } (None, None) => flags[m] |= Flags::MATE_UNMAPPED, } @@ -1359,7 +1409,7 @@ fn build_pair_records( if proper { f |= Flags::PROPERLY_SEGMENTED; if ph.reference == h.reference { - b = b.set_template_length(template_length(h, ph, m == 0)); + b = b.set_template_length(template_length(h, ph)); } } } @@ -1407,6 +1457,8 @@ struct Job<'a> { /// Read pairs seen, and of those how many resolved as proper pairs. pairs: AtomicU64, proper_pairs: AtomicU64, + /// Pairs whose mates both mapped, but to different references. + discordant_pairs: AtomicU64, } /// Whether `--max-read-len` keeps a read of `len` bases out of alignment. @@ -1605,6 +1657,7 @@ fn process_pair_chunk( max_ties: ctx.opts.max_ties, max_insert: cfg.max_insert, both_strands: ctx.opts.both_strands, + orientation: cfg.orientation, }; let mappings = ctx .aligner @@ -1638,6 +1691,11 @@ fn process_pair_chunk( ([place(0), place(1)], true) } PairMapping::Separate(ms) => { + if let (Some(a), Some(b)) = (ms[0].primary(), ms[1].primary()) + && a.reference != b.reference + { + ctx.discordant_pairs.fetch_add(1, Ordering::Relaxed); + } let place = |m: usize| { let me = &ms[m]; let others: Vec<&Hit> = me.hits.iter().skip(1).collect(); @@ -1836,6 +1894,7 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { } else { let given = [ ("--max-insert", args.max_insert.is_some()), + ("--pair-orientation", args.pair_orientation.is_some()), ("--adapter-r1", args.adapter_r1.is_some()), ("--adapter-r2", args.adapter_r2.is_some()), ("--min-adapter-overlap", args.min_adapter_overlap.is_some()), @@ -1860,8 +1919,11 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { }); if paired { info!( - "paired-end: {} mode, max insert {}, adapters {}", + "paired-end: {} mode, {:?} orientation, max insert {}, adapters {}", mode.name(), + pair_cfg + .as_ref() + .map_or(PairOrientation::Fr, |c| c.orientation), pair_cfg.as_ref().map_or(0, |c| c.max_insert), match pair_cfg.as_ref().and_then(|c| c.adapters.as_ref()) { Some(a) => format!("searched (min overlap {})", a.min_overlap), @@ -2039,6 +2101,7 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { pair: pair_cfg, pairs: AtomicU64::new(0), proper_pairs: AtomicU64::new(0), + discordant_pairs: AtomicU64::new(0), }; let busy = dispatch(&ctx, batches, res_tx, permit_rx); // Join every thread before deciding which error to report: a failed stage @@ -2123,6 +2186,14 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { 100.0 * proper as f64 / pairs.max(1) as f64, style::count(pairs - proper), ); + let discordant = ctx.discordant_pairs.load(Ordering::Relaxed); + if discordant > 0 { + info!( + "{} of the non-proper pairs are discordant: their mates map best to different \ + references, and each is reported at its own", + style::count(discordant) + ); + } } let wall = started.elapsed().as_secs_f64(); let mapped = summary.unique + summary.tied; diff --git a/crates/eschalign-cli/tests/align_paired.rs b/crates/eschalign-cli/tests/align_paired.rs index 2037c60..f95ced7 100644 --- a/crates/eschalign-cli/tests/align_paired.rs +++ b/crates/eschalign-cli/tests/align_paired.rs @@ -17,6 +17,7 @@ use sam::alignment::io::Write as _; use sam::alignment::record::Flags; use sam::alignment::record::data::field::Tag; use sam::alignment::record_buf::data::field::Value; +use std::io::Read; use std::path::{Path, PathBuf}; use std::process::{Command, Output}; @@ -459,7 +460,7 @@ fn align_single_input(dir: &Path, tag: &str, reads: &Path, extra: &[&str]) -> Ve } #[test] -fn pair_inputs_agree_fastq_interleaved_ubam() { +fn fastq_pair_matches_interleaved_and_ubam() { let dir = tempfile::tempdir().unwrap(); let (il, ubam) = other_forms(dir.path()); for extra in [&[][..], &["--strand", "both", "--secondary"][..]] { @@ -498,55 +499,117 @@ fn pair_inputs_agree_fastq_interleaved_ubam() { assert_eq!(lines(&suffixed), lines(&plain)); } +fn write_fastq(p: &Path, recs: &[(String, String, String)]) { + let t: String = recs + .iter() + .map(|(n, s, q)| format!("@{n}\n{s}\n+\n{q}\n")) + .collect(); + std::fs::write(p, t).unwrap(); +} + +/// Run eschalign expecting failure; its stderr. +fn fail(args: &[&std::ffi::OsStr]) -> String { + let o = eschalign(args); + assert!(!o.status.success(), "this run must fail"); + String::from_utf8_lossy(&o.stderr).into_owned() +} + #[test] -fn pair_orphan_is_error() { +fn mate_name_mismatch_is_an_error() { let dir = tempfile::tempdir().unwrap(); let r1 = fastq_records(&fx().join("pairs_R1.fq")); - let r2 = fastq_records(&fx().join("pairs_R2.fq")); - let write = |p: &Path, recs: &[(String, String, String)]| { - let t: String = recs - .iter() - .map(|(n, s, q)| format!("@{n}\n{s}\n+\n{q}\n")) - .collect(); - std::fs::write(p, t).unwrap(); - }; - let fail = |args: &[&std::ffi::OsStr]| -> String { - let o = eschalign(args); - assert!(!o.status.success(), "an orphan mate must fail the run"); - String::from_utf8_lossy(&o.stderr).into_owned() - }; + let mut r2 = fastq_records(&fx().join("pairs_R2.fq")); + let (full, sw) = (dir.path().join("full.fq"), dir.path().join("swapped.fq")); + write_fastq(&full, &r1); + r2.swap(0, 1); + write_fastq(&sw, &r2); let out = dir.path().join("o.bam"); - // Mate file one record short: the read left over is named. - let (short, full) = (dir.path().join("short.fq"), dir.path().join("full.fq")); - write(&short, &r2[..r2.len() - 1]); - write(&full, &r1); + // The first mismatch is named, with both files. let err = fail(&[ "-1".as_ref(), full.as_os_str(), "-2".as_ref(), - short.as_os_str(), + sw.as_os_str(), "-r".as_ref(), panel().as_os_str(), "-o".as_ref(), out.as_os_str(), ]); - assert!(err.contains("plain") && err.contains("no mate"), "{err}"); - // Names that do not match at the same record. - let mut swapped = r2.clone(); - swapped.swap(0, 1); - let sw = dir.path().join("swapped.fq"); - write(&sw, &swapped); + assert!(err.contains("paralog") && err.contains("pair_tie"), "{err}"); + assert!( + err.contains("full.fq") && err.contains("swapped.fq"), + "{err}" + ); + // Interleaved: the neighbour is not the mate. + let il = dir.path().join("il.fq"); + let recs: Vec<_> = r1 + .iter() + .zip(&r2) + .flat_map(|(a, b)| [a.clone(), b.clone()]) + .collect(); + write_fastq(&il, &recs); let err = fail(&[ - "-1".as_ref(), - full.as_os_str(), - "-2".as_ref(), - sw.as_os_str(), + il.as_os_str(), + "--interleaved".as_ref(), "-r".as_ref(), panel().as_os_str(), "-o".as_ref(), out.as_os_str(), ]); - assert!(err.contains("paralog") && err.contains("pair_tie"), "{err}"); + assert!(err.contains("paralog") && err.contains("no mate"), "{err}"); +} + +#[test] +fn non_adjacent_bam_mates_say_to_collate() { + let dir = tempfile::tempdir().unwrap(); + let (_, ubam) = other_forms(dir.path()); + let (header, mut recs) = read_bam(&ubam); + // Mate 1s first, then mate 2s: valid pairs, but not adjacent. + recs.sort_by_key(|r| r.flags().is_last_segment()); + let far = dir.path().join("far.bam"); + let mut w = bam::io::Writer::new(std::fs::File::create(&far).unwrap()); + w.write_header(&header).unwrap(); + for r in &recs { + w.write_alignment_record(&header, r).unwrap(); + } + w.try_finish().unwrap(); + let err = fail(&[ + far.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + dir.path().join("o.bam").as_os_str(), + ]); + assert!(err.contains("samtools collate"), "{err}"); +} + +#[test] +fn unequal_counts_is_an_error() { + let dir = tempfile::tempdir().unwrap(); + let r1 = fastq_records(&fx().join("pairs_R1.fq")); + let r2 = fastq_records(&fx().join("pairs_R2.fq")); + let (short, full) = (dir.path().join("short.fq"), dir.path().join("full.fq")); + write_fastq(&short, &r2[..r2.len() - 1]); + write_fastq(&full, &r1); + let out = dir.path().join("o.bam"); + let last = &r1.last().unwrap().0; + // Either file may be the short one; the read left over is named. + for (a, b) in [(&full, &short), (&short, &full)] { + let err = fail(&[ + "-1".as_ref(), + a.as_os_str(), + "-2".as_ref(), + b.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]); + assert!( + err.contains(last.as_str()) && err.contains("no mate"), + "{err}" + ); + } // Interleaved with a trailing orphan. let il = dir.path().join("il.fq"); let mut recs: Vec<_> = r1 @@ -559,7 +622,7 @@ fn pair_orphan_is_error() { "ACGTACGTACGTACGTACGT".into(), "I".repeat(20), )); - write(&il, &recs); + write_fastq(&il, &recs); let err = fail(&[ il.as_os_str(), "--interleaved".as_ref(), @@ -587,7 +650,10 @@ fn pair_orphan_is_error() { "-o".as_ref(), out.as_os_str(), ]); - assert!(err.contains("plain") && err.contains("no mate"), "{err}"); + assert!( + err.contains("discordant") && err.contains("no mate"), + "{err}" + ); } #[test] @@ -619,7 +685,7 @@ fn pair_ubam_tags_are_copied_on_both_mates() { } #[test] -fn pair_sort_coordinate_matches_samtools() { +fn sort_coordinate_matches_samtools_and_fixmate() { let dir = tempfile::tempdir().unwrap(); for extra in [&[][..], &["--secondary", "--strand", "both"][..]] { let (_, unsorted) = align_fixture(dir.path(), "u", extra); @@ -659,7 +725,7 @@ fn pair_sort_coordinate_matches_samtools() { // And samtools itself, when there is one. let Some(st) = samtools() else { eprintln!( - "[align_paired] pair_sort_coordinate_matches_samtools: no samtools; model only" + "[align_paired] sort_coordinate_matches_samtools_and_fixmate: no samtools; model only" ); continue; }; @@ -673,6 +739,46 @@ fn pair_sort_coordinate_matches_samtools() { assert!(o.status.success(), "{}", String::from_utf8_lossy(&o.stderr)); let (_, st_recs) = read_bam(&theirs); assert_eq!(lines(&sorted), lines(&st_recs), "samtools order {extra:?}"); + // fixmate round trip: name-sort our coordinate-sorted output, let + // samtools recompute every mate field, and compare. It must have + // nothing to change. + if extra.contains(&"--secondary") { + continue; // fixmate is defined for primary records + } + let sorted_path = dir.path().join("s.bam"); + let (byname, fixed) = (dir.path().join("n.bam"), dir.path().join("f.bam")); + let st = samtools().unwrap(); + let o = Command::new(&st) + .args(["sort", "-n", "--no-PG", "-o"]) + .arg(&byname) + .arg(&sorted_path) + .output() + .unwrap(); + assert!(o.status.success(), "{}", String::from_utf8_lossy(&o.stderr)); + let o = Command::new(&st) + .args(["fixmate", "-m", "--no-PG", "-O", "bam"]) + .arg(&byname) + .arg(&fixed) + .output() + .unwrap(); + assert!(o.status.success(), "{}", String::from_utf8_lossy(&o.stderr)); + let (_, fixed) = read_bam(&fixed); + let mate_fields = |r: &RecordBuf| { + ( + name(r), + r.flags().is_first_segment(), + flag(r) & (0x1 | 0x4 | 0x8 | 0x10 | 0x20 | 0x40 | 0x80), + r.mate_reference_sequence_id(), + r.mate_alignment_start(), + r.template_length(), + string_tag(r, b"MC"), + ) + }; + let mut ours: Vec<_> = sorted.iter().map(mate_fields).collect(); + let mut theirs: Vec<_> = fixed.iter().map(mate_fields).collect(); + ours.sort(); + theirs.sort(); + assert_eq!(ours, theirs, "samtools fixmate disagrees {extra:?}"); } } @@ -691,7 +797,7 @@ fn samtools() -> Option { /// record building, so the records must be identical. #[cfg(feature = "gpu")] #[test] -fn pair_gpu_matches_cpu() { +fn device_gpu_output_is_byte_identical() { let dir = tempfile::tempdir().unwrap(); for (k, extra) in [ &[][..], @@ -728,7 +834,7 @@ fn pair_gpu_matches_cpu() { stderr.contains("--device gpu cannot run"), "GPU run failed for a reason other than a missing device:\n{stderr}" ); - eprintln!("[align_paired] skipping pair_gpu_matches_cpu: {stderr}"); + eprintln!("[align_paired] skipping device_gpu_output_is_byte_identical: {stderr}"); return; } assert!(stderr.contains("on GPU"), "{stderr}"); @@ -803,3 +909,209 @@ fn pair_chunking_keeps_mates_together() { assert!(pair[0].flags().is_first_segment() && pair[1].flags().is_last_segment()); } } + +/// Mate 1 alone ties between the paralogs and mate 2 alone separates them +/// (`paralog` in the fixture: mate 1 sits where `P1` and `P2` are identical, +/// mate 2 covers the one base that differs). Together they are a unique pair. +#[test] +fn joint_score_breaks_tie_a_mate_cannot() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "joint", &[]); + let [a, b] = mates(&recs, "paralog"); + // Alone, each mate as the unpaired aligner sees it (both strands). + let alone = |reads: &str, tag: &str| { + let out = dir.path().join(format!("{tag}.bam")); + run_ok(&[ + fx().join(reads).as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + "--min-score".as_ref(), + "40".as_ref(), + "--mode".as_ref(), + "semiglobal".as_ref(), + "--strand".as_ref(), + "both".as_ref(), + ]); + let (_, r) = read_bam(&out); + r.into_iter().find(|r| name(r) == "paralog").unwrap() + }; + let m1 = alone("pairs_R1.fq", "m1"); + let m2 = alone("pairs_R2.fq", "m2"); + assert_eq!( + m1.mapping_quality().map(|q| q.get()), + Some(0), + "mate 1 ties" + ); + assert!(string_tag(&m1, b"XA").is_some()); + assert_eq!( + m2.mapping_quality().map(|q| q.get()), + Some(60), + "mate 2 decides" + ); + assert_eq!(ref_name(&h, &m2).as_deref(), Some("P1")); + // Jointly: both mates on P1, MAPQ 60, no XA, and each carries the other's + // MAPQ in MQ and CIGAR in MC. + for r in [a, b] { + assert_eq!(ref_name(&h, r).as_deref(), Some("P1")); + assert_eq!(r.mapping_quality().map(|q| q.get()), Some(60)); + assert!(string_tag(r, b"XA").is_none()); + assert_eq!(int_tag(r, b"MQ"), Some(60)); + assert_eq!(string_tag(r, b"MC").as_deref(), Some("60M")); + } +} + +#[test] +fn discordant_pair_flags_and_counts() { + let dir = tempfile::tempdir().unwrap(); + let out = dir.path().join("d.bam"); + let (r1, r2) = (fx().join("pairs_R1.fq"), fx().join("pairs_R2.fq")); + let stderr = run_ok(&[ + "-1".as_ref(), + r1.as_os_str(), + "-2".as_ref(), + r2.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + "--min-score".as_ref(), + "40".as_ref(), + ]); + let (h, recs) = read_bam(&out); + let [a, b] = mates(&recs, "discordant"); + // Each mate at its own best reference, never forced onto one. + assert_eq!(ref_name(&h, a).as_deref(), Some("U1")); + assert_eq!(ref_name(&h, b).as_deref(), Some("U2")); + assert_eq!(a.mate_reference_sequence_id(), b.reference_sequence_id()); + assert_eq!(b.mate_reference_sequence_id(), a.reference_sequence_id()); + assert_eq!(a.mate_alignment_start(), b.alignment_start()); + // Paired, mate 1 forward, mate 2 reverse, not proper, no TLEN. + assert_eq!(flag(a), 0x1 | 0x20 | 0x40); + assert_eq!(flag(b), 0x1 | 0x10 | 0x80); + assert_eq!((a.template_length(), b.template_length()), (0, 0)); + assert_eq!(string_tag(a, b"MC").as_deref(), Some("60M")); + assert_eq!(int_tag(a, b"MQ"), Some(60)); + // The run summary counts it (and only it: `wrong_orientation` is + // non-proper too, but on one reference). + assert!( + stderr.contains("1 of the non-proper pairs are discordant"), + "{stderr}" + ); + let [c, d] = mates(&recs, "wrong_orientation"); + assert!(!c.flags().is_properly_segmented()); + assert_eq!(c.reference_sequence_id(), d.reference_sequence_id()); +} + +/// `wrong_orientation` is two mates on the same strand: not a pair under the +/// default `fr`, a proper one under `ff`. +#[test] +fn pair_orientation_ff_makes_same_strand_pair_proper() { + let dir = tempfile::tempdir().unwrap(); + let (_, fr) = align_fixture(dir.path(), "fr", &[]); + let [a, _] = mates(&fr, "wrong_orientation"); + assert!(!a.flags().is_properly_segmented()); + let (h, ff) = align_fixture(dir.path(), "ff", &["--pair-orientation", "ff"]); + let [a, b] = mates(&ff, "wrong_orientation"); + assert_eq!(flag(a), 0x1 | 0x2 | 0x40); + assert_eq!(flag(b), 0x1 | 0x2 | 0x80); + assert_eq!(ref_name(&h, a), ref_name(&h, b)); + assert_eq!((a.template_length(), b.template_length()), (100, -100)); +} + +// --------------------------------------------------------------------------- +// Single-end input is unchanged +// --------------------------------------------------------------------------- + +/// The single-end path must write the same bytes it wrote before paired-end +/// support existed. The golden (`fixtures/unpaired_golden.tsv`) was measured +/// with the binary at the commit before `feat/paired-end`: for each flag +/// combination, the record count and an FNV-1a 64 digest of the *uncompressed* +/// BAM record stream (the header is skipped; its `@PG CL` carries the +/// temporary paths of the run). `ESCHALIGN_BLESS=1` rewrites it; do that only +/// when a single-end change is intended, and say so in the commit. +fn fixtures() -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures") +} + +const CONFIGS: &[(&str, &[&str])] = &[ + ("default", &[]), + ( + "semiglobal_both_secondary", + &["--mode", "semiglobal", "--strand", "both", "--secondary"], + ), + ("sort_coordinate", &["--sort", "coordinate"]), + ( + "min_score_max_ties", + &["--min-score", "60", "--max-ties", "1"], + ), + ( + "scoring_bwa", + &["--scoring", "1,-1,-2,-1", "--strand", "both"], + ), + ("max_read_len", &["--max-read-len", "60"]), +]; + +/// `(records, digest)` of a BAM's record stream, header skipped. +fn digest(path: &Path) -> (usize, u64) { + let mut raw = Vec::new(); + noodles_bgzf::io::Reader::new(std::fs::File::open(path).unwrap()) + .read_to_end(&mut raw) + .unwrap(); + let u32_at = |o: usize| u32::from_le_bytes(raw[o..o + 4].try_into().unwrap()) as usize; + assert_eq!(&raw[..4], b"BAM\x01"); + let mut off = 8 + u32_at(4); + let n_ref = u32_at(off); + off += 4; + for _ in 0..n_ref { + off += 4 + u32_at(off) + 4; + } + let body = &raw[off..]; + let (mut n, mut p) = (0, 0); + while p < body.len() { + p += 4 + u32::from_le_bytes(body[p..p + 4].try_into().unwrap()) as usize; + n += 1; + } + let mut h = 0xcbf2_9ce4_8422_2325u64; + for &b in body { + h ^= u64::from(b); + h = h.wrapping_mul(0x0000_0100_0000_01b3); + } + (n, h) +} + +#[test] +fn single_end_output_unchanged() { + let dir = tempfile::tempdir().unwrap(); + let mut got = String::new(); + for (label, extra) in CONFIGS { + let out = dir.path().join(format!("{label}.bam")); + let o = Command::new(env!("CARGO_BIN_EXE_eschalign")) + .arg(fixtures().join("trna_mappings_padded.bam")) + .arg("-r") + .arg(fixtures().join("trna_reference.fa")) + .arg("-o") + .arg(&out) + .args(*extra) + .output() + .expect("eschalign should launch"); + assert!( + o.status.success(), + "{label}: {}", + String::from_utf8_lossy(&o.stderr) + ); + let (n, h) = digest(&out); + got.push_str(&format!("{label}\t{n}\t{h:016x}\n")); + } + let golden = fixtures().join("unpaired_golden.tsv"); + // A test-only switch, not a runtime knob: any value blesses. + #[allow(clippy::disallowed_methods)] + let bless = std::env::var_os("ESCHALIGN_BLESS").is_some(); + if bless { + std::fs::write(&golden, &got).unwrap(); + return; + } + let want = std::fs::read_to_string(&golden).expect("fixtures/unpaired_golden.tsv"); + assert_eq!(got, want, "the single-end output changed"); +} diff --git a/crates/eschalign-cli/tests/align_unpaired.rs b/crates/eschalign-cli/tests/align_unpaired.rs deleted file mode 100644 index c171bf0..0000000 --- a/crates/eschalign-cli/tests/align_unpaired.rs +++ /dev/null @@ -1,100 +0,0 @@ -// SPDX-License-Identifier: MIT - -//! `unpaired_output_unchanged`: the unpaired path must write the same bytes it -//! wrote before paired-end support existed. -//! -//! The golden (`fixtures/unpaired_golden.tsv`) was measured with the binary at -//! the commit before `feat/paired-end`: for each flag combination, the record -//! count and an FNV-1a 64 digest of the *uncompressed* BAM record stream (the -//! header is skipped — its `@PG CL` carries the temporary paths of the run). -//! `ESCHALIGN_BLESS=1` rewrites it; do that only when an unpaired change is -//! intended, and say so in the commit. - -use std::io::Read; -use std::path::{Path, PathBuf}; -use std::process::Command; - -fn fixtures() -> PathBuf { - Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures") -} - -const CONFIGS: &[(&str, &[&str])] = &[ - ("default", &[]), - ( - "semiglobal_both_secondary", - &["--mode", "semiglobal", "--strand", "both", "--secondary"], - ), - ("sort_coordinate", &["--sort", "coordinate"]), - ( - "min_score_max_ties", - &["--min-score", "60", "--max-ties", "1"], - ), - ( - "scoring_bwa", - &["--scoring", "1,-1,-2,-1", "--strand", "both"], - ), - ("max_read_len", &["--max-read-len", "60"]), -]; - -/// `(records, digest)` of a BAM's record stream, header skipped. -fn digest(path: &Path) -> (usize, u64) { - let mut raw = Vec::new(); - noodles_bgzf::io::Reader::new(std::fs::File::open(path).unwrap()) - .read_to_end(&mut raw) - .unwrap(); - let u32_at = |o: usize| u32::from_le_bytes(raw[o..o + 4].try_into().unwrap()) as usize; - assert_eq!(&raw[..4], b"BAM\x01"); - let mut off = 8 + u32_at(4); - let n_ref = u32_at(off); - off += 4; - for _ in 0..n_ref { - off += 4 + u32_at(off) + 4; - } - let body = &raw[off..]; - let (mut n, mut p) = (0, 0); - while p < body.len() { - p += 4 + u32::from_le_bytes(body[p..p + 4].try_into().unwrap()) as usize; - n += 1; - } - let mut h = 0xcbf2_9ce4_8422_2325u64; - for &b in body { - h ^= u64::from(b); - h = h.wrapping_mul(0x0000_0100_0000_01b3); - } - (n, h) -} - -#[test] -fn unpaired_output_unchanged() { - let dir = tempfile::tempdir().unwrap(); - let mut got = String::new(); - for (label, extra) in CONFIGS { - let out = dir.path().join(format!("{label}.bam")); - let o = Command::new(env!("CARGO_BIN_EXE_eschalign")) - .arg(fixtures().join("trna_mappings_padded.bam")) - .arg("-r") - .arg(fixtures().join("trna_reference.fa")) - .arg("-o") - .arg(&out) - .args(*extra) - .output() - .expect("eschalign should launch"); - assert!( - o.status.success(), - "{label}: {}", - String::from_utf8_lossy(&o.stderr) - ); - let (n, h) = digest(&out); - got.push_str(&format!("{label}\t{n}\t{h:016x}\n")); - } - let golden = fixtures().join("unpaired_golden.tsv"); - // A test-only switch, not a runtime knob: any value blesses. - #[allow(clippy::disallowed_methods)] - let bless = std::env::var_os("ESCHALIGN_BLESS").is_some(); - if bless { - std::fs::write(&golden, &got).unwrap(); - return; - } - let want = std::fs::read_to_string(&golden).expect("fixtures/unpaired_golden.tsv"); - assert_eq!(got, want, "the unpaired output changed"); -} diff --git a/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py b/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py index debe94d..7f4039c 100644 --- a/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py +++ b/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py @@ -65,6 +65,8 @@ def add(name, r1, r2): add("one_unmapped", u2[20:80], rand(60)) # An ordinary pair. add("plain", u2[10:70], rc(u2[80:140])) +# Discordant: mate 1 comes from U1, mate 2 from U2 (each an unambiguous best). +add("discordant", u1[10:70], rc(u2[80:140])) def qual(n, shift): diff --git a/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq b/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq index 8cd2326..e97f2db 100644 --- a/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq +++ b/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq @@ -34,3 +34,7 @@ CTTGCGTATGGCACCTTTCCGGTATTGCTACAACCATCTGGGGGGCCGGATCTATTCTGG AAACGTCCGTCTTGCGTATGGCACCTTTCCGGTATTGCTACAACCATCTGGGGGGCCGGA + 58;>ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69 +@discordant +CACGGATGCAGTACAAAGTTGTATTCTATCAGTAGAGATATTCAACGTACCTCCCGTACG ++ +ADG69ADG69ADG69 diff --git a/crates/eschalign/src/lib.rs b/crates/eschalign/src/lib.rs index 17595d0..0da48d8 100644 --- a/crates/eschalign/src/lib.rs +++ b/crates/eschalign/src/lib.rs @@ -63,7 +63,8 @@ pub mod simd; pub use error::AlignError; pub use mapper::{ - Aligner, Hit, MapOptions, PairHit, PairMapping, PairOptions, ReadMapping, ScoreMatrix, + Aligner, Hit, MapOptions, PairHit, PairMapping, PairOptions, PairOrientation, ReadMapping, + ScoreMatrix, }; pub use panel::{Panel, Reference}; pub use scalar::{Alignment, CigarOp}; diff --git a/crates/eschalign/src/mapper.rs b/crates/eschalign/src/mapper.rs index 615edc7..0351149 100644 --- a/crates/eschalign/src/mapper.rs +++ b/crates/eschalign/src/mapper.rs @@ -24,7 +24,7 @@ use crate::scoring::{Mode, Scoring}; use crate::simd::{self, Backend, Kernel, Profile}; mod pair; -pub use pair::{PairHit, PairMapping, PairOptions}; +pub use pair::{PairHit, PairMapping, PairOptions, PairOrientation}; /// Per-read policy. The default: minimum score 0, every tie traced, forward /// strand only. diff --git a/crates/eschalign/src/mapper/pair.rs b/crates/eschalign/src/mapper/pair.rs index a218757..a859517 100644 --- a/crates/eschalign/src/mapper/pair.rs +++ b/crates/eschalign/src/mapper/pair.rs @@ -8,7 +8,9 @@ //! forward template strand mate 1 is read as sequenced and mate 2 is //! reverse-complemented (an `FR` pair, mate 1 upstream); with //! [`PairOptions::both_strands`] the reverse template strand is a candidate -//! too (mate 1 reverse-complemented, mate 2 as sequenced). A candidate's score +//! too (mate 1 reverse-complemented, mate 2 as sequenced). [`PairOrientation`] +//! picks the other layouts: `RF` keeps the strands but wants the mates facing +//! away from each other, `FF` puts both mates on the template strand. A candidate's score //! is `score(mate 1) + score(mate 2)`, each against the reference, and it is //! eligible only when **both** mates are at or above //! [`PairOptions::min_score`] — the sum alone does not let one good mate carry @@ -40,6 +42,52 @@ use super::*; +/// How the two mates of a proper pair face each other. +/// +/// Every candidate has a *template strand*: mate 1's strand. The orientation +/// fixes mate 2's strand relative to it, and the geometry a pair must have on +/// the reference (`--pair-orientation` on the CLI). +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)] +pub enum PairOrientation { + /// Mates on opposite strands facing each other (Illumina's usual + /// "innie"): the forward-strand mate starts before the reverse-strand + /// mate ends. + #[default] + Fr, + /// Mates on opposite strands facing away from each other ("outie"): the + /// reverse-strand mate starts at or before the forward-strand mate. + Rf, + /// Mates on the same strand, mate 1 upstream of mate 2 along it. + Ff, +} + +impl PairOrientation { + /// Whether mate `m` (0 or 1) is read reverse-complemented when mate 1 + /// follows the template strand `template_reverse`. + fn mate_reverse(self, m: usize, template_reverse: bool) -> bool { + match self { + PairOrientation::Fr | PairOrientation::Rf => template_reverse ^ (m == 1), + PairOrientation::Ff => template_reverse, + } + } + + /// The geometry check on two traced mates (`al[0]`, `al[1]`). + fn geometry_ok(self, al: &[Alignment; 2], template_reverse: bool) -> bool { + // Mate order along the template strand: mate 1 first on the forward + // template strand, mate 2 first on the reverse one. + let (first, second) = if template_reverse { + (&al[1], &al[0]) + } else { + (&al[0], &al[1]) + }; + match self { + PairOrientation::Fr => first.ref_start < second.ref_end, + PairOrientation::Rf => second.ref_start <= first.ref_start, + PairOrientation::Ff => first.ref_start <= second.ref_start, + } + } +} + /// Per-pair policy. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub struct PairOptions { @@ -52,6 +100,8 @@ pub struct PairOptions { pub max_insert: usize, /// Also consider the reverse template strand (mate 1 reversed). pub both_strands: bool, + /// How the mates must face each other. + pub orientation: PairOrientation, } impl Default for PairOptions { @@ -61,6 +111,7 @@ impl Default for PairOptions { max_ties: None, max_insert: 500, both_strands: false, + orientation: PairOrientation::Fr, } } } @@ -130,12 +181,6 @@ struct PairState<'a> { } impl PairState<'_> { - /// Mate `m`'s orientation in a candidate on template strand `reverse`: - /// mate 1 follows the template, mate 2 is its opposite. - fn mate_reverse(m: usize, template_reverse: bool) -> bool { - template_reverse ^ (m == 1) - } - fn open(&self) -> bool { self.found.is_empty() && self.next_tier < self.tiers.len() } @@ -194,7 +239,7 @@ impl Aligner { } for &t in strands { for m in 0..2 { - let rv = PairState::mate_reverse(m, t); + let rv = opts.orientation.mate_reverse(m, t); if st.rows[m][usize::from(rv)].is_empty() { let mut row = Vec::new(); self.score_or_take( @@ -205,8 +250,8 @@ impl Aligner { st.rows[m][usize::from(rv)] = row; } } - let r0 = &st.rows[0][usize::from(PairState::mate_reverse(0, t))]; - let r1 = &st.rows[1][usize::from(PairState::mate_reverse(1, t))]; + let r0 = &st.rows[0][usize::from(opts.orientation.mate_reverse(0, t))]; + let r1 = &st.rows[1][usize::from(opts.orientation.mate_reverse(1, t))]; for r in 0..n_refs { let (a, b) = (r0[r], r1[r]); if a >= opts.min_score && b >= opts.min_score { @@ -255,7 +300,7 @@ impl Aligner { for c in &st.cands[a..b] { let reference = &self.panel.get(c.reference).codes[..]; for m in 0..2 { - let rv = PairState::mate_reverse(m, c.reverse); + let rv = opts.orientation.mate_reverse(m, c.reverse); jobs.push((&st.codes[m][usize::from(rv)], reference)); } } @@ -309,7 +354,7 @@ impl Aligner { continue; } for (m, sub) in suboptimal.iter_mut().enumerate() { - let s = st.rows[m][usize::from(PairState::mate_reverse(m, t))][r]; + let s = st.rows[m][usize::from(opts.orientation.mate_reverse(m, t))][r]; *sub = Some(sub.map_or(s, |x: i32| x.max(s))); } } @@ -328,7 +373,7 @@ impl Aligner { let reference = self.panel.get(c.reference); let [a0, a1] = al; let mut mates = [a0, a1].into_iter().enumerate().map(|(m, alignment)| { - let rv = PairState::mate_reverse(m, c.reverse); + let rv = opts.orientation.mate_reverse(m, c.reverse); let letters: &[u8] = if rv { &rc[m] } else { st.letters[m] }; let (md, nm) = sam::md_nm( letters, @@ -393,14 +438,7 @@ impl Aligner { if al[0].is_empty() || al[1].is_empty() { return false; } - // The mate read as sequenced on the template strand is the "forward" - // one: mate 1 on the forward template strand, mate 2 on the reverse. - let (fwd, rev) = if template_reverse { - (&al[1], &al[0]) - } else { - (&al[0], &al[1]) - }; - if fwd.ref_start >= rev.ref_end { + if !opts.orientation.geometry_ok(al, template_reverse) { return false; } let lo = al[0].ref_start.min(al[1].ref_start); @@ -539,6 +577,62 @@ mod tests { } } + #[test] + fn orientation_rf_and_ff() { + let (al, a, _) = aligner(Mode::SemiGlobal); + let rf = PairOptions { + min_score: 30, + orientation: PairOrientation::Rf, + ..PairOptions::default() + }; + let ff = PairOptions { + orientation: PairOrientation::Ff, + ..rf + }; + // An FR pair (mate 1 forward upstream, mate 2 reverse downstream) is + // not an RF pair: it is the same alignments, so only the geometry + // decides, and an innie is not an outie. + let r1 = a[..30].to_vec(); + let r2 = alphabet::reverse_complement(&a[30..]); + assert!(matches!( + al.map_pairs( + &[[&r1, &r2]], + &PairOptions { + orientation: PairOrientation::Fr, + ..rf + } + )[0], + PairMapping::Proper { .. } + )); + assert!(matches!( + al.map_pairs(&[[&r1, &r2]], &rf)[0], + PairMapping::Separate(_) + )); + // An outie: mate 1 forward downstream, mate 2 reverse upstream. + let o1 = a[30..].to_vec(); + let o2 = alphabet::reverse_complement(&a[..30]); + match &al.map_pairs(&[[&o1, &o2]], &rf)[0] { + PairMapping::Proper { hits, .. } => { + assert!(!hits[0].mates[0].reverse && hits[0].mates[1].reverse); + } + other => panic!("not proper: {other:?}"), + } + // Same strand, mate 1 upstream. + let f1 = a[..30].to_vec(); + let f2 = a[30..].to_vec(); + match &al.map_pairs(&[[&f1, &f2]], &ff)[0] { + PairMapping::Proper { hits, .. } => { + assert!(!hits[0].mates[0].reverse && !hits[0].mates[1].reverse); + } + other => panic!("not proper: {other:?}"), + } + // Mate 2 upstream of mate 1 on the same strand is not FF. + assert!(matches!( + al.map_pairs(&[[&f2, &f1]], &ff)[0], + PairMapping::Separate(_) + )); + } + #[test] fn both_mates_must_clear_min_score() { let (al, a, _) = aligner(Mode::SemiGlobal); diff --git a/docs/align.md b/docs/align.md index 35a915b..765a0bb 100644 --- a/docs/align.md +++ b/docs/align.md @@ -58,6 +58,7 @@ reverse-strand record is refused. | `-t, --threads ` | Threads for parallel processing | | `-1, --r1 `, `-2, --r2 ` | Paired input: mate 1 and mate 2 files (plain or gzip). Mates are matched by name, `/1` and `/2` stripped — see [Paired-end](#paired-end) | | `--interleaved` | `READS` is an interleaved FASTQ/FASTA (mate 1, mate 2, …). A BAM's pairing comes from its flags and needs no option | +| `--pair-orientation ` | Paired only: how the mates face each other (default `fr`: opposite strands facing each other; `rf`: opposite strands facing away; `ff`: same strand, mate 1 upstream) | | `--max-insert ` | Paired only: longest insert a proper pair may have (default `500`) | | `--adapter-r1 `, `--adapter-r2 ` | Paired only: 3' adapter of each mate (default TruSeq Read 1 / Read 2) | | `--min-adapter-overlap ` | Paired only: shortest 3' match counted as adapter (default `5`) | @@ -313,8 +314,9 @@ Records are in input order with the two mates adjacent, mate 1 then mate 2; |-------|--------------| | `FLAG` | 0x1; 0x2 only for a proper pair; 0x10 / 0x20 for each mate's / its partner's strand; 0x40 / 0x80; 0x4 / 0x8 for an unmapped mate / partner; QC-fail and duplicate (0x200, 0x400) from the input | | `RNEXT` / `PNEXT` | the partner's reference and position (`=` is the same reference); for an unmapped partner, the mate's own | -| `TLEN` | the true insert — leftmost aligned base to rightmost, over both mates — positive on the leftmost mate (mate 1 if they start together), negative on the other. Set for any pair with both mates on one reference, proper or not | -| `MC:Z` | the partner's CIGAR (absent when the partner is unmapped) | +| `TLEN` | as `samtools fixmate` computes it: the partner's 5' end minus this mate's (a forward mate's 5' end is its start, a reverse mate's its end), so for a proper pair the insert, positive on the forward-strand mate and negative on the other. Set for any pair with both mates on one reference, proper or not | +| `MC:Z` | the partner's CIGAR (`*` when the partner is unmapped) | +| `MQ:i` | the partner's MAPQ (60 unique, 0 tied), when the partner is mapped | | `AS:i` | that mate's own score (the pair score is the two summed) | | `XS:i` | that mate's best score over the candidates outside the pair's tie set | | `XA:Z` | `ref,±pos,CIGAR,NM;` — see below |