fix(quant): Aligned.toTranscriptome.out.bam record for record as STAR (primary draw, pair-level soft-clip budget, MAPQ, mate order, merged CIGAR, deletions, alignment order)

[BenjaminDEMAILLE]

Stacked on #311 (genomic alignments matching STAR), with #306 (transcriptome BAM carries only NH HI [RG]) and #313 (transcript tables in STAR's order) merged into the branch; their changes show in this diff until they land. The transcriptome BAM header (no @HD/@pg) is #297's.

Finding

On the yeast 200k pairs (annotated index, --quantMode TranscriptomeSAM GeneCounts) the transcriptome BAM differed from STAR 2.7.11b in primary flags, dropped/extra alignments, MAPQ, mate order, CIGARs (25M125M vs 150M), HI order and (with indels allowed) missing deletion alignments.

Causes and fixes

• Primary flag: STAR draws int(rngUniformReal0to1(rngMultOrder)*nAlignT) once per mapped read (also when nAlignT is 0), from a per-thread mt19937 seeded runRNGseed*(iChunk+1) (ReadAlign.cpp, ReadAlign_quantTranscriptome.cpp). Workers now emit records flagged secondary and the ordered writer stage draws from the libc++ Mt19937 (src/solo/libcxx_rng.rs) seeded with runRNGseed, so the result equals STAR --runThreadN 1 for any rustar thread count. STAR with several threads is not reproducible (chunk to thread assignment varies); documented, not reproducible by design.
• Soft-clip extension budget: STAR sums the extension mismatches of both mates and the pair nMM and compares with min(outFilterMismatchNmaxTotal, outFilterMismatchNoverLmax*(Lread-1)), Lread = len1+len2+1. rustar tested each mate alone. New filter_and_project_pair.
• Indel ban uses the CIGAR (n_gap misses insertions of stitched pairs).
• MAPQ is from nAlignT (alignBAM nTrOut), not the genomic multimapper count.
• Mate order: STAR writes the left segment of the projected alignment first (mate2 first when the projected mate1 is reverse).
• CIGAR: blocks either side of a junction are one block in transcript space (alignToTranscript, canonSJ>=0 merges), so adjacent M merge. Deletions and insertions (canonSJ -1/-2) are not junctions: no junction match required, block stays in the exon. Gap kinds are read off the CIGAR (block_gaps).
• Alignment order: STAR walks trMult (window order). The genomic output sorts by score and position (unchanged); pairs now carry star_order (rank before that sort) and the transcriptome follows it.

Results (rustar vs STAR, samtools view, identical/total records)

• default: 20k 27468/27468, 200k 274542/274542 (1 and 8 threads, in order)
• BanSingleEnd: 20k 31806/31806, 200k 317598/317598
• BanSingleEnd_ExtendSoftclip: 20k 28198/28198, 200k 282338/282338
  STAR 8 threads vs STAR 1 thread: 244946/274542 (primary flags only), confirming the nondeterminism.
  Genomic Aligned.out.sam on the 200k pairs: records byte-identical to the starting build.
  --quantTranscriptomeBan does not exist in 2.7.11b; its Singleend is --quantTranscriptomeSAMoutput BanSingleEnd.

Single-end and BySJout

• Single-end: Transcript gets star_order (rank before the genomic score/position sort, which is unchanged); the SE transcriptome follows it. Single-end r1_200k vs STAR 1 thread: default 20k 16147/16147, 200k 161592/161592; BanSingleEnd_ExtendSoftclip 20k 16392/16392, 200k 164574/164574.
• --outFilterType BySJout: STAR holds a read with an unannotated junction (outFilterBySJout), writes every other read first, then maps the held reads again with only the filtered novel junctions (stitchWindowAligns.cpp, outFilterBySJoutStage==2) and draws for them then. The writer now writes non-held reads in read order, then the held reads, re-mapped with the second-stage junction filter (thread-local novel set checked in passes_finalization_filters; the genomic output is untouched). A junction counts as annotated if the stitcher flagged it or its CIGAR coordinates are in the sjdb. PE vs STAR 1 thread: 20k 27466/27470, 200k 274526/274606 (default mode unchanged at 274542/274542).
• Remaining BySJout gap: reads where rustar's outSJfilter* survivor set differs from STAR's (10 reads of 200k), and the primary of 20 reads that follow them; the genomic BySJout output has the same cause.

Tests

cargo test --release (609 unit tests, new: deletion kept in one exon, pair-level budget, merged CIGAR across a junction), clippy -D warnings, fmt.

Checked again by the main session on 20k yeast reads: paired-end 27,468 and single-end 16,147 transcriptome records identical to STAR --runThreadN 1, in order, rustar on 8 threads.

Part of #298.

🤖 Generated with Claude Code