solo: --soloFeatures Transcript3p, with --soloClusterCBfile by BenjaminDEMAILLE · Pull Request #158 · scverse/rustar-aligner

BenjaminDEMAILLE · 2026-07-29T08:28:06Z

--soloFeatures Transcript3p: quantify transcripts rather than genes, from where each read's 3' end sits relative to each transcript's.

What changed

In a 3'-biased assay that distance is what separates isoforms: a read 200 bases from the end of one and 4000 from the end of another is evidence for the first. The feature records, per read, every transcript the alignment is concordant with and the spliced distance to that transcript's 3' end. The distribution of those distances is estimated from the run's own histogram, smoothed with a running average, cut where the 3' peak decays into the body, and used as the likelihood in an EM over UMIs.

--soloClusterCBfile (new) assigns cells to clusters.

Output under Solo.out/Transcript3p/raw/: matrix.mtx (transcripts × clusters), features.tsv, and transcriptEndDistanceDistribution.txt.

Two behaviours that are not the obvious ones

Output is per cluster, not per cell, and --soloClusterCBfile is required. One cell does not have enough UMIs to resolve isoforms, so the EM would be fitting noise. Asking for the feature without a clustering is refused rather than run.

A UMI seen on several reads contributes the intersection of their transcript sets, not the union. Those reads came from one molecule, so a transcript missing from any of them cannot be its source. The union would let one stray read resurrect an isoform every other read excluded.

What this reuses rather than reimplements

Concordance needed no new code. align_to_transcripts already refuses to project an alignment that leaves the transcript, touches an intron, or crosses a junction the transcript does not have, which is what STAR's Transcriptome_classifyAlign.cpp calls Concordant. A projection that survives is concordant; one that does not, is not. The projection also puts the 5' end at coordinate zero for both strands, so the distance to the 3' end is one expression rather than two.

Faithful to STAR including its quirks

Two are reproduced rather than corrected, because the cut point and every weight depend on them (SoloFeature_quantTranscript.cpp):

the running-average divisor is min(2N+1, i + N), not the number of elements actually summed;
the per-transcript length factor reads the cumulative distribution at trLen - 1.

Numbers are formatted as C++'s default stream prints them (six significant digits, fixed inside [1e-4, 1e6) and scientific outside) because the normalised distribution runs down to ~1e-4, where Rust's {} and C++'s default disagree on both notation and digit count.

Verification

Nine unit tests in src/solo/transcript3p.rs:

the_distance_distribution_is_normalised_and_cut_at_the_peak
short_transcripts_get_a_length_correction, and a long one gets none
a_unique_umi_goes_entirely_to_its_transcript
an_ambiguous_umi_follows_the_evidence: nine unambiguous UMIs pull the tenth
a_umi_seen_twice_keeps_only_the_shared_transcripts: the intersection rule
a_cell_outside_every_cluster_is_skipped
cluster_file_parsing_skips_unknown_barcodes, a_trailing_barcode_without_a_cluster_ends_the_parse
numbers_are_formatted_the_way_c_plus_plus_prints_them

Gate: 569 lib + 26 integration tests, cargo clippy --all-targets -- -D warnings, cargo fmt --check, MSRV 1.89, all green.

Output-neutral unless --soloFeatures Transcript3p is requested; it adds a feature directory and changes nothing existing.

Two corrections to the previous revision of this description

"That needs a STARsolo oracle I do not have available" was wrong. STARsolo is the same binary as STAR, and STAR 2.7.11b was available the whole time. I tried to run the differential rather than repeat the claim, and it does not run, for a reason worth stating precisely rather than as "no oracle".

STAR marks this feature as under development. In parametersDefault, and in STAR --help, both Transcript3p and --soloClusterCBfile sit between:

#####UnderDevelopment_begin : not supported - do not use
                            Transcript3p    ... quantification of transcript for 3' protocols
#####UnderDevelopment_end

The port follows STAR's code, so it inherits whatever that code does, including its unfinished parts. This is not a stable interface upstream, and a reviewer should weigh that before taking it. The module doc now says so; it did not before.

Why the differential does not run

Transcriptome.cpp:17 picks the directory it loads gene and transcript info from:

trInfoDir = P.pGe.sjdbGTFfile=="-" ? P.pGe.gDir : P.sjdbInsert.outDir;

An index built with --sjdbGTFfile records that path in genomeParameters.txt, and STAR reads it back, so sjdbGTFfile != "-" on any later run against that index and trInfoDir becomes P.sjdbInsert.outDir. That is empty when the constructor runs, so STAR opens /geneInfo.tab and aborts:

Transcriptome.cpp:18:Transcriptome: exiting because of *INPUT FILE* error: could not open input file /geneInfo.tab

