~40% of every aging_SS3 worm is discarded as "unmapped: too short", and the recipe has no trimming step

[lhqing]

Roughly 40% of every worm in aging_SS3 is being thrown away at the aligner, and the pipeline
has no trimming step at all. This is a measurement, not yet a diagnosis — the primary-source review
of what to trim is still in flight — but the number is large enough and consistent enough to put
up now.

It surfaced incidentally: I was benchmarking --outSAMtype BAM SortedByCoordinate against
map-then-samtools sort (separate write-up to follow) and noticed the mapping speed was implausibly
low for a 1.3 GB ce11 index.

The measurement

Two independent sources, agreeing.

43 samples from the production run (ss3-ce11-ws298-ce321d3fc6d4, still in flight on fat02),
read straight out of each cell's Log.final.out:

metric	min	median	mean	max
input reads	372,963	10,195,271	11,102,148	25,580,704
uniquely mapped %	19.61	40.72	39.20	49.06
mapped to multiple loci %	5.06	20.35	19.42	23.81
unmapped: too short %	32.07	39.10	41.26	64.01
mapped to too many loci %	0.01	0.03	0.04	0.14
mismatch rate per base %	0.26	0.32	0.32	0.37

One median worm on an idle node (day9_N2_3, 12,590,188 pairs, cpu02, 8 threads, nothing else
running): 41.98% unique, 21.20% multi, 36.74% too short, mapping speed 27.74 M reads/hour.

The production medians are worse than the idle-node figure, which is what 20-way concurrency on one
node does to the speed column; the fractions are a property of the reads and match across both.

Why this is not a quality problem

Three columns rule out the boring explanations:

• Mismatch rate per base is 0.32%. The bases that do align, align cleanly. These are not bad reads.
• "Mapped to too many loci" is 0.03%. Repeat content is not eating them.
• The loss is in too short specifically — STAR's category for a read whose best alignment
  failed a length-relative filter, not for a read it could not place anywhere.

So ~40% of the reads carry enough clean sequence to be worth sequencing, and are being discarded on
a length test.

The leading hypothesis, labelled as one

STAR's --outFilterScoreMinOverLread and --outFilterMatchNminOverLread both default to 0.66
and are, as their names say, relative to the full read length — including any part of the read
that is not genomic. A read that is half adapter, TSO, polyA or oligo-dT cannot clear 66% of its
own length in matched bases no matter how good the genomic half is, and lands in too short.

aging_SS3 is exactly the shape that would produce this:

• The library structure resolved from the bytes is R1:ATTGCGCAATG@0:umi@11+8:GGG@19:cdna@22 — R1
  opens with 22 bp of non-genomic tag + UMI + motif.
• umi_extract trims R1 to offset 22 to lift the UMI out, so that prefix is handled.
• R2 is never touched at all, and neither read is checked for Tn5/Nextera read-through, polyA, or
  residual oligo-dT — SS3 is tagmented, so read-through into the adapter is expected on short
  fragments.

I want to be plain that the mechanism above is a hypothesis consistent with the numbers, not
something I have verified against the SS3 protocol or the STAR source yet. A primary-source review
(SS3 methods, zUMIs' actual behaviour, STAR's filter semantics, exact adapter sequences) is running
and I will append its conclusions here rather than guess at them now.

What it costs

Two different bills, and the second is the larger one:

1. Compute. STAR does full alignment work on every one of those reads before rejecting it.
   Some fraction of the plate's runtime is being spent to produce nothing.
2. Data. ~40% of a 784-worm sequencing run is not reaching the matrix. If even part of that is
   recoverable by trimming, it is recoverable for every cell, permanently, and it changes the
   per-cell depth every downstream analysis sees.

The second is why this is worth doing properly rather than quickly.

What I am not claiming

• Not that trimming will make STAR faster. Trimming shortens reads and removes failed alignment
  attempts, but it adds a pass over every FASTQ. Which way that nets out is a measurement nobody here
  has made, and I would not put it in a recipe on intuition.
