[BenjaminDEMAILLE]

Part of #354 (finding 19).

Stacked on #365 (chore/dead-code), itself stacked on #346. Merge after those.

Problem

src/lib.rs has four read loops (single-end, paired-end, STARsolo, STARsolo paired) of 450 to 760 lines each, which are near copies. A rule added to one did not reach the others.

Change

Plain functions in src/lib.rs, before the loops. No trait, no framework, one shared stage.

Shared now	Used by	Replaces
OutputRules::build_records: which records a read may produce (--outSAMtype None, --outSAMunmapped, multimap limit)	all 4	4 copies of the gate
spawn_decoder, read_limit, output_read_name	all 4	4 copies
apply_rg_per_file_batch, now generic over the result	all 4	4 copies
record_single_stats	SE, solo	2 copies
PairOutcome, record_pair_stats, build_aligned_pair_records (with a hook for what only the loop knows: WASP and splicing tags, or gene tags)	PE, solo paired	2 copies
run_writer_stage: chimeric, held reads, signal, SAM/BAM, transcriptome, unmapped, BySJout second stage, signal tracks, trailer, flushes	SE, PE	2 copies
open_chimeric_writer, load_wasp_context	SE, PE	2 copies
SoloReadProduct, SoloRules, consume_solo_products, solo_tag_out, add_gene_tags_in_output_order, add_solo_read_tags	solo, solo paired	2 copies

Function	Lines before	Lines after
align_reads_single_end	635	444
align_reads_solo	450	277
align_reads_solo_pe	522	297
align_reads_paired_end	760	518
Shared functions	0	640

src/lib.rs goes from 4319 to 4190 lines (43 of them a new test), so the gain is modest in total lines: the loops lose 831, the shared stage gains 640.

Still per loop. The solo loops keep their consumer on the dispatcher thread: moving them to a writer thread would change threading, so they do not use run_writer_stage. The clip-to-empty early returns also remain in each loop, because they differ (below); they read rules.output_unmapped.

Hot path. No new per-read allocation: build_records takes closures, and the gate pushes into the same buffer. PairOutcome::both_mapped is the Vec that was already built per pair, now of &PairedAlignment, and the PE loop no longer builds a second one for quant. AlignmentBatchResults and SoloReadProduct keep their fields, so batch_size_for gives the same batch sizes.

One difference, on stderr only: PE now logs Chimeric junction output complete at info level, as SE did. Seen in a diff of the two stderr logs of a chimeric PE run; it was the only line that differs in that run.

Differences between the loops that look like bugs (kept, not fixed here)

Read from the code; none of them changes in this PR.

