Tracking: fixes from the whole-repo audit (main at 32fbbc1) #354

Tracking the fixes from a whole-repo audit of main at 32fbbc1. The audit report is below, word for word. Numbers 1-20 are its findings; the 10 smaller items it left out are tracked as 21-30. A new row is a difference found while fixing these, not in the report.

Status

Open PRs touch the same code: alignment fixes (stitch.rs, read_align.rs) stack on #346, and FASTQ input fixes stack on #353. Each PR states its base.

Audit report

What this repo does: rustar-aligner is a Rust copy of STAR, a tool that finds where RNA sequencing reads come from in a genome. It builds a genome index, aligns reads, and writes alignments, gene counts and single-cell (STARsolo) matrices, aiming to match STAR 2.7.11b exactly.

Load I assumed: one person or pipeline running one process with many threads, on hundreds of millions of reads and genomes of about 3 billion bases.

All 813 tests pass and clippy reports no warnings, but these bugs get past them. "STAR" below means real runs of the Homebrew 2.7.11b binary on small made-up genomes.

Must fix

1. Reads are lost on genomes with many small contigs (src/params/mod.rs:1732-1757)

   • What this is: This code sets the size of the genome "bins" used to group a read's seeds (short exact matches).
   • Problem: It returns before applying the size cap, and the cap reads the command-line value, not the index's. Contigs sharing a bin overwrite each other. With 400 contigs and --genomeChrBinNbits 11, STAR maps 3998 unique and 7 multi; this tool maps 3945 and 1, with 59 "too short".
   • Fix: Always apply the cap, using the value stored in the index.
   • If we skip it: Draft genomes silently lose reads, and repeat reads get reported as unique.

2. Half-mapped pairs report one location as unique (src/align/read_align.rs:1316-1376)

   • What this is: This is the fallback for a pair where only one mate aligns.
   • Problem: It keeps one location and skips the multimapper filters. A mate in a 3-copy repeat gets NH:3 MAPQ:1 on 3 loci from STAR, but one record with NH:1 MAPQ:255 here, on a different copy. It also adds a record for the unmapped mate, which STAR left out on the same options.
   • Fix: Send single-mate results through the same selection and filter steps as pairs.
   • If we skip it: Downstream tools treat repeat reads as high-confidence unique hits.

3. --alignEndsType EndToEnd is wrong for pairs (src/align/stitch.rs:1231-1304, src/align/read_align.rs:884)

   • What this is: This mode must align every read to its very end, with no clipping.
   • Problem: Inner mate ends are never forced, and outer ends don't charge their penalty to the pair score. With 4 mismatches at each inner end, STAR gives 100M/100M AS:182; this tool gives 96M4S/4S96M AS:190.
   • Fix: Pass the end-to-end flags at the mate-boundary and outer extensions, so the pair score carries the penalty.
   • If we skip it: Paired runs in this mode get clipped reads and inflated scores, even though the user asked for no clipping.

4. Solo runs silently write empty or all-zero side outputs (src/lib.rs:906-913, src/lib.rs:1021-1029, src/lib.rs:2461)

   • What this is: Solo runs use their own read loops, which don't receive the quant or unmapped writers.
   • Problem: With --quantMode GeneCounts TranscriptomeSAM --outReadsUnmapped Fastx, I got exit 0, all zeros in ReadsPerGene.out.tab, an empty transcriptome BAM and an empty unmapped file. BySJout, chimeric and WASP options have no code in those loops either. Only a code comment mentions part of this.
   • Fix: Reject these combinations at startup now, then add support where needed.
   • If we skip it: A user can analyse zeros without knowing it.

5. A failed decompress command still exits 0 (src/io/fastq.rs:167-181)

   • What this is: --readFilesCommand runs a program that decompresses the reads.
   • Problem: Its exit status is never checked. --readFilesCommand gzcat --readFilesIn missing.gz gave exit 0, 0 reads and a clean log. With r1.gz,missing.gz the second file was dropped, while STAR exits 102.
   • Fix: Wrap the child's output so that end of file calls wait() and returns an error on failure.
   • If we skip it: A pipeline can finish "successfully" with reads missing.

6. The GTF parser fails the wrong way in both directions (src/junction/gtf.rs:155-159, src/junction/gtf.rs:190-213, src/index/mod.rs:311-319)

   • What this is: This code reads the gene annotation file.
   • Problem: A GTF whose chromosome names don't match the genome (X vs c1) exits 0 with no junctions, while STAR exits 104. GFF3 Parent= attributes and a single exon with no transcript ID stop the build here, while STAR accepts both.
   • Fix: Fail if no exon lands on a known chromosome, split attributes on a space or =, and warn and skip on a missing ID.
   • If we skip it: The common chr naming mistake gives an unannotated index and all reads as N_noFeature.

