Skip to content

Refuse a junction or annotation table cut short - #374

Open
BenjaminDEMAILLE wants to merge 1 commit into
fix/junction-index-safetyfrom
fix/annotation-table-safety
Open

BenjaminDEMAILLE wants to merge 1 commit into
fix/junction-index-safetyfrom
fix/annotation-table-safety

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

Part of #354 (truncated junction and annotation tables in the index).

Stacked on #357. The base is fix/junction-index-safety; only the last commit is this PR.

Problem

#357 checks the sizes of Genome, SA, SAindex and the chromosome tables. The text tables were not covered:

  • A last line cut mid-value in sjdbInfo.txt, exonInfo.tab or geneInfo.tab loaded, with a wrong last value (for example an exLenCum of 20 instead of 200). Exit 0.
  • sjdbList.out.tab and exonGeTrInfo.tab were never read, so a short one went unnoticed. Exit 0.
  • The junction count in sjdbInfo.txt was not checked against the Genome.

Whole lines missing from sjdbInfo.txt, transcriptInfo.tab, exonInfo.tab and geneInfo.tab were already refused by the row-count checks.

What STAR does (2.7.11b source)

  • sjdbInfo.txt: Genome::loadSJDB reads the count, the overhang and the rows with >> and checks nothing (Genome_genomeLoad.cpp:487-519). Missing rows read as zeros.
  • sjdbList.out.tab: read only when junctions are inserted at mapping (--sjdbGTFfile, --sjdbFileChrStartEnd, --twopassMode Basic) to reload the index's junctions (sjdbInsertJunctions.cpp:15-18). sjdbLoadFromStream.cpp reads lines until EOF, with no count.
  • transcriptInfo.tab: refused when a read fails, exit 105 (Transcriptome.cpp:51-56). A last line without its newline also fails: the last >> hits EOF.
  • geneInfo.tab, exonInfo.tab, exonGeTrInfo.tab: read with no check (Transcriptome.cpp:18-29, 62-75, 78-97).
  • STAR reads the transcript tables only with --quantMode or STARsolo (Transcriptome.cpp:9).

Change

  • read_index_table (src/index/io.rs): a table that does not end with a newline is refused, "is empty or its last line is cut: the index is damaged or was not fully copied". Used for sjdbInfo.txt, sjdbList.out.tab, transcriptInfo.tab, exonInfo.tab, geneInfo.tab and exonGeTrInfo.tab. Every one is written with a final newline by STAR and by rustar-aligner.
  • read_sjdb_info_tab: the junctions, 2*sjdbOverhang+1 bytes each, have to fill the Genome after the chromosomes.
  • check_sjdb_list: when junctions are inserted, sjdbList.out.tab has to hold as many junctions as sjdbInfo.txt. Without insertion it is not read, as in STAR.
  • check_exon_ge_tr_info: exonGeTrInfo.tab, when present, has to hold the rows its header says. rustar-aligner rebuilds this table from the others, but STAR counts genes from it, so a short one gives other numbers in STAR on the same index.
  • sjdb_insert(params) factors out the sjdbInsert.yes test that check_genome_version already had.

Where STAR exits 0 this is a divergence: DIVERGENCE.md §1.9, needs maintainer sign-off.

Measured against STAR 2.7.11b

Synthetic genome of 60 kb, GTF of 20 genes, 40 transcripts, 100 exons, 60 junctions, --sjdbOverhang 99; 400 spliced 100 bp reads. Each tool built its own index; all six tables are byte-identical between the two. STAR is the Linux 2.7.11b build for --quantMode. Mapping was plain for sjdbInfo.txt, with --sjdbFileChrStartEnd (one new junction) for sjdbList.out.tab, and with --quantMode TranscriptomeSAM GeneCounts for the transcript tables. "line" drops the last line, "mid" the last 3 bytes, "10 rows" the last 10 lines.

Intact index, both tools: rc 0, 400 unique, 40 junctions in SJ.out.tab, same SAM records, 400 reads in Aligned.toTranscriptome.out.bam, same ReadsPerGene.out.tab.

Table, cut STAR rustar before rustar after
sjdbInfo.txt line rc 0 rc 1, header 60 found 59 rc 1, same
sjdbInfo.txt mid rc 0 rc 1, bad strand rc 1, last line is cut
sjdbInfo.txt 10 rows rc 139 (segfault) rc 1 rc 1
sjdbInfo.txt overhang 99 → 98 in header rc 0, 210 of 400 unique, 73 junctions not measured rc 1, do not fill
sjdbList.out.tab line, insertion 100% CPU, killed after 10 min rc 0 rc 1, has 59 junctions but sjdbInfo.txt has 60
sjdbList.out.tab mid, insertion rc 0 rc 0 rc 1, last line is cut
sjdbList.out.tab line, no insertion not read rc 0 rc 0
transcriptInfo.tab line / mid rc 105 rc 1 / rc 1 rc 1 / rc 1
exonInfo.tab line rc 0 rc 1 rc 1
exonInfo.tab mid rc 0 rc 0 rc 1, last line is cut
exonInfo.tab 10 rows rc 0, 360 reads in toTranscriptome rc 1 rc 1
geneInfo.tab line rc 0, other ReadsPerGene.out.tab rc 1 rc 1
geneInfo.tab mid rc 0 rc 0 rc 1, last line is cut
exonGeTrInfo.tab line / mid rc 0 / rc 0 rc 0 / rc 0 rc 1 / rc 1
exonGeTrInfo.tab 10 rows rc 0, 40 reads N_noFeature instead of 0 rc 0 rc 1

A cut of 2 bytes (the newline and one digit) gives the same "after" row as "mid" for all six tables.

Tests

tests/index_load.rs, on an index built from a one-transcript GTF:

  • table_with_a_cut_last_line_is_refused: each of the six tables loses 2 bytes, mapping with junction insertion fails naming the file; the restored index maps.
  • short_sjdb_list_is_refused_when_junctions_are_inserted: one junction of two left; plain mapping succeeds, --twopassMode Basic fails.
  • sjdb_info_that_disagrees_with_the_genome_is_refused: overhang 50 → 49 in the header.

All three fail without the src/ change and pass with it. read_sjdb_info_tab_roundtrip now gives its genome the junction bytes. cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

  • A cut in the middle of a table (a missing row that is not the last) whose row count is still right. Neither tool can see it.
  • A table edited by hand that lost its final newline is now refused.
  • rustar-aligner loads the transcript tables, and so refuses a damaged one, even without --quantMode or STARsolo, where STAR does not read them. That was already so for missing lines before this PR.
  • STARsolo runs and Velocyto with a damaged table were not run.
  • Why STAR hangs on a short sjdbList.out.tab was not traced.

🤖 Generated with Claude Code

sjdbInfo.txt, sjdbList.out.tab, transcriptInfo.tab, exonInfo.tab,
geneInfo.tab and exonGeTrInfo.tab must end with a newline and hold the
rows their header (or sjdbInfo.txt) says. sjdbInfo.txt's junctions must
fill the Genome after the chromosomes. sjdbList.out.tab is checked only
when junctions are inserted, the one time STAR reads it. STAR reads most
of these silently and maps or counts wrongly, or hangs: DIVERGENCE.md 1.9.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant