How the STARsolo pipeline is parameterised, and why

[lhqing]

Working notes from the session that produced #198. The immediate trigger was narrow — "our CRAM has no barcode" — but answering it forced us to justify every STARsolo flag we pass, and to reject two we were tempted by. Recording the reasoning here because the why behind a parameter is the part that evaporates, and next time someone will be staring at the same shell block wondering whether they can change it.

Everything below was measured on one real lane of GSE208154 / SAMN29720279 (L001, 12,993,522 reads, ce11 + WS298, in the pinned align-rna image), not estimated. Where a number is extrapolated to a full sample it is marked as such; the extrapolation was validated at 0.3% error against the real 4.51 GB deliverable.

---

1. The frame: three owners, and which one owns a flag

This is the single most useful thing to internalise before touching starsolo.smk. Every parameter we pass STAR belongs to exactly one of three owners, and the answer determines where you write it:

owner	owns	lives in	gated by
KB (kb/specs/<tech>/spec.yaml)	how to parse reads — byte-decided facts	backend.params → config["solo"]	params gate: KB keys verbatim from spec
processing manifest	what to count — instructable choices	processing.yaml → config	params gate: processing keys verbatim from manifest
the module (starsolo.smk)	output format and alignment quality — neither parse nor count	a literal in the shell block	nothing; it is an implementation detail

compose/params.py enforces that the emitted key set is exactly union(KB keys, processing keys) — disjoint, no orphans. And required_config is computed from module source (workflows/__init__.py), not typed beside it. Those two facts combine into a rule with teeth:

Writing SOLO["soloFoo"] in the shell block silently makes solo.soloFoo a required config key that every KB spec must then declare. Writing --soloFoo bar as a literal costs nothing.

So the question "should this be configurable?" has a real price attached, and the default answer for a quality/format knob is no — hardcode it. But "costs nothing" is about effort, not correctness: if the value genuinely varies by chemistry, a literal (or a branch) is not cheap, it is simply wrong, and wrong in a way that stays invisible until a new chemistry arrives. §2b is that case, caught just in time. The precedent was already in the tree: --outSAMtype used to be a KB param in star.smk and a literal in starsolo.smk. Two modules, two answers, one detail. It is now a literal in both.

2. The parameters, and why each is what it is

Parsing — KB-owned, byte-decided

flag	source	why it is not a choice
--soloType	spec	The chemistry is this. CB_UMI_Simple (10x) vs CB_UMI_Complex (SPLiT-seq, BD Rhapsody) is a fact about the molecule.
--soloCBstart/CBlen/UMIstart/UMIlen	spec (Simple)	Contiguous barcode → offsets locate it.
--soloCBposition/UMIposition	derived from element model (Complex)	Barcodes scattered between linkers → no start/length exists to give. Derived at compose from the spec's element coordinates rather than declared, so the two cannot disagree.
--soloCBwhitelist	spec	Materialised on demand by rule onlist as a temp() — 10x v3 is 6,794,880 barcodes / 111 MB, and compiling one dataset three ways used to write it three times, forever.
--soloStrand	spec	Byte-decided. Also the one thing the params gate cannot check — only a kb e2e count-matrix run catches an inverted strand.
--soloBarcodeReadLength	spec, optional	10x R1 is routinely sequenced longer than the 28 nt the barcode occupies; STAR's default check then FATALs. Read with SOLO.get(...), deliberately — a subscript would make it required for SPLiT-seq too, whose spec must not declare it.
--soloCBmatchWLtype	spec	How hard to try to rescue a barcode that isn't on the list. The value tracks each vendor's own reference pipeline, so it is a per-chemistry fact. See §2b — this one moved out of the module, and the reasoning generalises.

2b. The flag that changed owners, and the test that caught it

--soloCBmatchWLtype lived in the module as a soloType branch: 1MM_multi_Nbase_pseudocounts for Simple, 1MM for Complex. That looked fine and it was wrong, for a reason worth stealing:

chemistry	soloType	value	tracks
10x (7 specs)	Simple	1MM_multi_Nbase_pseudocounts	CellRanger ≥3
BD Rhapsody (3), original SPLiT-seq (1)	Complex	1MM	#43 / STAR #1607
Parse Evercode (planned)	Complex	EditDist_2	Parse's own pipeline