Reproduced on a fixture built for this: 20 cells, 5000 reads drawn 3'-biased from 200 yeast transcripts, --soloType CB_UMI_Simple. It fails the same way on --soloFeatures Gene alone, so it is not specific to Transcript3p, and passing --sjdbGTFfile - does not help because genomeParameters.txt is applied after the command line. With a GTF-less index plus --sjdbGTFfile at mapping, STAR does write _STARgenome/geneInfo.tab, but the Transcriptome constructor runs before sjdbInsert.outDir is set and aborts anyway.

So the oracle exists but I could not get a Solo.out out of it on this machine and this fixture. I would rather report that with the source line and the command than call it "no oracle available".

Still not run, then: the new matrix has not been compared against STARsolo's. The unit tests below pin the arithmetic and the output format; they do not pin agreement with STAR.

Split out of #152 following the one-theme rule in CONTRIBUTING.md.

Quantifies transcripts rather than genes, from where each read's 3' end sits relative to each transcript's. In a 3'-biased assay that distance is what separates isoforms: a read 200 bases from the end of one and 4000 from the end of another is evidence for the first. The distribution of those distances is estimated from the run's own histogram, smoothed and cut where the 3' peak decays into the body, and used as the likelihood in an EM over UMIs. Concordance needed no new code. `align_to_transcripts` already refuses to project an alignment that leaves the transcript, touches an intron, or crosses a junction the transcript does not have — which is exactly STAR's `Concordant` (`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant; one that does not, is not. The projection also puts the 5' end at coordinate zero for both strands, so the distance to the 3' end is one expression rather than two. Two behaviours worth stating because they are not the obvious ones: Output is per cluster, not per cell, and `--soloClusterCBfile` is required. A single cell does not have enough UMIs to resolve isoforms, so the EM would be fitting noise. Asking for the feature without a clustering is refused rather than run. A UMI seen on several reads contributes the *intersection* of their transcript sets. Those reads came from one molecule, so a transcript missing from any of them cannot be its source. Taking the union would let a single stray read resurrect an isoform every other read excluded. Two of STAR's quirks are reproduced rather than corrected, because the cut point and every weight depend on them: the running-average divisor is `min(2N+1, i + N)` rather than the number of elements actually summed, and the transcript length factor is taken from the cumulative distribution at `trLen - 1` (`SoloFeature_quantTranscript.cpp`). Numbers are formatted the way C++'s default stream prints them — six significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s default disagree on both notation and digit count.

CONTRIBUTING.md requires the description to match the code; the entries for the other themes split out of scverse#152 belong to their own PRs.

Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it.

parametersDefault puts both Transcript3p and --soloClusterCBfile between "#####UnderDevelopment_begin : not supported - do not use" and "#####UnderDevelopment_end", and STAR --help prints that banner around them. The module said none of this. It matters for how the port is read: it follows STAR's code, so it inherits the unfinished parts of that code, and a differential against STAR compares two implementations of something STAR does not support. A reviewer should be told that before deciding to take it.

