[Psy-Fer]

STARsolo Velocyto and stats files as STARsolo

Velocyto's spliced / unspliced / ambiguous matrices were computed with a rule of our own: per read, a junction meant spliced, an intronic block unspliced, anything else ambiguous. On the 10x v2 mouse chr19 fixture that gave 52 / 100 / 362 UMIs against STARsolo's 369 / 57 / 32. This ports STAR's algorithm. Raw and filtered matrices are now byte-identical to STARsolo 2.7.11b.

What STAR does, and what is ported

• Per read (Transcriptome_classifyAlign.cpp:227-258):
  • It applies only to uniquely mapped reads (nAlignG == 1) whose barcode and UMI pass.
  • The read is compared with every transcript that contains it on the strand --soloStrand keeps, using alignToTranscriptMinOverlap with velocyto's 5-base flank. feat(quant): unspliced targets in TranscriptomeSAM, splicing-status tag and GeneSplicing counts for bulk total RNA-seq #283 already ported that for bulk use (quant/splice_status.rs), and it's reused unchanged. I re-checked it line by line, including the binarySearch1 exon lookup.
  • Each compatible transcript is kept with STAR's reAnn1 type bits. STAR drops a transcript only when a block of a spliced read runs into one of its introns; it does not compare the read's junctions with the transcript's, and neither do we.
  • The read's blocks are its exons merged across indels, mate by mate. A mate pair is one alignment, as in E2.
• Per molecule (SoloFeature_countVelocyto.cpp):
  • The cell and UMI come from the Gene feature's readInfo (:33-38); reads with either undefined are skipped.
  • A molecule keeps the transcripts all its reads share, their type bits OR-ed (:43-77). Once that set is empty, later reads are skipped.
  • A molecule whose transcripts span two genes is dropped. Otherwise it is classified by STAR's model flags (:95-136), using the collapse feat(quant): unspliced targets in TranscriptomeSAM, splicing-status tag and GeneSplicing counts for bulk total RNA-seq #283 ported.
• Output (SoloFeature_outputResults.cpp): every (cell, gene) row appears in all three matrices, zeros included.
• Velocyto/filtered/, which wasn't written before (SoloFeature_cellFiltering.cpp:16-30): the Gene feature's called cells that Velocyto detected. A cell is detected when it's one of the whitelist matches of a read with a Velocyto record; under 1MM_multi, every candidate counts. With no such cell, nothing is written (:106-108).
• Memory. STAR keeps the records in a temporary file. Here each read's record is its index, its transcript count and one packed u32 per transcript (tr << 4 | bits) in one flat store, and the detected cells are a bitset. That's 12 bytes for a one-transcript read, kept for the whole run, instead of a 56-byte struct plus two heap allocations.

readInfo per feature

STAR keeps a readInfo per feature: CB/UB read the first feature's (samAttrFeature), Velocyto the Gene feature's (ParametersSolo.cpp:420). E2 kept one, for the first feature; it's now a slot per feature.

• The no-feature records that fill it go to their own feature's recorder.
• Read indices are tracked whenever any feature keeps readInfo, so Velocyto works with --outSAMtype None.
• They are stamped onto BAM records only for the CB/UB pass.

Refusals

• --soloFeatures Velocyto without Gene: STAR segfaults, since it indexes the Gene feature's readInfo unconditionally; we refuse it. Before, we silently added a Gene feature and wrote a Gene/ directory nobody asked for.
• Velocyto with --soloType SmartSeq: refused with STAR's message (ParametersSolo.cpp:238-241).

Validation against STARsolo 2.7.11b (10x v2 mouse chr19 fixture)

Raw and filtered spliced/unspliced/ambiguous.mtx, features.tsv and barcodes.tsv are byte-identical in all 10 configurations:

• default (Gene, GeneFull, SJ, Velocyto);
• --soloStrand Unstranded, and Reverse;
• --soloCBmatchWLtype 1MM_multi;
• --soloFeatures GeneFull Gene Velocyto, so the two readInfos differ;
• CB UB GX GN in a sorted BAM; the BAM is identical too;
• --soloUMIdedup 1MM_Directional, and Exact;
• --soloCellFilter EmptyDrops_CR;
• paired-end cDNA with a separate barcode read (synthetic 45+45 pairs).