Parse Evercode and BD Rhapsody are both CB_UMI_Complex and need different values. A yes/no branch yields two answers; we need three. So the branch was not merely inelegant — it was incapable, and would have failed silently the day someone added Parse (they'd hit the wall, not understand it, and settle for 1MM).

The test: can the discriminator you're branching on distinguish every case you'll ever need? If two chemistries agree on the branch variable but disagree on the value, the branch is the wrong mechanism and the fact belongs with the chemistry.

It is now a backend.params key in all 11 starsolo specs. The price is real and worth naming: a new SOLO["..."] subscript becomes a required config key (required_config is computed from module source), so every spec must declare it or compose refuses.

Legality is itself chemistry-dependent, measured by running STAR rather than reading --help:

value	Simple	Complex
Exact, 1MM	legal	legal
1MM_multi, 1MM_multi_pseudocounts, 1MM_multi_Nbase_pseudocounts	legal	rejected
EditDist_2	rejected	legal

An illegal pair is a hard STAR FATAL on a compute node, after the genome loads. That is now a params-gate assertion instead — compose refuses, which is that gate's whole purpose. Note STAR's global default 1MM_multi is illegal for Complex, so the old 1MM pin was load-bearing rather than stylistic.

⚠️ Every wrong answer in this area passes a 10x-only test suite. The failures land on Complex chemistries exclusively.

Counting — processing-owned, instructable

flag	value	why
--soloFeatures	Gene GeneFull GeneFull_ExonOverIntron GeneFull_Ex50pAS Velocyto	The "compute every alternative in one pass" law: STAR produces all five for roughly one pass of work and the outputs are small, so we compute them all and let the consumer choose.

Output format and quality — module-owned literals

flag	value	why
--outSAMtype	BAM SortedByCoordinate	Forced. STAR: "CB and/or UB attributes can only be output in the sorted BAM." Happy side effect: seqforge io cram no longer needs samtools sort, which deletes a whole re-sort pass and a temp-file leak (§4).
--limitBAMsortRAM	from config["mem_mb"]	Default 0 means "reuse the genome allocation", which is too small on a tiny fixture and FATALs. Not optional once you sort.
--outSAMattributes	NH HI AS nM CB UB	The whole point of #198 — see §3.
--clipAdapterType	CellRanger4	CellRanger ≥4 polyA/TSO clipping.
--outFilterScoreMin	30	Drops low-score alignments.
--soloUMIfiltering	MultiGeneUMI_CR	Removes UMIs shared across genes.
--soloUMIdedup	1MM_CR	CellRanger's UMI collapsing (STAR's default is 1MM_All).
--soloCellFilter	EmptyDrops_CR	§5.
--soloMultiMappers	Unique (default, but now a decision)	§6.
--readFilesIn	cDNA mate first, barcode second; each mate's runs comma-joined	A sample pooled across N runs maps in one pass. Joining with spaces instead — the old bug — makes STAR read the extra files as extra mates and crash.

The five CellRanger-flavoured flags come from scRecounter (Arc Institute, same problem at 10⁵ scale) and are the documented "CellRanger ≥4 equivalent" set from Kaminow, Yunusov & Dobin 2021. We were emitting STARsolo defaults, which are not CellRanger-comparable — a real problem for a corpus meant to be compared against published matrices.

3. The CRAM: what we were shipping, and what we now ship

We passed no --outSAMattributes at all, so STAR used its Standard default NH HI AS nM. STARsolo puts only the cDNA mate in the BAM, and the barcode lives solely in R1 — so the barcode was irrecoverably absent from 920 GiB of retained CRAMs. They could not be recounted under a new GTF, could not give transcript-level counts, could not feed alevin-fry. They were browsable and nothing else.

Tag costs, measured (added naively, nothing else changed):

added	Δ
+ CB UB	+28.0%
+ CR UR	+29.7%
+ CY UY	+9.7%
+ GX GN	+0.3%

The surprise was on the other side of the ledger:

lever	effect
short read names	−16.2%
primary alignments only	−17.8%
archive codecs (level=9,lzma,bzip2)	−3.6%, slower
--output-fmt-option lossy_names=1	0% — it does nothing. Measured.

38-char Illumina read names (K00125:217:HCL2YBBXY:8:2111:24637:43374) were 16% of the file, and they are meaningless once CB/UB exist — the name was only ever the join key back to R1. So:

Adding the barcode and dropping the now-redundant names makes the file 12% smaller than what we shipped before. 4.21 → 3.70 GiB per sample.