* fix(solo): MultiGeneUMI_CR gives a tied UMI to nobody, not to everybody `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest read count. CellRanger's rule is the opposite on exactly that case: the gene with the *strictly* highest count takes the UMI, and a tie means no gene counts it. STAR walks the genes keeping a running maximum and clears its winner whenever it meets an equal count (`SoloFeature_collapseUMIall.cpp:212-224`): if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; } else if (ig.second==maxu) { maxg=-1; }; ... if ( maxg+1==0 ) continue; // not counted for any gene One read per gene is the ordinary shape of a multi-gene UMI, and it is always a tie, so the old rule made the flag inert in practice rather than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared between two genes), against STAR 2.7.11b with the same flags: identical entries STAR counts rustar counts before 13 749 / 14 806 15 423 16 465 after 13 902 / 13 967 15 423 15 414 The flag removed nothing at all before; STAR removes 1 030 counts. The gap goes from +1 042 to -9. The outcome does not depend on the order the genes are visited — a strict maximum always ends as the winner, a tie always ends with none — so iterating a `HashMap` here stays deterministic. `multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests missed: they only covered 3 reads against 1, where both rules agree. Not yet implemented, and stated so rather than left to be discovered: STAR applies a second condition, that the winning gene must also hold the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file, lines 226-232). That needs the pre-correction counts, which this code does not keep. The 65 entries still differing out of 13 967 are the place to look for its effect. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * docs(changelog): record the MultiGeneUMI_CR tie fix Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs STAR corrects UMIs within each gene *before* deciding which gene owns a UMI, and applies two conditions, not one (`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`): 1. one gene must hold a strictly higher read count than every other, on the **corrected** UMI map — that is #173, already landed; 2. and that winner must not be beaten in the **uncorrected** map at the same key. The second condition exists because correction moves reads between UMIs: a gene can win only because correction folded a neighbouring UMI onto it, and STAR rejects that win rather than counting it. Reproducing it needs the order STAR uses. The generic path here filters multi-gene UMIs first and corrects afterwards, which cannot express either condition: by the time correction happens the ownership decision is already made. `MultiGeneUMI_CR` therefore takes its own path, which is also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`, so there is no combination this bypasses. `cellranger_1mm_map` exposes the correction mapping that `cellranger_1mm` already computed and threw away. Measured against **CellRanger 10.0.0** on the 20 000-read fixture from #172, with #165 and #173 also applied: identical entries CellRanger rustar #165 + #173 13 651 / 13 709 15 111 15 091 plus this change 13 676 / 13 709 15 111 15 116 Entries CellRanger has and we do not go from 29 to 7, and the count gap from -20 to +5, which is 0.03%. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * docs(divergence): drop a reference to a test that no longer exists * feat(solo): bit-exact libc++ mt19937, generate_canonical and discrete_distribution STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to doubles with `std::generate_canonical<double, 53>`, and picks categories with `std::discrete_distribution`. Two of those three are implementation-defined in the parts that matter: the standard fixes mt19937's output but not how `generate_canonical` consumes it, and says nothing about how `discrete_distribution` maps a uniform onto categories. So porting "the algorithm" is not enough — it has to be libc++'s algorithm, because that is what STAR is built against and where its numbers come from. libc++ accumulates two 32-bit draws in *ascending* significance and divides by 2^64; a most-significant-first accumulation, or one draw scaled to 53 bits, both give perfectly good uniforms and neither reproduces STAR. Every expected value in the tests came out of a C++ program compiled against the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp` is that program, kept so the values can be regenerated rather than trusted. `generate_canonical` is compared as bit patterns, since a difference in the last place changes which category a sample lands in. Not yet wired into the EmptyDrops path. `solo::count` samples with a `SplitMix64` stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree on an arbitrary number of draws. It is true of whatever was checked and unknown elsewhere. Replacing it moves cell calls, so it belongs in its own change with the solo differential run against it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * feat(solo): EmptyDrops_CR uses Simple Good-Turing and libc++'s sampler Two approximations in the CellRanger cell-calling path are replaced by what CellRanger and STAR actually compute. Both move cell calls, which is the point: the previous numbers were plausible rather than right. The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson, via Elworthy's implementation, which is what STAR vendors). The ambient counts come from a small sample of empty droplets, so a gene seen twice there is not twice as likely as one seen once, and a gene seen zero times is not impossible — it is one the sample was too small to show. SGT fits the frequency spectrum and reserves mass for the unseen from the singleton rate, then smooths the rest along a log-log line. What was here before had the right shape and the wrong numbers: it reserved mass the same way but distributed the remainder in proportion to raw counts, with no smoothing at all. The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and `std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as STAR seeds it. The previous sampler was a SplitMix64 stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree over an arbitrary number of draws. The libc++ types were ported and checked against real libc++ in the previous commit on this branch; this wires them in. One generator per simulation, no shared state, so the walks still run in any order on any number of threads and give the same p-values. D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer than five distinct frequencies and `analyse()` bails, so it reads whatever the stack held. Here it is zero from construction, which is what "no basis for reserving unseen mass" means. Recorded in docs-old/dev/divergences.md. * docs: record the EmptyDrops SGT divergence in DIVERGENCE.md Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact / Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier version of this work carried. * docs(divergence): file the EmptyDrops entry under section 1, note the second RNG * fix(params): refuse MultiGeneUMI_CR without --soloUMIdedup 1MM_CR * feat(solo): --soloFeatures Transcript3p, with --soloClusterCBfile Quantifies transcripts rather than genes, from where each read's 3' end sits relative to each transcript's. In a 3'-biased assay that distance is what separates isoforms: a read 200 bases from the end of one and 4000 from the end of another is evidence for the first. The distribution of those distances is estimated from the run's own histogram, smoothed and cut where the 3' peak decays into the body, and used as the likelihood in an EM over UMIs. Concordance needed no new code. `align_to_transcripts` already refuses to project an alignment that leaves the transcript, touches an intron, or crosses a junction the transcript does not have — which is exactly STAR's `Concordant` (`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant; one that does not, is not. The projection also puts the 5' end at coordinate zero for both strands, so the distance to the 3' end is one expression rather than two. Two behaviours worth stating because they are not the obvious ones: Output is per cluster, not per cell, and `--soloClusterCBfile` is required. A single cell does not have enough UMIs to resolve isoforms, so the EM would be fitting noise. Asking for the feature without a clustering is refused rather than run. A UMI seen on several reads contributes the *intersection* of their transcript sets. Those reads came from one molecule, so a transcript missing from any of them cannot be its source. Taking the union would let a single stray read resurrect an isoform every other read excluded. Two of STAR's quirks are reproduced rather than corrected, because the cut point and every weight depend on them: the running-average divisor is `min(2N+1, i + N)` rather than the number of elements actually summed, and the transcript length factor is taken from the cumulative distribution at `trLen - 1` (`SoloFeature_quantTranscript.cpp`). Numbers are formatted the way C++'s default stream prints them — six significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s default disagree on both notation and digit count. * refactor(solo): drop Transcript3pAcc::merge, which nothing calls Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it. * docs(solo): note that STAR marks Transcript3p under development parametersDefault puts both Transcript3p and --soloClusterCBfile between "#####UnderDevelopment_begin : not supported - do not use" and "#####UnderDevelopment_end", and STAR --help prints that banner around them. The module said none of this. It matters for how the port is read: it follows STAR's code, so it inherits the unfinished parts of that code, and a differential against STAR compares two implementations of something STAR does not support. A reviewer should be told that before deciding to take it. * fix(solo): implement MultiGeneUMI_All instead of aliasing it to MultiGeneUMI `--soloUMIfiltering MultiGeneUMI_All` resolved to the same variant as `MultiGeneUMI`, which is neither what STAR does nor what the option is documented to do. Of the three available behaviours it was the only one nobody had asked for. In STAR the option is a no-op: it is parsed and stored, but its consumption site tests only the `MultiGeneUMI` flag, so selecting it leaves the filter entirely off. Documented, it removes a UMI seen in more than one gene from *all* of them, rather than from the losers only. `UmiFiltering::MultiGeneUmiAll` now exists and does the documented thing: a UMI appearing in several genes is evidence of a collision or of chimeric amplification, so it is discarded outright rather than attributed to whichever gene happened to read deepest. Single-gene UMIs are untouched, which the test checks across every mode. Raised upstream as #144 before changing it, since "be faithful to STAR" and "do what the flag says" genuinely point in opposite directions here. Also adds `docs-old/dev/divergences.md`, recording this and the homopolymer-UMI rule, so deliberate differences are written down rather than rediscovered as surprises in a differential run. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * docs(changelog): record the MultiGeneUMI_All fix Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * docs(divergence): correct the MultiGeneUMI_All entry, defer the homopolymer one --------- Co-authored-by: Benjamin Demaille <benjamin.demaille@icloud.com> Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