7. A transcript ID on two chromosomes invents junctions (src/junction/gtf.rs:187-225)

   • What this is: Exons are grouped into transcripts to derive junctions.
   • Problem: Grouping uses the ID only, then places everything on the first exon's chromosome. My test gave three false c1 junctions and lost the real c1 1201-1999 and c2 1701-2499. STAR keeps both real ones, though it also writes one odd junction of its own.
   • Fix: Group by chromosome and ID. This won't copy STAR's odd junction, so add a DIVERGENCE.md entry.
   • If we skip it: False junctions get annotated scoring and affect alignment.

8. A whitelist of the wrong length gives an empty matrix and exit 0 (src/solo/whitelist.rs:498-559)

   • What this is: This code matches cell barcodes against the allowed list.
   • Problem: Entry length is never compared with barcode length. A 14-base whitelist with --soloCBlen 16 gave 0 valid barcodes and exit 0, while STAR exits 104. Leading As pack to zero, so some reads can falsely match.
   • Fix: Check the length at startup and in match_cb.
   • If we skip it: Empty results look like a failed experiment, and false matches give wrong cells.

9. CB_UMI_Complex handles UMIs wrongly (src/solo/mod.rs:211-216, src/solo/count.rs:1537, src/solo/count.rs:1878)

   • What this is: A UMI is a short tag that identifies one original molecule.
   • Problem: A missing --soloUMIposition becomes an empty UMI, so 200 reads counted as 1 in my run, with exit 0. A bad CB position is dropped the same way. UMI comparison also uses --soloUMIlen (default 10), not the position's actual length, so 12-base UMIs merge wrongly.
   • Fix: Make bad positions a fatal error, and take the UMI length from the position.
   • If we skip it: Counts are silently far too low.

10. Junctions near a genome edge stop the run (src/junction/sjdb_insert.rs:530-540)

    • What this is: Inserting a junction copies the sequence on each side of it.
    • Problem: An exon at 1-50 on the first contig with overhang 100 exits 1 with "donor underflows", while STAR succeeds. A reviewer also reproduced this after pass 1 of a two-pass run, and with a junction near the genome's end.
    • Fix: Leave the out-of-range part as spacer bytes; the buffer already starts that way.
    • If we skip it: A long two-pass run can die after its first pass.

11. SAM output can lose its last bytes with exit 0 (src/lib.rs:294-301, src/lib.rs:336-343)

    • What this is: BAM writers have a finish step; SAM writers use an empty default.
    • Problem: The last buffer is written when the writer is dropped, and that drop ignores errors. With a file-size limit, 20 reads gave a file with 14 records, exit 0, and a log saying 20 mapped.
    • Fix: Add a finish that flushes and returns the error, to both SAM writers.
    • If we skip it: A full disk at the end of a run gives a cut file that looks complete.

12. SJ.out.tab is missing when no junctions are found (src/lib.rs:1035, src/lib.rs:1103)

    • What this is: This file lists splice junctions.
    • Problem: It is skipped when empty, while STAR writes an empty file. I confirmed this by running both.
    • Fix: Delete && !sj_stats.is_empty() in both places.
    • If we skip it: Pipelines that expect the file fail on small or unspliced samples.

Should fix

13. Chimeric detection can hang on one read (src/chimeric/detect.rs:665-668, src/chimeric/detect.rs:734-743)

    • What this is: place() trims a segment at a chimeric (fusion) junction.
    • Problem: A negative length is cast to usize, giving a huge number. With --chimJunctionOverhangMin 1 and permissive options, one read hung past a 20 s limit. The default overhang of 20 looks safe by reading the check.
    • Fix: Return None when jr_max is 0 or the new length would not be positive.
    • If we skip it: The job stalls with no error until the scheduler kills it.

14. The seed lookup table goes wrong when an old junction is dropped (src/index/sjdb_otf.rs:354-386)

    • What this is: The SAindex is a lookup table into the suffix array (the sorted list of genome positions), updated when junctions are added on the fly.
    • Problem: Removed entries aren't subtracted. A reloaded index holding the same repeat junction on opposite strands, plus one new junction, panicked. The trigger is rare.
    • Fix: Rebuild the SAindex when dropped_old is true.
    • If we skip it: A rare crash, and seeds can be misdirected before it happens.

15. A damaged index loads without complaint (src/index/io.rs:263-284, src/index/io.rs:364-422, src/index/packed_array.rs:152-163)

    • What this is: This code loads the index files from disk.
    • Problem: The sizes recorded in genomeFileSizes are ignored. An SA file 40 bytes short gave exit 0 with entries missing, and one cut in half panicked mid-alignment. Garbage in chrLength.txt also panics.
    • Fix: Compare the file sizes on load, and return errors instead of unwrapping.
    • If we skip it: An interrupted copy gives silently incomplete results.

