feat: STARlong binary + long-read window-coverage filter

[BenjaminDEMAILLE]

STARlong: the rustar-aligner-long binary, the long-read window-coverage filter, and two places where the long-read path did not match STARlong.

Correction to the previous revision of this description

I wrote: "Not run: validation against real STARlong output. [...] I do not have a STARlong binary and a long-read dataset available."

The binary claim was wrong. STARlong 2.7.11b builds from the same source tree as STAR and was available the whole time. The differential is below, and running it found two defects that reading the code had not.

What changed

rustar-aligner-long, mirroring STAR's STARlong: identical CLI, forcing the long-read stitcher via Parameters::long_read. Native STAR selects this by binary rather than by flag, and this matches that.

window_read_coverage, STAR's WC accumulation from ReadAlign_stitchPieces.cpp:208-220. The --winReadCoverageRelativeMin / --winReadCoverageBasesMin cluster filter is gated on long_read, since standard STAR does not filter clusters by seed coverage. --winReadCoverageBasesMin is new; --winReadCoverageRelativeMin already existed with default 0.5 and becomes effective.

The per-motif stitch mismatch cap no longer applies to long reads. stitchAlignToTranscript.cpp:308-315 gates the junction differently in the two binaries:

#ifdef COMPILE_FOR_LONG_READS
    if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal )
#else
    if ( (trA->nMM + nMM) <= outFilterMismatchNmaxTotal
         && ( jCan<0 || (jCan<7 && nMM <= P.alignSJstitchMismatchNmax[(jCan+1)/2]) ) )
#endif

The alignSJstitchMismatchNmax term is compiled out of STARlong. Applying it anyway did not cost one junction, it cost whole reads: the chaining DP split into a chain either side of the rejected junction and neither half cleared --outFilterScoreMinOverLread. Traced on sim_12_Q0075, where the killed junction is the annotated 35-base intron Mito 24871-24905, rejected as non-canonical because of a single mismatch beside it.

A traceback edge that fails to re-stitch no longer discards its window. Phase A scores an edge with one function and the traceback re-stitches with stitch_align_to_transcript, so the two can disagree. This returned an empty transcript vector on a disagreement. STAR does not: stitchWindowSeeds (ReadAlign_stitchWindowSeeds.cpp:130-138) adds whatever stitchAlignToTranscript returns unconditionally, with the failure return commented out and marked "this should not happen".

One faithful defect, kept

STAR initialises rLast = 0 before the coverage loop (ReadAlign_stitchPieces.cpp:208) and then tests r1 + L1 > rLast + 1, so a seed starting at read position 0 is counted one base short. That is reproduced rather than corrected: the filter thresholds are calibrated against STAR's own numbers, so fixing the count in isolation would change which windows survive relative to the oracle.

Verification against STARlong 2.7.11b

The bundled tiers are 150 bp Illumina reads and never reach this path: a short read produces one window with near-total coverage, so the coverage filter can discard nothing. This PR adds test/simulate_long_reads.py, which emits the full spliced sequence of each multi-exon transcript as one read, with substitution errors from a seeded RNG so the FASTQ is byte-reproducible. 500 reads, mean 2.1 kb, longest 5.6 kb, every read multi-exon.

Both tools run single-threaded with STAR's documented long-read parameter set. STAR maps 430 of 500.

	primaries	exact records	same locus	AS lower than STAR	unmapped by rustar
before	419	360	405	32	13
after	432	372	423	9	1

Exact records compare FLAG, RNAME, POS, MAPQ, CIGAR, NH and AS.

python3 test/simulate_long_reads.py "$DATA/reference/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa" \
  "$DATA/reference/Saccharomyces_cerevisiae.R64-1-1.110.gtf" "$DATA/reads/sim_long_500.fastq" 500 1000
LONG="--outFilterMismatchNmax 1000 --seedSearchStartLmax 30 --seedPerReadNmax 100000 \
      --seedPerWindowNmax 100 --alignTranscriptsPerWindowNmax 10000"
STARlong --genomeDir "$STAR_IDX"   --readFilesIn "$DATA/reads/sim_long_500.fastq" --runThreadN 1 \
  --outSAMtype SAM --outFileNamePrefix star/   $LONG --alignTranscriptsPerReadNmax 100000
./target/release/rustar-aligner-long --genomeDir "$RUSTAR_IDX" \
  --readFilesIn "$DATA/reads/sim_long_500.fastq" --runThreadN 1 \
  --outSAMtype SAM --outFileNamePrefix rustar/ $LONG

Two notes on that command, both found by running it.

--alignTranscriptsPerWindowNmax 10000 with the default --alignTranscriptsPerReadNmax maps zero reads, in STAR as well: STAR's headroom break is trNtotal + alignTranscriptsPerWindowNmax >= alignTranscriptsPerReadNmax, which is already true on the first window, and STARlong logs 435 such warnings and writes an empty SAM. STAR's own recommended set raises both, so it does not hit this. It is not a rustar-aligner defect, and I checked before assuming it was.

STARlong at --runThreadN 8 wrote an empty SAM here while exiting 0. Single-threaded it is fine. That is the STAR binary or my harness, not rustar-aligner, and both sides were run single-threaded so the comparison is like for like.

What is left

One read STAR maps and this does not, sim_403_YER179W: STAR emits 132M92N873M, and the chain built here covers only 240 of 1005 bases. The last two exons are at the right locus (V:548432, V:548645), but the window holds no seed for read positions 240 onward; those seeds went to unrelated clusters. Disabling the coverage filter changes nothing, so it is seed distribution rather than the stitcher or the filter. That is the same storeAligns piece-ordering work already named in #146, not something this PR can close.

Three reads this maps and STAR does not. STAR reports them unmapped: other, which is a failure to find a window rather than a rejection, and they land on the mitochondrial transcript they were simulated from with the exon structure STAR itself produces for other reads of that transcript. Recorded in DIVERGENCE.md §2.2 rather than corrected.

Nine reads where AS is lower than STAR's, and all nine belong to #146, not here: seven are the one-base annotated intron on YOR239W written as 1D where STAR writes 1N, and two are a terminal micro-exon soft-clipped instead of spliced. Both classes are fixed in #146. On a local merge of the two branches, exact records rise from 372 to 402; that merge is mine and not either PR, so the number is indicative rather than a claim about this diff.

Verification

Unit tests in src/align/read_align.rs:

• window_read_coverage_single_seed: a 10-base seed at read position 0 counts 9, pinning the undercount above
• window_read_coverage_seed_not_at_zero: the same seed at position 5 counts 10, showing the undercount is confined to position 0
• window_read_coverage_union_of_overlapping_seeds: overlapping seeds count their union, not their sum
• window_read_coverage_disjoint_seeds

In src/align/score.rs:

• the_stitch_mismatch_cap_does_not_apply_to_long_reads
• the_scorer_takes_long_read_from_the_parameters

Gate: 568 lib + integration tests, cargo clippy --all-targets -- -D warnings, cargo fmt --check, all green. Rebased onto current main, which also restores DIVERGENCE.md: the branch predated that file, so the diff previously read as deleting it.