[BenjaminDEMAILLE]

Part of #354 (finding 2).

Moved from BenjaminDEMAILLE#1. The base, perf/hotpaths, is a copy of #346's branch at aeba641, pushed here so this PR can target it. If #346 gets new commits, the copy needs them too. Once #346 merges, retarget this PR to main and delete perf/hotpaths.

Stacked on #346 (perf/hotpaths, tip aeba641). Only the one commit is this PR. It does not overlap with the PR for finding 3: the two merge without a conflict, and the merge is measured there.

Problem

A pair with one aligned mate kept one location and skipped the multimapper filters. A mate in a 3-copy repeat, with its other mate unaligned, was reported as uniquely mapped at one copy:

3-copy mate, other mate unaligned	STAR	perf/hotpaths	This branch
Alignment records	3, NH:3	1	3, NH:3
SAM records, no --outSAMunmapped	3	2	3
Log: unique / multiple loci	0 / 1	1 / 0	0 / 1

It also wrote an unmapped-mate record that was not asked for.

What STAR does

Observed output of Homebrew STAR 2.7.11b; source not read.

• A 3-copy mate gets NH:3, MAPQ 1, flags 73/329/329, in the order 5051, 15651, 10351, and no unmapped-mate record unless --outSAMunmapped Within is set.
• With Within: one unmapped-mate record after the alignments, tied to the primary (flags 133 / 69 / 165).
• With Within KeepPairs: one after each single-mate alignment (389 / 325 / 421 on secondaries, with the AS and nM of its transcript).
• A pair and single-mate alignments of one read share a score range: with --outFilterMultimapScoreRange 30, the pair's records say NH:3 HI:1 and the singles follow.
• NH counts alignments past --outSAMmultNmax.
• Mates that align but cannot pair are two alignments, NH:2, each a single.
• --outFlagOR and --outFlagAND apply to the mapped record only.
• A mate held by the mismatch limit is "too short", since extensions stop at that limit.
• --outReadsUnmapped Fastx writes both mates of a half-mapped pair.

Change

• src/align/read_align.rs: PairedAlignmentResult gets score, output_rank, g_length and mates. One Vec<PairedAlignmentResult> in window order goes through one dedup, score range, order, filter_paired_transcripts and too-many-loci check. The half-mapped fallback is gone, which is why the file shrinks.
• src/io/sam.rs: build_paired_records takes the whole list. build_half_mapped_records writes the mapped mate first and returns the unmapped mate separately, so the caller places it for Within and KeepPairs.
• src/lib.rs: record_pair_alignment_stats and the output branches count results.len(). The half-mapped special cases are removed.
• src/solo/mod.rs, src/wasp/mod.rs, src/quant/mod.rs, src/quant/splice_status.rs: take the list or its length. More than one alignment is a multimapper; no pair is unmapped.

Measured against STAR

Fields 1-11 plus NH, HI, AS, nM.

Fixtures: the audit's h1.fq/h2.fq on perf/hotpaths had 2 records against STAR's 3, unique 1 against 0, multiple loci 0 against 1. Here they are identical, in the same order. A crafted 11-pair fixture under 19 option sets is identical in the same order, with identical Log.final.out counts. On that fixture perf/hotpaths had 20 records against STAR's 22, unique 10 against 4, and multiple loci 0 against 6.

Fuzz: 4000 repeat-rich pairs. Reads differing from STAR, and records:

Options	Differing, perf/hotpaths	Differing, this branch	Records STAR / perf/hotpaths / this branch
defaults	523	0	7182 / 7086 / 7182
score and match filters 0.3	1311	0	12510 / 12490 / 12510
same, --outFilterMultimapScoreRange 30	1588	0	13686 / 13522 / 13686
same filters, --outSAMunmapped Within	754	0	14024 / 12698 / 14024
filters 0, Within KeepPairs	887	0	16130 / 13274 / 16130
chimeric detection on	1312	0
--outFilterType BySJout	1309	0
--outSAMmultNmax 2	1770	0	8878 / 9692 / 8878
--peOverlapNbasesMin	524	1

The one --peOverlapNbasesMin read (r2383) has flag 163 where STAR has 161, the proper-pair bit. perf/hotpaths has the same difference.

Unmapped.out.mate* contains the same reads as STAR; the only difference is STAR's name comment ( 0:N: 10), already a row in #354. BAM Unsorted passes samtools quickcheck with 22 records, and I did not compare its fields.

Deliberate divergences

None, so no DIVERGENCE.md entry and no sign-off needed for this. The differences listed next are not new choices, except the shortcut: one, which I flag for your review.

Not the same as STAR, or not fixed here

• Log.final.out says Mismatch rate 0.00% where STAR says nan% when no read is unique. Older than this PR.
• --outMultimapperOrder Random differs from STAR, also on perf/hotpaths: another random number sequence.
• Single-end NH under --outSAMmultNmax still reports the output count: NH:2 where STAR has NH:3. The pair path now follows STAR; the single-end builders were out of scope. It belongs in Tracking: fixes from the whole-repo audit (main at 32fbbc1) #354.
• A read with any single-mate alignment gets no WASP (vW), sp, or solo gene tags. Marked shortcut: in src/lib.rs. I did not compare these with STAR.
• Single-mate alignments are not added to bedGraph signal. Not compared with STAR.
• --outFilterMismatchNoverReadLmax is refused by this binary, though STAR accepts it. Seen while testing, out of scope.

Tests

tests/half_mapped_multimap.rs, all seven failing on the old src/:

• mate_in_a_repeat_is_a_multimapper_on_every_copy
• unmapped_mate_record_only_on_request
• single_mates_are_held_to_the_loci_and_mismatch_limits
• pair_and_single_mates_share_the_score_range
• nh_counts_alignments_past_the_output_limit
• mates_that_cannot_pair_are_two_alignments
• repeat_mate_counts_as_multimapping_in_gene_counts: its expected value is not from STAR, see below.

Adapted in src/io/sam.rs: test_build_half_mapped_flags, test_build_half_mapped_rnext_pnext, test_build_half_mapped_mate_order.

cargo test --release (826 passed), cargo clippy --all-targets and cargo fmt --check are clean on macOS.

Not checked

• test/parity_suite.sh was not run.
• --quantMode and STARsolo cannot be compared here: on macOS STAR fails with could not open input file /geneInfo.tab. So ReadsPerGene.out.tab and solo counts on half-mapped reads are not checked against STAR.
• No real dataset: all numbers are from synthetic reads.
• Speed was not measured.

🤖 Generated with Claude Code