Psy-Fer · 2026-08-06T11:13:05Z

Closed by #212

BenjaminDEMAILLE added 2 commits July 29, 2026 10:27

docs(changelog): keep only this PR's entry

9d1f787

CONTRIBUTING.md requires the description to match the code; the entries for the other themes split out of scverse#152 belong to their own PRs.

BenjaminDEMAILLE mentioned this pull request Jul 29, 2026

solo: --soloUMIfiltering MultiGeneUMI_All (closes #144) #152

Closed

BenjaminDEMAILLE added 2 commits July 29, 2026 11:53

refactor(solo): drop Transcript3pAcc::merge, which nothing calls

295968a

Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it.

This was referenced Jul 30, 2026

solo: two CellRanger-matching flags are inert, and the count matrix is 8.9% off #172

Open

PR review queue: what is ready, what is stacked, what is waiting on a decision #195

Open

Psy-Fer mentioned this pull request Aug 6, 2026

Batch PRs #173, #175, #156, #158, #152 #212

Merged

Psy-Fer closed this Aug 6, 2026

Navigation Menu

Uh oh!

solo: --soloFeatures Transcript3p, with --soloClusterCBfile - #158

solo: --soloFeatures Transcript3p, with --soloClusterCBfile#158
BenjaminDEMAILLE wants to merge 4 commits into
scverse:mainfrom
BenjaminDEMAILLE:bd/solo-transcript3p

BenjaminDEMAILLE commented Jul 29, 2026 •

edited

Loading

Uh oh!

Psy-Fer commented Aug 6, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

BenjaminDEMAILLE commented Jul 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What changed

Two behaviours that are not the obvious ones

What this reuses rather than reimplements

Faithful to STAR including its quirks

Verification

Two corrections to the previous revision of this description

Why the differential does not run

Uh oh!

Psy-Fer commented Aug 6, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

BenjaminDEMAILLE commented Jul 29, 2026 •

edited

Loading