fix(sj): outFilterType BySJout and SJ.out.tab filtering as STAR (outputSJ, two-stage held reads)

[BenjaminDEMAILLE]

Stacked on #315.

Finding

With --outFilterType BySJout the genomic output and SJ.out.tab differed from STAR 2.7.11b on the yeast 200k pairs (e.g. 802 mates STAR keeps that rustar dropped, 502 diffs). SJ.out.tab also differed with Normal filtering (missing and extra rows, wrong overhang and annotation columns).

Cause (STAR source)

• outputSJ.cpp: the count filter casts thresholds to uint (-1 never passes); the overhang test is per junction; outSJfilterIntronMaxVsReadN is indexed by the read count and compared with the gap (end - start + 1); the distance filter works on the sorted list of junctions that passed the count filter (annotated ones included as neighbours, exempt themselves), donor and acceptor. rustar used fixed indices, gap off by one, a donor/acceptor test on intron end vs next start, and -1 as "no limit".
• ReadAlign_outputAlignments.cpp (outFilterBySJout), ReadAlignChunk_processChunks.cpp: reads with an unannotated junction are held, the others are output at once. After the first stage the novel junctions (count + distance filter over the junctions of ALL reads, chunkOutSJ1) are fixed, and the held reads are mapped again (stitch rejects other unannotated junctions) and output after all other reads. The final SJ.out.tab (outFilterBySJoutStage==2) uses the junctions of the output reads only, with no distance filter. rustar instead dropped reads whose primary alignment had a non-surviving junction.
• ReadAlign_outputTranscriptSJ.cpp: overhang is the shorter adjacent exon, and indels and insertions end an exon; annotation is by sjdb coordinates whatever the strand (non-canonical annotated junctions, e.g. HAC1, were missed).

Fix

• junction/sj_output.rs: surviving_junctions ported from outputSJ; first_stage stats and bysjout_novel_junctions; no distance filter in the last BySJout stage.
• lib.rs (SE and PE): the per-batch aligner is now process(base, reads, stage2); first stage holds reads and records their junctions in the first-stage stats; the writer then maps held reads again with the novel set (same code path, so stats, SJ, WASP, transcriptome draws follow STAR order). The temp-file buffering and the post-hoc drop are gone (also saves disk and RAM).
• Overhang by exon (indels split exons) and sjdb annotation by coordinates or stitcher flag.

Results (yeast 200k pairs, sidx, vs STAR 2.7.11b, --runThreadN 1)

• BySJout genomic: same 332790 ties 4052 diff 0 only_rustar 0 only_star 0 (was diff 502, only_star 802); read order equal to STAR (held reads last).
• BySJout and Normal SJ.out.tab: identical to STAR (0 diff lines; before 40 and 21).
• BySJout Aligned.toTranscriptome.out.bam: identical records to STAR.
• Log.final.out mapped counts equal STAR. Normal Aligned.out.sam byte-identical to the starting build.

Tests

Unit tests for the intron-length filter, the uint cast of -1 and the distance filter (annotated neighbours, last stage); BySJout integration test now checks the second-stage log. cargo test --release, clippy -D warnings, fmt pass.

Checked again by the main session

200k yeast pairs, annotated index, rustar on 8 threads against STAR 2.7.11b --runThreadN 1:

	result
BySJout, genomic SAM	0 mates differ (0 only-rustar, 0 only-STAR)
BySJout, SJ.out.tab	identical (265 lines)
BySJout, transcriptome BAM	274,606 records identical, in order
Normal, genomic SAM	byte-identical to #311
Normal, SJ.out.tab	identical to STAR (it was not before)

Not covered: the two-pass novel-junction filter (--twopassMode Basic) still has its own filter.

Part of #298.

🤖 Generated with Claude Code