Skip to content

feat: pair orientation, MQ, discordant count, fixmate-exact TLEN #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jayhesselberth wants to merge 1 commit into
base: main
Choose a base branch
from
Open

feat: pair orientation, MQ, discordant count, fixmate-exact TLEN #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -22,15 +26,15 @@
- 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

- 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

Expand Down
4 changes: 2 additions & 2 deletions CLAUDE.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand Down
103 changes: 87 additions & 16 deletions crates/eschalign-cli/src/align.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::*;

Expand Down Expand Up @@ -160,6 +160,13 @@ pub struct AlignArgs {
#[arg(long, value_name = "N")]
pub max_insert: Option<usize>,

/// 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<OrientationArg>,

/// Paired input only: mate 1's 3' adapter (its prefix is enough)
/// [default: TruSeq Read 1, AGATCGGAAGAGCACACGTCTGAACTCCAGTCA]
#[arg(long, value_name = "SEQ")]
Expand Down Expand Up @@ -278,6 +285,23 @@ pub enum StrandArg {
Both,
}

#[derive(Clone, Copy, Debug, PartialEq, Eq, clap::ValueEnum)]
pub enum OrientationArg {
Fr,
Rf,
Ff,
}

impl From<OrientationArg> 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,
Expand Down Expand Up @@ -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)
);
Expand Down Expand Up @@ -885,6 +909,7 @@ fn read_pair_batches(
#[derive(Clone)]
struct PairCfg {
max_insert: usize,
orientation: PairOrientation,
/// `None` under `--no-adapter-trim`.
adapters: Option<AdapterSet>,
}
Expand Down Expand Up @@ -959,6 +984,9 @@ fn pair_cfg_from(args: &AlignArgs) -> anyhow::Result<PairCfg> {
};
Ok(PairCfg {
max_insert: args.max_insert.unwrap_or(500),
orientation: args
.pair_orientation
.map_or(PairOrientation::Fr, Into::into),
adapters,
})
}
Expand Down Expand Up @@ -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> {
Expand All @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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,
}
Expand Down Expand Up @@ -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));
}
}
}
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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()),
Expand All @@ -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),
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
Loading
Loading