1. Unmapped SAM records keep the full read in SE and PE (build_unmapped_record(&read.sequence…)), but use the clipped sequence and quality in solo and solo paired.
2. A read clipped to nothing gets an unmapped record in SE and PE; in solo it gets none, even with --outSAMunmapped Within.
3. Solo paired has no check for a mate clipped to nothing; SE, PE and solo do.
4. Solo ignores the adapter clip (clip_mate) unless CellRanger4; solo paired has no CellRanger4 path and uses fixed clips only.
5. The solo loops do not support quant, unmapped FASTQ, the transcriptome writer, WASP, the splicing tag, chimeric, BySJout or bedGraph (finding 4 of Tracking: fixes from the whole-repo audit (main at 32fbbc1) #354).
6. Unmapped FASTQ: SE writes the full read, PE writes the clipped mates.
7. PE's transcriptome projection counts results.len() towards the multimap limit, where signal_contrib and SAM output count both_mapped.
8. Under WithinBAM, PE records a chimeric read without checking chim_segment_min, SE checks it. Probably harmless, as chimeras only exist when it is set.
9. Solo paired counts SJ-feature junctions only for one both-mapped pair; half-mapped pairs add none.
10. The read index for solo tags is cast to u32, in both solo loops: it wraps past 4 294 967 295 reads.

Measured

macOS, aarch64, 16 cores. Baseline: clean build of aeba641 (before the dead-code PR). Also compared with the dead-code PR's bcddd31. Comparison: this branch.

Output. 84 scenarios, each at --runThreadN 1, 4 and 8, every binary under the same argv[0] and relative prefix so that command lines in headers match. md5 of every output file, raw BAM and SAM included, plus decoded BAM records and Log.final.out without its clock lines. Log.out, Log.progress.out and stderr are not in the md5s.

Set	Scenarios	Rows compared	Rows differing from baseline
Generated fixture (20 000 reads) and tiny fixtures	75	2367	4 (below)
Earlier fixtures, including rejected inputs	9	132	0

The 4 rows are Solo.out/{Gene,GeneFull}/CellReads.stats at 4 and 8 threads in one scenario. This is existing behaviour: two runs of the baseline differ from each other in that file, and baseline, baseline and this branch have the same md5 once sorted (64136fbd), so only the row order changes between runs. At 1 thread it is identical. Not touched here.

Scenarios cover all four loops: SAM, BAM unsorted and sorted, --outStd, Within, TranscriptomeSAM, GeneCounts, Fastx, chimeric junctions, WithinBAM, BySJout, fixed clip, clip to empty, adapter, CellRanger4, bedGraph, WASP, splicing tag, per-file RG, NoQS, None, 2-pass, read limits, multimap limits 0 and 1, half-mapped pairs, STARsolo Gene, GeneFull, SJ, Velocyto, EM, CB_samTagOut, a third barcode file and --soloBarcodeMate 1.

The fixture was checked to reach these paths on the baseline: SE 15 965 unique, 1614 multi, 1417 too short, 1003 other, of 20 000 reads; a second fixture gives 2151 (SE) and 2152 (PE) chimeric reads; relaxed --outFilter*OverLread gives 92 records of half-mapped pairs in PE and in solo paired.

Scenarios the program refuses: sp_tag_m (an invalid parameter combination of mine) and 5 of the 9 earlier scenarios (input files the reader rejects). Same exit code 1 in all binaries; the error text is the same for the 3 I compared, which exercises the decoder's error path.

Decoded records and every non-header file are also identical across 1, 4 and 8 threads, apart from the same CellReads.stats rows.

Speed. --runThreadN 8. The machine was not idle: unrelated bwa-mem2 and salmon jobs were running, so wall time is noisy (the same binary varies between 4 and 16 s). Instructions retired and CPU time are the comparison; wall time is given as the minimum of 5. From /usr/bin/time -l, binaries alternated.

Case	Input	Instructions, dead-code PR	this PR	Change	CPU s (median of 5)	this PR	Wall min, dead-code PR	this PR
SE, largest fixture	B200.fq, 801 000 reads, 176 MB	93.32 G	93.24 G	−0.08 %	10.29	10.55	4.92	4.21
SE	800 000 reads	344.33 G	344.26 G	−0.02 %	28.67	28.75	6.30	6.70
PE	800 000 pairs	1400.2 G	1400.0 G	−0.01 %	106.69	105.30	20.35	30.27
Solo	800 000 reads	350.08 G	350.23 G	+0.04 %	28.52	28.57	9.40	9.25
Solo paired	800 000 pairs	1407.1 G	1406.8 G	−0.02 %	106.33	107.52	30.65	29.73

Instructions are the median of 3 runs for SE and solo, 1 run for PE and solo paired. No measurable change in instructions; the CPU and wall differences are within the noise of this loaded machine. Median peak RSS, dead-code PR then this PR: 0.069 and 0.069 GB (B200.fq), 0.053 and 0.056 (SE), 0.063 and 0.064 (PE), 0.098 and 0.090 (solo), 0.109 and 0.119 (solo paired). Thread timing makes it vary; the solo figures of the baseline are 0.104 and 0.105.

Tests

cargo test --release: 816 tests pass, tests/determinism.rs and tests/phase9_threading.rs included. New: output_rules_tests::records_follow_the_output_rules for the gate. cargo clippy --all-targets has no warnings; cargo fmt --check passes.

Not checked

• I did not run test/parity_suite.sh, and did not compare against STAR here; the claim is "same output as before".
• Real sequencing data, gzip input, Linux.
• Log.out content.
• Wall time on an idle machine.
• Conflict with Write SJ.out.tab when empty, report a failed final SAM write #355: git merge-tree reports conflicts in CHANGELOG.md, src/io/sam.rs and src/lib.rs, as the dead-code PR describes. I did not carry out the resolution.

🤖 Generated with Claude Code