[BenjaminDEMAILLE]

Part of #354 (findings 1 and 13).

Claims about STAR below are observed output of Homebrew STAR 2.7.11b; source not read. Each tool used its own index.

1. Window bins on genomes with many small contigs (finding 1)

Problem

redefine_window_params had two faults:

• It returned before the genomeChrBinNbits cap when --alignIntronMax and --alignMatesGapMax were both 0, which is the default.
• The cap read the command-line value, not the index's. --genomeChrBinNbits is normally given only to genomeGenerate, so at mapping time the cap saw the default 18. loaded_genome_params in src/index/mod.rs reads the index's value, but only into a local copy.

So an index built with small chromosome bins kept 2^16 window bins. Contigs then shared a bin: reads were lost as "too short", and repeats were reported as unique.

What STAR does

The comment "STAR: no redefinition when both are 0" is half true, and I removed it. Log.out on the 400-contig index:

Options	STAR Log.out	winBinNbits
none	winBinNbits=16 > pGe.gChrBinNbits=11 redefining:	11
--genomeChrBinNbits 18 at mapping (index has 11)	same lines	11
--winBinNbits 5	no redefinition	5
--alignIntronMax 100000	redefined winBinNbits=15, then capped; winFlankNbins=49 and winAnchorDistNbins=98	11
--alignIntronMax 1450	redefined winBinNbits=8; winFlankNbins=6 and winAnchorDistNbins=12	8
--alignMatesGapMax 1450	winBinNbits=8; winFlankNbins=6 and winAnchorDistNbins=12	8
--alignIntronMax 1000, nGenome 14,450,688, index bins 12	winBinNbits=8; winFlankNbins=4 and winAnchorDistNbins=8	8

So with both limits 0, only the cap applies; winFlankNbins and winAnchorDistNbins stay as given. The last two rows also show integer divisions: 1450/4 and 14450688/40000 computed as floats would round to 9. main used floats, so I changed that too. It is in the same formula and is one more reason the bins differed from STAR.

STAR's own line alignIntronMax=alignMatesGapMax=0, the max intron size will be approximately determined by (2^winBinNbits)*winAnchorDistNbins=589824 is printed before the cap, so it can state 589824 when the real figure is 18432.

Change

• redefine_window_params always applies the cap, and flanks are derived after it, only when a limit is set.
• run_align sets genome_chr_bin_nbits from the index before calling it, using loaded_genome_params, now pub(crate). This is the only caller.
• win_bin_window_dist had no caller; the log line now uses it. I removed its comment, which described uses that do not exist.
• Comments in src/align/read_align.rs and tests/alignment_features.rs said winBinNbits*winAnchorDistNbins (~589 kb). It is 2^winBinNbits, and the cap can lower it.

Measured against STAR

400-contig fixture, --genomeChrBinNbits 11, 4005 single-end reads. "Differing" is the count of lines in comm -3 of sorted SAM fields 1-9 against STAR.

Mapping options	Tool	Unique	Multi	Too short	Records	Differing
none	STAR	3998	7	0	4024
	main	3945	1	59	3947	89
	this branch	3998	7	0	4024	0
--alignIntronMax 100000	STAR	3998	7	0	4024
	main	3977	1	27	3979	57
	this branch	3998	7	0	4024	0
--alignIntronMax 1450	STAR	3998	7	0	4024
	main	3998	7	0	4024	0
	this branch	3998	7	0	4024	0
--alignMatesGapMax 1450	STAR	3998	7	0	4024
	main	3998	7	0	4024	0
	this branch	3998	7	0	4024	0
--winBinNbits 5	STAR	3999	6	0	4023
	main	3999	6	0	4023	0
	this branch	3999	6	0	4023	0

The target numbers are reached. The rounding change does not move any count on this fixture; it is only checked through winBinNbits in the table above, which is what window_params_follow_star_log asserts.

2. Chimeric detection hang (finding 13)

Problem

In place(), new_l was computed as i64 and cast to usize. When the trailing segment's first block ends before the leading segment's last block, the value is negative and becomes a huge length. Loops over the block length then do not finish in any useful time.

Reproduced on main with the audit's read pair (a 153-base mate 1 and a 12-base mate 2) on its 29 kb genome, --chimSegmentMin 10 --chimMultimapNmax 20 --chimScoreDropMax 200 --chimNonchimScoreDropMin 0 --chimScoreSeparation 0, under timeout:

--chimJunctionOverhangMin	STAR	main	This branch
1	bus error, exit 138	killed at the 20 s limit, exit 124	exit 0, no chimera
2	exit 0, no chimera	exit 0, no chimera	exit 0, no chimera

Chimeric.out.junction counts exclude its header and comment lines.

Change

place() returns None only when the trailing block's first block ends before the leading segment's last block starts (roStart1 + L <= roStart0), the case whose new length wraps in STAR's uint (ReadAlign_chimericDetectionOld.cpp:154-155, :253-256). A block ending exactly there, or jr_best == jr_max because the scan skips offset readLength[0], gives length 0, which STAR keeps; so does place() since 74d8a02 (the first version rejected it too, which review against STAR's source caught). The overhang filters decide it later, so this only matters with --chimJunctionOverhangMin 0.

Measured against STAR

Random multi-segment reads on the same genome, with repeats and canonical introns: 4000 SE, 4000 PE of which about 30% have a 12-base mate 2. Chimeric.out.junction lines, sorted:

Run	STAR	main	This branch	Differing from STAR
SE, old detector, overhang 1	1766	1766	1766	0
SE, --chimMultimapNmax 20, overhang 1	11234	11234	11234	0
PE, old detector, overhang 1	1376	1376	1376	0
PE, --chimMultimapNmax 20, overhang 1	bus error, exit 138	killed at the 60 s limit, exit 124	9576, exit 0	not comparable
PE, --chimMultimapNmax 20, overhang 2	9556	9556	9556	0
PE, --chimMultimapNmax 20, overhang 3	9523	9523	9523	0
PE, --chimMultimapNmax 20, default overhang	9133	9133	9133	0

Deliberate divergence, needs maintainer sign-off

DIVERGENCE.md §1.9: where STAR dies with a bus error, this skips the candidate and exits 0. There is no STAR output to match in that case, so the 9576 lines above cannot be checked against STAR.

Tests

Each failed when I took out the corresponding fix, and passes with it:

• src/params/mod.rs: window_params_follow_star_log, with the expected values from the STAR Log.out table above.
• tests/alignment_features.rs: test_window_bins_capped_by_the_index_chr_bins, 200 contigs, index built with --genomeChrBinNbits 11, not repeated at mapping. Expected: 800 unique, 1 multi with 4 records. STAR gives those numbers on the same data, with 0 differing lines in SAM fields 1-9; main gives 801 and 0, with 1 record.
• src/chimeric/detect.rs: place_rejects_a_trailing_block_ending_before_the_last_leading_block, place_keeps_a_junction_leaving_no_trailing_bases.

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

Not checked

• test/parity_suite.sh was not run.
• The 10k yeast SE and PE harnesses in CONTRIBUTING.md were not run. Their index uses the default 18-bit bins and their runs set no intron or gap limit, so by reading the code I expect no change, but that is not measured.
• --chimOutType WithinBAM was not compared on the reads that hang.
• The end-to-end hang is not a test in this PR. The regression tests call place() directly and assert None.
• Only macOS was run.

🤖 Generated with Claude Code