[BenjaminDEMAILLE]

Part of #354 (findings 10, 14 and 15). Independent of the other open PRs apart from CHANGELOG.md and DIVERGENCE.md.

Three defects in how the index is built, updated and loaded. One commit, because all three share tests/index_load.rs; I can split it if that is easier to review.

1. A junction near a genome edge stopped the run

Problem

build_gsj (src/junction/sjdb_insert.rs) copies --sjdbOverhang bases on each side of a junction. It returned an error when the donor window started before position 0 or the acceptor window ran past the genome's end. On any later contig, the same flank just reads the spacer bytes between chromosomes, and the rest of the pipeline already handles those.

So the failure depended on which contig the junction was on:

• A GTF with an exon at 1-50 on the first contig and overhang 100 failed at genomeGenerate.
• --twopassMode Basic --sjdbOverhang 100 with a novel junction at 81-900 on the first contig finished pass 1, then exited 1 with donor underflows. That loses all of pass 1.

What STAR does

Observed with Homebrew STAR 2.7.11b: both cases succeed. A donor flank before the start reads the spacer padding before G; for the end, see §3.5 below.

Change

Each flank copies only the bases that exist; the missing part keeps the spacer byte (GSJ_SPACING) the buffer starts with. An intron that itself runs past the genome is still an error. As before, overhangs above 200 are not handled specially.

Measured

Two contigs of 4094 bases, --genomeChrBinNbits 12 --sjdbOverhang 100, each tool with its own index:

Case	main	This branch	Against STAR
Exon 1-50 and 500-700 on c1	exit 1	exit 0, 200/200 reads map	Genome, SA, SAindex, sjdbList.out.tab, sjdbInfo.txt byte-identical
Exon 4000-4094 on the last contig	exit 1	exit 0, 200/200 reads map	sjdbList.out.tab, sjdbInfo.txt identical; Genome differs by one byte, see below
Two-pass, novel junction at c1 81-900	exit 1 after pass 1	exit 0	Pass 1 SJ.out.tab lists c1 81 900 in both

Deliberate divergence: DIVERGENCE.md §3.5

The first commit wrote spacers for an acceptor flank past the genome's end, where STAR's Genome held 5 5 0. Reading STAR's source explained the 0: at genomeGenerate STAR keeps the reverse strand right after G (Genome_genomeGenerate.cpp:176-181) and copies flanks straight from G (sjdbPrepare.cpp:210-211), so the flank reads comp(G[2n-1-p]). At mapping, that strand is not in memory and STAR reads spacers.

c18748f does the same (build_gsj, mirror_past_end). Last contig of 4094 bases, exons 3000-3100 and 4000-4094, --sjdbOverhang 100 --genomeChrBinNbits 12:

Last base	main	first commit	c18748f
T	exit 1	Genome, SA, SAindex differ	Genome, SA, SAindex, sjdbList.out.tab, sjdbInfo.txt identical
G	exit 1	Genome, SA, SAindex differ	the same identical, except 2 bytes of SAindex

I did not find what the 2 SAindex bytes come from; §3.5 now describes only them, and still needs maintainer sign-off.

2. On-the-fly insertion left a wrong SAindex when an old junction was dropped

Problem

When a junction list is reloaded, two entries can collapse into one. For example, c1 5001 6000 . and c1 5000 5999 + on a repeat are distinct at genomeGenerate, and the same junction once the first gets its strand from its motif. insert_junctions (src/index/sjdb_otf.rs) then removes that block's suffixes from the SA, but update_sa_index, a literal port of STAR's in-place update, only ever moves entries up. Every SAindex value after a removed block was left too high.

The result was seeds pointing at the wrong SA range. At the end of the table, it was a panic at index < self.length, reported as alignment task panicked. Without a new junction, nothing was inserted and mapping was fine.

Change

When dropped_old is true, the SAindex is rebuilt with SaIndex::build from the new SA, keeping sparse_d. The common path, no dropped junction, is unchanged and still updates in place. I did not compare this with STAR, because I did not establish what STAR does when it drops a block.

Measured

Reproduced on a 40 kb genome with 4000 reads: main panics. On this branch, SA and SAindex under <prefix>_STARgenome/ are byte-identical to a fresh genomeGenerate of the surviving junctions. The new test fails without the rebuild and passes with it; I checked that by disabling the branch.

A rebuild costs one pass over the SA, only on this rare path.

3. A damaged index loaded without complaint

Problem

load_suffix_array derived the SA's length from its file size and ignored the size recorded in genomeFileSizes. load_sa_index trusted its header, and PackedArray::read returns 0 past the end of the data. Chromosome tables were parsed with unwrap.

Damage	main
SA 40 bytes short	exit 0, entries silently missing
SA cut in half	panic mid-alignment
SAindex cut to 70 %	exit 0
Garbage in chrLength.txt	unwrap panic
chrStart.txt too short	index panic

The Genome file was already checked. An interrupted copy to a cluster's scratch disk is the realistic way to get here.

Change

All in src/index/io.rs:

• read_genome_file_size becomes read_genome_file_sizes, returning both values. A recorded SA size of 0 means unknown, since genome/mod.rs writes that placeholder before the SA exists.
• SA must match its recorded size exactly. If none is recorded, nothing is checked, as before.
• SAindex data must be at least PackedArray::data_byte_len_for bytes.
• parse_u64_lines replaces the unwraps, and the chromosome tables must agree in length.

Errors name the file and both sizes, and say to copy again or regenerate.

Checked for false refusals

These still load and map: a STAR-built index (genomeFileSizes as STAR writes it), and a <prefix>_STARgenome/ written with --sjdbInsertSave All, whose recorded size equals its SA.

Tests

tests/index_load.rs:

• short_sa_is_refused, short_saindex_is_refused, garbage_chr_length_is_refused_without_a_panic, short_chr_start_is_refused_without_a_panic
• junction_at_the_genome_start_builds_and_maps
• insertion_that_drops_an_old_junction_keeps_a_valid_saindex

src/junction/sjdb_insert.rs: build_gsj_pads_a_donor_flank_before_the_genome_start, build_gsj_pads_an_acceptor_flank_past_the_genome_end, build_gsj_errors_when_the_intron_overruns_the_genome. These replace build_gsj_errors_when_flank_underflows and build_gsj_errors_when_flank_overruns_genome, which locked in the old error.

cargo test --release (820 passed), cargo clippy --all-targets and cargo fmt --check are clean on macOS. I did not run test/parity_suite.sh. Output is unaffected unless a run hits one of these cases.

🤖 Generated with Claude Code