That is the whole lesson of §3: the naive framing was "how much do we pay for the barcode", and the answer was "nothing, you pay for what the barcode makes obsolete".

The road not taken, recorded because it is cheap to revisit

R1 in 10x v3 is exactly 28 nt = 16 (CB) + 12 (UMI). Nothing in it is discarded. So CR+UR+CY+UY, plus --outSAMunmapped Within, makes the CRAM a byte-exact replacement for both FASTQs. Verified rather than assumed — reconstructed R1 and R2 (sequence and quality, un-reverse-complementing reverse-strand records) for 200,000 reads and diffed against the source:

compared 200000 reads: EXACT=200000  MISMATCH=0

7.00 GiB vs 15.39 GiB of FASTQ — 45% — and it feeds anything via samtools fastq. We did not take it, because for accessioned data SRA already holds the FASTQ. For in-house data with no accession the calculus flips, and the measurements are already done.

4. The temp-file leak, and why the fix is structural

cram.py built samtools sort with no -T, so temp files landed in the CWD as samtools.<pid>.<tid>.tmp.NNNN.bam. A killed sort leaked them; snakemake could not clean them because they were undeclared outputs. 41.4 GiB across 5 pipeline dirs before we noticed.

The tempting fix is to pass -T into a snakemake-owned directory. The better fix is --outSAMtype BAM SortedByCoordinate, which removes samtools sort from the pipeline entirely — the leak becomes impossible rather than configured correctly.

Prefer the mechanism that cannot fail over the one that must be set up right. This is the same instinct as rule onlist using temp() instead of a documented cleanup step, and _STARtmp being rm -rf'd in the shell block rather than left to a convention.

5. --soloCellFilter: decided on reproducibility, not biology

The instinct is to argue about which cell caller is better. The decisive fact turned out to be architectural and is already written in our own code:

• h5ad.py:27 — "Everything is read from raw/, never filtered/. Cell calling is a downstream decision."
• qc.py:127-130 — the filtered/ tree is consumed only to record default_filtered_barcodes in the QC bundle.

So --soloCellFilter changes the QC bundle and not a single count in any .h5ad. Blast radius: one QC field. That reframes the question from "which is biologically right" to "which is cheap and reproducible".

The real risk was that EmptyDrops_CR is Monte-Carlo (simN = 10000), and a nondeterministic QC bundle is a problem for a content-addressed compiler. Tested with STAR --runMode soloCellFiltering (which re-filters an existing raw matrix in ~1 s — very handy):

run	cells	md5
CellRanger2.2 (old default)	441	f7439745760d9f4e
EmptyDrops_CR, 24 threads	486	23ac8197f0b75226
repeat, 24 threads	486	23ac8197f0b75226
4 threads	486	23ac8197f0b75226
--runRNGseed 12345	486	23ac8197f0b75226

Bit-identical across repeats, thread counts and seeds. EmptyDrops_CR is also a strict superset here: 45 barcodes rescued, 0 dropped.

Caveat kept deliberately: seed-independence is stronger than expected, and the benign explanation — no candidate sitting near the FDR boundary on this sample — cannot be excluded. So #198 carries an action to encode the check as a test rather than trust the observation. Cheap, because soloCellFiltering runs on an existing matrix.

6. --soloMultiMappers: the headline number was a lie

This is the most instructive one. scRecounter passes EM Uniform; we rejected it.

Multi-gene UMIs measured at 18.17% (Gene) / 19.56% (GeneFull) — high enough to look like we were throwing away a fifth of the data. Then the biotype breakdown:

biotype	multi-gene UMIs	% of all multi-gene	genes
rRNA	557,527	86.65%	3
protein_coding	70,137	10.90%	1,273
pseudogene	8,949	1.39%	213
ncRNA	5,314	0.83%	59