16. --readFilesCommand takes only one word (src/params/mod.rs:537-538, src/io/fastq.rs:168)

    • What this is: STAR accepts a command with arguments here.
    • Problem: gunzip -c gives a usage error, and - tries to run a program called -. STAR accepts gzcat -c and -.
    • Fix: Use Vec<String> with num_args = 1.., and treat - as no command.
    • If we skip it: Commands copied from existing STAR pipelines fail, though the error is at least loud.

17. FASTQ edge cases fail badly (src/io/fastq.rs:199-243, src/io/fastq.rs:472)

    • What this is: This is the read parser.
    • Problem: A blank line after the last record kills the run at the very end, while STAR finishes. A quality line shorter than its sequence panics with no read name. FASTA reads are also rejected, while STAR accepts them.
    • Fix: Treat a leading blank line as end of file, check the lengths, and name the read in the error. Then support FASTA or document that it is rejected.
    • If we skip it: Hours of work are lost to one trailing newline.

18. Sort temp files land in the wrong folder (src/io/bam.rs:164-173)

    • What this is: Coordinate sorting spills to temp files next to the output.
    • Problem: Path::new("out/sub/").parent() is out, so --outFileNamePrefix /data/run1/ spills into /data.
    • Fix: Return the prefix itself when it ends in /.
    • If we skip it: Spills fill the wrong disk, and a read-only parent should fail the run at the first spill (reasoned, not run).

19. Four read loops are near-copies (src/lib.rs:1828-4216)

    • What this is: Single-end, paired-end, solo and solo paired each have their own loop of about 450 to 800 lines.
    • Problem: A rule added to one doesn't reach the others, which is how finding 4 happened. Solo and solo paired share 218 distinct lines; single-end and paired-end share 310.
    • Fix: Pull the shared per-read output and writer stage into one function, by job, without adding a framework.
    • If we skip it: Every new feature must be added up to four times, and the copies keep drifting apart.

20. Dead code on the hot path and elsewhere (src/align/stitch.rs:3317-3452, src/io/bam.rs:17-35, src/io/sam.rs:172-231, src/io/bam.rs:646-708)

    • What this is: This code doesn't affect output, or runs only in tests.
    • Problem:
      • A 135-line per-seed scoring block ends in let _ = best_pre_score. Deleting it in a scratch copy kept SAM md5 identical on 800k and 300k reads.
      • The reverse-strand branches are unreachable from the production path, because stitch.rs:3266-3273 forces the flag to false.
      • BufferedBamRecords has no users, and SamWriter::write_alignment is only called by its own test.
      • find_paired_seeds and get_genome_positions (src/align/seed.rs) are test-only.
      • The stdout writers repeat the file writers.
    • Fix: Delete the dead code, and make the writers generic over their output.
    • If we skip it: People keep reading and maintaining code that does nothing, and one block runs for every read.

Left out (10 smaller items):

• --soloUMIlen above 32 corrupts UMIs.
• Extra --soloUMIdedup values are dropped.
• soloCellFiltering writes gzip bytes into features.tsv.
• UMI deduplication is quadratic (200k UMIs in one gene took 104 s).
• bamRemoveDuplicates holds the whole BAM in memory.
• The gene name tag is ignored when the GTF is given at mapping.
• CellReads.stats mixes features.
• .1/.2 mate names are rejected.
• --outBAMcompression -1 writes uncompressed BAM.
• The docs are stale: docs-old/, and the README still says v0.1.

Verdict: Not healthy yet for general use: the main flows match STAR, but less-tested options give wrong results with exit 0. Fix 1-3 first, then 4, 5 and 11.

Lean: about -400 lines (my estimate, more with the reverse branches and finding 19), -0 dependencies possible.

Not checked:

• Fixes: none has been run, except deleting the dead block in finding 20.
• Real data: nothing was compared on a real genome. Everything used small made-up genomes, and I did not run test/parity_suite.sh.
• Scale: genomes over 4 billion bases were not run. An open question remains on GstrandBit near 2^32 (src/genome/mod.rs:568).
• Reproduced by me: findings 1-9 and 12, except 8's false matches and 9's UMI-width bug, plus the GTF part of 10.
• Reproduced only by the reviewers I sent out: 11, 13-18, the two-pass and genome-end parts of 10, and the dead-block test in 20. For 13, the hang was run, but its cause comes from reading the code.
• Not read:
  • the src/io/sam.rs record builders
  • src/wasp/mod.rs from line 480
  • src/junction/sj_output.rs
  • src/align/transcript.rs and src/chimeric/segment.rs
  • src/index/sa_build.rs after line 900
  • src/solo/sgt.rs after line 40, and src/solo/libcxx_rng.rs
  • src/degnorm/mod.rs
  • most test modules, and the bz2, zstd and xz input paths at runtime
• Solo counts: a reviewer could not get STAR's solo runs going, so solo counts were never compared with STAR.

🤖 Generated with Claude Code