solo: --soloFeatures Transcript3p, with --soloClusterCBfile

[BenjaminDEMAILLE]

--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.