Excluding rRNA: 3.13%. The three genes are rrn-1.1, rrn-1.2, rrn-3.1 — the tandem rDNA array at I:15,062,083-15,071,033. rrn-1.1 and rrn-1.2 have 1 and 0 unique UMIs and EM hands each ~243,220: the copies are identical, EM has no information, so it splits evenly and emits a large arbitrary number that reads exactly like data. (The same array is why chromosome I holds 87.5M of the CRAM's 198.8M alignment records — a detail that looked like noise at the start of the session and turned out to be the whole story.)

Check the biotype composition before believing an aggregate percentage. A single repeat family can carry an entire headline statistic.

Also against: all four multimapper matrices are fractional (integer=False), which breaks DESeq2-style pseudobulk and NB likelihoods; and it costs 8 extra matrices (one per gene-axis feature per method — not Velocyto), each ~18% denser than the unique matrix, i.e. ~3.4× the matrix payload of every h5ad.

And the ecosystem agrees: multi-gene reads are discarded by default in Cell Ranger (MAPQ 255 only), STARsolo, alevin-fry (cr-like), and bustools count --genecounts. EM exists everywhere and is default nowhere. The legitimate exceptions are studies where multi-mapping regions are the biology — transposable elements (scTE, TEtranscripts), immune loci, recent paralogs.

The best part: we already had the diagnostic for free. Every shipped QC bundle carries features_stats.<feature>.MultiFeature and subMultiFeatureMultiGenomic under plain Unique — on the real sample, 25,037,969 and 24,463,852 (97.7% multi-genomic, the repeat signature). So we can identify a dataset with a genuine multimapping problem without paying for a single extra matrix. And matrix.mtx is byte-identical with and without the flag, so the decision is reversible per-dataset.

Before adding an output to answer a question, check whether the stats you already bundle answer it.

7. The unrelated 280 GiB

Chasing "could we just use a loose UMI cutoff instead of a cell caller" turned up something bigger than the CRAM work. STARsolo's raw/barcodes.tsv is the entire whitelist — 6,794,880 barcodes — and we write all of them into every h5ad. Only 845,694 (12.4%) have a nonzero count in any feature; obs alone is 108.7 MB of barcode strings indexing nothing.

	barcodes	size
as shipped	6,794,880	633.8 MB
drop all-zero barcodes (lossless)	845,694	95.5 MB
≥200 UMI (lossy — do not)	1,960	13.3 MB

−85%, provably lossless. ~280 GiB across the 350 h5ad files currently on disk.

Two traps found on the way:

1. The mask must be the union across all features, not X.sum() > 0. X is Gene (exonic) but GeneFull counts introns, so a barcode can be zero in Gene and nonzero in GeneFull. Filtering on X alone silently drops 22,319 barcodes.
2. Never push a UMI cutoff into the deliverable. At ≥200 UMI you discard 69% of UMIs — and those are precisely the low-count droplets that define the ambient profile. SoupX / CellBender / DecontX would have nothing to estimate from, and EmptyDrops could not be re-run downstream at all, since it builds its null from empty droplets. A permissive cutoff is a fine label; it is not a filter.

8. Reusable lessons

• Measure the thing, then extrapolate, then check the extrapolation. One lane × 16.567 predicted the real deliverable to 0.3%. That single check is what made every other number in STARsolo CRAM carries no barcode: add CB/UB (measured: 12% SMALLER than today), or stop writing it #198 quotable.
• Ask what a retained artifact is for. scRecounter's answer to the disk question is --outSAMtype None — keep no alignment at all. That is the honest baseline. A barcode-less CRAM was the one option nobody would defend.
• The cheapest byte is the one you stop storing. Read names (16% of a CRAM) and all-zero barcodes (85% of an h5ad) were both pure overhead nobody had priced.
• Check the biotype/category breakdown before believing a percentage. §6.
• Prefer the mechanism that cannot fail. §4.
• Host samtools 1.13 silently truncates CRAM 3.1. samtools ≥1.21 writes 3.1 by default, so anything reading our CRAMs must run inside align-rna. This bit us mid-session and produced a confusing "Failure to decode slice".
• Check that your discriminator can distinguish every case you'll need. A soloType branch gave two answers where three chemistries needed three — see §2b. If two cases agree on the branch variable but disagree on the answer, the branch is the wrong mechanism.
• Run the tool to learn its rules; don't read the help text. The soloCBmatchWLtype × soloType legality matrix came out of 12 actual STAR invocations. The help text was inconsistent (an option named EditDist_2 described as "edit distance of 3", plus a typo), and a comment in our own tree had recorded a restriction that was real but narrower than stated.
• STAR's error strings are documentation. Grepping the binary for EXITING because of surfaced both hard constraints on CB/UB and the soloCBmatchWLtype/Complex incompatibility, faster and more reliably than the manual.

---

Tracking issue: #198 (labelled ready-for-agent, three scoped work items).