In every run, Gene, GeneFull and SJ are unchanged and identical to STARsolo's.

E2's comparisons are re-run and unchanged: every tag configuration (CB_samTagOut included) is identical to STARsolo. Non-solo yeast SE/PE output is byte-identical to main.

--soloVelocytoAmbiguous removed

A rustar-only option, from #90 via #133, that folded ambiguous molecules into spliced. It rested on the old rule putting exon-only reads in "ambiguous"; under STAR's model they are spliced, and "ambiguous" means exonic in one isoform and intronic in another. It had no STAR equivalent and no DIVERGENCE entry, so it's gone; anyone who wants the old merge can add spliced and ambiguous downstream.

Stats files as STAR writes them

Missing before this PR:

• Barcodes.stats;
• each feature's Features.stats and UMIperCellSorted.txt;
• SJ's and Velocyto's Summary.csv.

The Gene/GeneFull Summary.csv was there, but it lacked the Q30 and "Unique+Multiple" lines, printed six fixed decimals, and computed several values its own way.

What's ported:

• Barcodes.stats is STAR's SoloReadBarcodeStats. Each read counts once under its final cbMatch code, with a UMI rejection replacing the barcode's code (addStats, SoloReadBarcode_getCBandUMI.cpp:93-130); codes STAR doesn't list fall to its default. With --soloCBwhitelist None STAR bins nothing (:95-96) and the file is all zeros; it's now written in that mode too, though our no-whitelist matrices are still unimplemented. A missing adapter or UMI is counted directly, as STAR counts it (:340-350). Extraction now reports why there's no barcode instead of a bare None.

• Features.stats is STAR's SoloReadFeatureStats, from its three stages:

  • per read (SoloReadFeature::record): unmapped, no feature, multi-feature;
  • at counting (inputRecords): exact, oneExact and 1MM_multi posterior rejections. For SJ, these count only records whose junction is a matrix row, as STAR skips the others first (SoloReadFeature_inputRecords.cpp:23-27);
  • at UMI collapse: reads, unique reads, cells, UMIs.

  Velocyto keeps only the UMIs and cells: STAR's countVelocyto never sums the per-read counts into the feature (only countCBgeneUMI calls addStats).

• Summary.csv is STAR's statsOutput. Each line's definition and arithmetic follow STAR's:

  Line	Ported as
  Reads With Valid Barcodes	subtracts both barcode and feature rejections
  Q30 lines	barcode-read and clipped cDNA qualities, at ≥ Phred+33+30
  Reads Mapped to …: Unique / Unique+Multiple	counted reads, or NoMulti
  means	integer division
  medians	the element at nCells/2
  UMI median	taken from every barcode's UMIs sorted descending
  doubles	C++'s default ostream format (0.0357143, 1, -inf)

• UMIperCellSorted.txt: written when cells are called.

• CB_samTagOut writes Barcodes.stats, as STAR does.

• CB_UMI_Complex: two kinds of read that both tools reject can land in different Barcodes.stats bins, because our Complex matching uses the product of the segment whitelists rather than STAR's per-segment matching. Recorded as DIVERGENCE §3.4; nothing else changes.

More STAR differences fixed along the way:

• Exact-match counts. A barcode's exact-match count, which oneExact and the 1MM_multi posterior use, counted reads whose UMI was then rejected. STAR counts only reads whose final code is 0 (the MATCH_CellRanger build counts the others too; the default build doesn't).
• UniqueAndMult-*.mtx wrote 0.50000 where STAR writes 0.5. E1's comparison was numeric and didn't catch it.
• 1MM_multi barcodes on SJ records were dropped; STAR resolves them by the posterior, as for Gene, and now so do we. A read with no junction is SJ noNoFeature whatever its barcode match; a 1MM_multi read with junctions was counted there too.
• --soloCellFilter TopCells N took the N barcodes with the most UMIs. STAR keeps every barcode with at least the UMIs of the one at index N, so ties are kept (SoloFeature_cellFiltering.cpp:46-59). With TopCells 100 on the fixture STAR calls 106 cells and we called 100; now identical. This was already wrong on main.