• Not that 40% is anomalous for Smart-seq3. It may well be normal for untrimmed SS3 and I have no
  published baseline in hand to compare against. That is one of the things the review is for.
• Not that the current manifest is wrong. Nothing here is a compiler defect — seqforge
  resolved the chemistry and the read structure correctly from the bytes. This is a question about
  what the recipe should do, which is processing.yaml's side of the line.

Where this would land if it holds up

Trimming is a processing.yaml concern, not a manifest one — same data, different recipe, and the
dataset hash should not move because we decided to trim. Worth noting there is an ordering hazard
that any implementation has to respect: the UMI lives in the first 19 bp of R1, so anything that
trims R1 before umi_extract runs destroys it. Trimming has to happen after extraction, or be
taught the structure.

Related

• One recipe figure sizes both the genome residency and the BAM sort budget, and they scale with different things #349 — the memory figure. Same run, different knob; the benchmark that turned this up also
  produced the two numbers that issue says do not exist, and I will post those with it.
• First real in-house dataset: 784 C. elegans Smart-seq3 worms, compiled end to end #350 — the aging_SS3 compile this is measured against.

Numbers above are reproducible from Log.final.out in the pipeline's results/ and from
/share/lhqlab/wormbase/inhouse/aging_SS3/script/bench_star_sort/ on ircbc.

[lhqing]

Resolved: it was the Tn5 read-through, and clipping it recovers most of the loss

The hypothesis in this thread is confirmed, and the fix has landed (#356 → #357, measurement #358 → #359).

The mechanism was exactly as posted here. STAR's --outFilterScoreMinOverLread / --outFilterMatchNminOverLread are relative to the full read length; a clipped base leaves that length while a soft-clipped one does not, so a read half of which is adapter cannot clear 66% of itself however cleanly its genomic half aligns.

Controlled before/after, four published GSE207085 cells (mouse, 150 bp PE, mm10) — same uBAM from seqforge's own extractor, same STAR command, one flag different:

	uniquely mapped %	too short %	mismatch/base %
unclipped	42.39	54.43	0.57
clipped	63.59	22.75	0.54

+21 points of unique mapping; 58% of the too short category eliminated; mismatch rate flat, so this is not junk being waved through. Dose-dependent, too: the cell that lost more adapter (294 → 198 bp per fragment) gained more (+26.6) than the one that lost less (293 → 213, +23.2).

Your "not claiming 40% is anomalous for Smart-seq3" caution was right to make, and the answer is that it is not. The unclipped arm above — different lab, species and tissue — lands on 42.39% uniquely mapped against this thread's 42.0% median for aging_SS3. Two independent deposits, same pathology. Published comparison: Deserranno et al. (BMC Genomics 2026) report 34.03% for Smart-seq3xpress through zUMIs at stock settings.

Where it landed differs from what this thread proposed, and that is recorded rather than quietly reversed. You wrote that trimming is a processing.yaml concern and the dataset hash should not move. ADR-0048 takes the other view for this case and states the dividing line: the test is terminality, not the presence of trimming. Past the mosaic end a read carries adapter, index and flowcell primer — no recipe could rationally keep it — so it is a fact about the molecule and belongs to the chemistry entry. A quality floor or a length cut have real alternatives and remain recipe knobs. It bumps kb_version and re-keys run_id; dataset_hash does not move.

Your ordering hazard is respected and checked, not assumed. The clip rides the aligner, which reads the uBAM the extractor already wrote, so it cannot reach R1 before the UMI is lifted out. There is a test asserting the extractor's command line carries no clip flag.

One caveat worth carrying forward: you named three suspects — Tn5 read-through, polyA, and residual oligo-dT. Only the first is fixed. The fate budget shows 13.4% of the freed reads land in unmapped: other (stubs that were almost entirely adapter), which is a real ceiling; whatever the other two cost is still on the table and unmeasured.

Full tables and method: docs/research/smartseq3-tn5-read-through.md.