Validation: STAR's own recorded command line was replayed for 30 STARsolo configurations of the fixture. They cover all barcode match types, every multimapper method, both dedups, all strands, EmptyDrops, MultiGeneUMI_CR, CellRanger4, two CB_UMI_Complex adapter setups, CB_samTagOut, the sorted-BAM tags, Velocyto, and paired-end.

453 of 457 Solo.out files are byte-identical to STARsolo's. The other four are in CB_UMI_Complex, where DIVERGENCE §1.6 rejects homopolymer UMIs of every base and STAR only poly-A. 43 reads land in noUMIhomopolymer instead of STAR's bins; the matrices are unchanged. §1.6 now says so.

Seven more configurations, all byte-identical (112 files):

• --soloFeatures Gene SJ with 1MM_multi, and with 1MM_multi_pseudocounts;
• --soloCellFilter TopCells 100, and None;
• Velocyto with --outSAMtype None;
• the synthetic fixture of test_starsolo_stats_files (below), run through STAR: all 14 Solo.out files;
• --soloCBwhitelist None --soloCBmatchWLtype Exact: Barcodes.stats, the only file we write in that mode.

The fixture has only 39 SJ reads, so the synthetic test is what exercises SJ's 1MM_multi path.

Yeast SE/PE alignments and SJ.out.tab are byte-identical to main.

CellRanger.summary.csv is opt-in

A rustar-only extension, from #90 via #133. It was written whenever Gene and GeneFull both ran, had no STAR equivalent and no DIVERGENCE entry. It's now written only with the new --soloCellRangerSummary Yes, and is documented in DIVERGENCE §1.8 and the CLI reference. By default Solo.out holds only STARsolo's files; a test checks both ways.

Not in this PR (all already the case on main)

• SmartSeq gets none of these files. A SmartSeq run of the chr19 reads against STAR showed more: STAR also writes GeneFull, SJ and filtered/, which we skip, and our Gene counts differ (same 315 entries, different values). That path needs a port of its own.
• --soloCBwhitelist None: no matrices, and it requires --soloCBmatchWLtype Exact where STAR accepts the default.
• sM and the CR/UR tags for CB_UMI_Complex reads with no adapter, UMI or CB.
• 1MM_multi barcodes on multi-gene records are dropped; STAR resolves them by the posterior. No fixture difference showed.

Tests

• Integration:
  • Velocyto, rewritten to STAR's semantics: against a single transcript, a read inside an exon is spliced, not ambiguous, and every matrix gets the gene's row, zeros included. Raw and filtered/ are both checked, once with Gene alone and once with GeneFull first and CB/UB, so the tags and Velocyto use different readInfos;
  • with a retained-intron isoform, an intronic read is ambiguous;
  • test_starsolo_stats_files: exact Barcodes.stats and SJ Features.stats contents for a fixture with a 1MM_multi read, an exact barcode with an N in its UMI, a barcode off the whitelist, and a read with no junction. If the bad-UMI read counted as an exact match, the posterior would drop the multi read. STAR gives the same files;
  • soloCellFiltering with TopCells now expects STAR's threshold.
• Unit:
  • the molecule collapse: intersection with OR-ed types, an emptied molecule, a multi-gene drop, undefined readInfo, and the detected-cell bitset;
  • the per-read strand filter; a spliced block entering one isoform's intron; a mate pair, whose gap is not a junction;
  • TopCells ties, and N past the end;
  • the two Velocyto refusals, and --soloCellRangerSummary validation.
• 808 tests in all. fmt and stable clippy are clean. Nightly clippy reports three lints in src/io/sam.rs and src/chimeric/bam.rs, files this PR doesn't touch; they're on main too.