The bacterial fraction, a 14x age effect, and the check that refused to write it

[lhqing]

A 16-worm pilot of #425, run 2026-08-18 on ircbc under main @ 5231b69. Both arms compiled and
submitted; one finished, one refused. The refusal is the most useful thing in this report.

Before spending ~1600 mapping jobs on the full 784-worm plate, we ran 16 worms through both arms of
#425 — plain ce11 and the chimera ce11_ecHT115 — under the same code, on the same day. The plain
arm finished. The chimeric arm mapped every cell correctly and then refused to write its matrices,
on all 16 of 16 cells, for a reason nobody had been able to observe before because nobody had ever
run this split on real data.

Every number #425 asked for exists anyway. Here they are, then the bug.

The batch

day1_N2_{1..8} and day17_N2_{2,6,7,8,9,10,11,12} — one strain, two ages, eight replicates each.
The two ages were chosen to make the per-sample metadata path testable: a defect that fanned one age
across every sample would have been invisible on a batch whose samples all agree. It wasn't fanned —
the manifest came out age: {day1: 8, day17: 8}, strain: {N2: 16} — and the two ages turned out to
matter for a reason we did not anticipate. See the bacterial fraction below.

Both arms were built from the same 32 files, so dataset_hash came out byte-identical
(88194d4d…) with different run_ids. Same data, different recipe, exactly what the two-artifact
split promises. The chimera arm selected map/star-umi-chimera off the assembly name alone, with no
second flag: 69 jobs against the plain arm's 52, the difference being split_chimera ×16 and
umi_count ×2 instead of ×1.

1. Did the Tn5 clip engage? Yes — and this plate had never been observed post-clip

Average input read length is what STAR was left with after its own clipping, summed over both
mates. 150 bp paired-end with nothing clipped arrives as ~300.

	value
sequenced	~300 bp per fragment
observed, median of 16	224.5 bp (min 198, max 242)

The clip engaged, and this is the first time that has been seen on this plate — #360 fixed it on
published mouse cells and the worms were never re-run. Identical in both arms, as it must be:
clipping happens before alignment and the reference has no say in it.

This was the number that could have invalidated everything else. If it had sat at 300, every other
figure here would have been computed on a denominator full of adapter.

2. What stays unmappable: 12%, down from 38%

unmapped: too short	median
this plate, unclipped (archived)	38.25%
published unclipped Smart-seq3	38–54%
clipped mouse cells (#360)	22.75%
this pilot, plain arm	12.03% (5.10–18.80)
this pilot, chimeric arm	10.30% (5.01–15.48)

Uniquely mapped rose to 58.26% from the 42.0% this plate showed unclipped. The gain is larger
than #360's mouse cells got, which is the opposite of what that write-up predicted ("expect a smaller
absolute gain") — worth someone's attention, though the two datasets differ in species, tissue and
depth, so this is an observation and not a contradiction.

3. The bacterial fraction — and a 14× age effect nobody was looking for

Every figure below is a lower bound. A read ambiguous across worm and bacterium is dropped as a
multimapper rather than assigned to either, so the true share is this plus some part of §4.

share_ecHT115, of alignment records	median	range
all 16	0.769%	0.074–6.107%
day 1 (n=8)	0.144%	0.074–0.706%
day 17 (n=8)	2.081%	0.831–6.107%

Old worms carry roughly 14× the bacterial share of young ones, and the two groups do not overlap:
every day-17 worm is above every day-1 worm. That is consistent with what is known about aged C.
elegans accumulating live E. coli in the gut, and it is the kind of result the whole exercise
existed to make visible — for years this signal has been sitting in unmapped and too many loci
with nothing saying which.

Read this as a pilot, not a result. n=8 per group, no significance test was run and none should
be quoted from this. The pilot was designed to test a compiler, not to control for library depth or
batch, and the two age groups differ in sequencing depth. It is a strong enough signal to justify the
full 784 — which is exactly what a pilot is for.

4. How much the multimapper filter hides: almost nothing

#425 asked for this because it sets how much to trust §3. Comparing the two arms' multimapper shares
measured the same way on both — over alignment records from each arm's own CRAM:

increase (chimera − plain)	value
median	−0.13 pp
range	−1.22 to +0.18 pp

Essentially zero, and occasionally negative because the chimeric arm gains newly-mappable,
unambiguous bacterial records that dilute its own ratio. The lower bound in §3 is tight and the
per-component figures need no louder caveat than the one they already carry.

One methodological note worth keeping: the obvious version of this comparison — STAR's read-pair
percentages against the split's record counts — gives a median increase of +3.1 pp. That entire
number is a denominator artifact. STAR's percentages are over read pairs; the split counts alignment
records. Anyone repeating this must count both arms the same way.

5. What a chimeric index costs: ~4%, and the recipe is ~20× oversized

sacct was unavailable (slurmdbd refusing connections), so this was measured directly on one cell,
8 threads, --genomeLoad NoSharedMemory so the two arms are self-contained and comparable.

	plain ce11	chimeric ce11_ecHT115	delta
peak RSS, one mapping job	2.345 GB	2.430 GB	+3.6%
shared genome segment	1.157 GiB	1.214 GiB	+5.0%
index on disk	1.2 G	1.3 G	—

The chimeric twin's header says its memory figures are the base's, carried over identical and
honestly unmeasured. They now have a number, and the honest carry-over was the right call: a worm
plus ~4 Mb of bacterium costs about 4%.

The more interesting finding is the other direction. resources.mem_gb defaults to 48 GB and
measured peak RSS is ~2.4 GB — roughly 20× headroom. The measured cell (9.04M pairs) is not the
pilot's largest (65.6M, ~7.3×), so the true worst case is higher, but on the shape observed —
fixed genome cost plus a per-read buffer — still far under 48. That figure also sets
--limitBAMsortRAM to three quarters of itself, so it is not free to change and #425 was right to
call it its own ticket. It now has a measurement to argue from.

The bug: a correct refusal, from an assumption that was true and didn't imply what it seemed to

The chimeric arm mapped all 16 cells, wrote all 16 CRAMs, and then failed split_chimera on every
one:

these outputs kept a different number of first and second mates, as component: (read1, read2) —
{'ce11': (1030336, 1033258), 'ecHT115': (9569, 9944)}. Both mates of a template carry one NH and
sit on one chromosome, so a stateless filter keeps them together; a difference means that
assumption is false for this aligner, and the split is not the thing to trust

That check exists because split.py's own docstring says the two facts under it "were read off the
aligner's source and nobody has yet watched them hold on a real chimera." Nobody had. Now someone
has.

Both stated premises are true. Measured over ~13M kept records across two cells: zero templates
whose mates carry different NH, zero templates with mates on different components. The inference
drawn from them is what fails — neither premise says anything about whether the mate is in the file
at all.

Root cause: --outSAMunmapped is never set, anywhere in src/, so STAR uses its default None
and omits the unmapped mate of a pair where only one mate aligned. The split then keeps the surviving
mate — correctly, it is a mapped, unique, primary alignment — and there is no partner beside it.

The arithmetic closes exactly, with no residual, on both cells and both components: kept records
carrying flag 0x8 equal query names appearing exactly once.

cell	component	read1	read2	gap	singletons
day1_N2_7	ce11	1,030,336	1,033,258	2,922	33 + 2,955
day1_N2_7	ecHT115	9,569	9,944	375	0 + 375
day17_N2_2	ce11	4,978,307	4,995,009	16,702	886 + 17,588
day17_N2_2	ecHT115	145,608	146,881	1,273	36 + 1,309

The data is healthy and the check is too strict. The per-component BAMs were valid and complete
on disk when the assertion fired; snakemake then deleted them as possibly-corrupt. Good data, thrown
away by a false positive.

There is a sharper way to put it. The module already knew the first half: its header notes that of
the four drop categories, "only the first can occur under the flags the aligner runs with today" —
that is, it knew unmapped records never reach the split under --outSAMunmapped None. What did not
follow through is that a pair can therefore arrive half-present.

And it lands hardest exactly where the ticket cares. ecHT115's singleton rate runs 2.5–13×
ce11's, and the reason is measurable: bacterial fragments are shorter and much more 3′
soft-clipped (39.4 bp vs 20.0 bp on day1_N2_7; read1 aligned-fraction 0.49 vs 0.74). Read1 already
gives up ~20 bp to the tag/UMI/motif structure, so a short bacterial fragment runs out of margin
against STAR's 0.66 length filters and fails while its full-length read2 still clears — which is
precisely the observed "read2 present, read1 missing" asymmetry.

What the fix could be

Not applied — this is a decision with a real trade-off, and it belongs in a record.

• Relax the check to tolerate the known mechanism — verify the gap is fully explained by
  0x8-flagged records instead of asserting equality. Right in principle, but it needs per-template
  state, against a module whose stated design is "no name sort, no buffer, streams in constant
  memory."
• Downgrade it to informational — count singletons as their own category beside kept/dropped.
  Cheapest, preserves the streaming design, and keeps the fullest bacterial signal. Costs the safety
  net that would catch a genuinely different future problem, unless something else takes that job.
• Require both mates present — restores parity by construction, and discards real uniquely-mapped
  evidence, disproportionately from the rarer organism (~0.7–1.9% of ecHT115 against ~0.15–0.2% of
  ce11).

Worth doing regardless: set --outSAMunmapped Within KeepPairs so the split's unmapped bucket
reports reality instead of being structurally always-zero. On its own it does not fix the parity
check — the newly emitted mate lands in a component-blind drop bucket and the kept-side counts don't
move.

Two smaller things this run surfaced

seqforge report renders the broken arm as healthy. The chimeric report says
samples_finished: 16/16, alerts: [] while both headline deliverables are missing. Its completion
check is scoped to the mapping module and never asks whether rule all's artifacts exist. Read on
its own, that page is misleading.

A failed split strands the shared genome segment. --genomeLoad LoadAndKeep never reached its
unload step, leaving ~1.35 GiB resident on the node with nattch=0 after the job ended. Released by
hand. On a shared cluster that is somebody else's memory.

Where things are

pilot16/plain/seqforge/report.html                 # finished: combined.h5ad, 10.6 MiB
pilot16/chimera/seqforge/report.html               # mapped 16/16; no matrices (split refused)
pilot16/script/metrics/{per_cell.tsv,summary.json} # every number above
pilot16/script/diag/                               # the diagnosis, reproducible

all under /share/lhqlab/wormbase/inhouse/aging_SS3/.

Verdict on the full 784

Not yet. The plain arm would run today and produce a corpus-ready matrix. The chimeric arm would
burn ~800 mapping jobs, map every one of them correctly, and refuse to write a single matrix — the
same failure, 784 times.

Fix the parity check first. Everything else in #425 is measured, and the pilot cost two hours to
learn it.

[lhqing]

The read-fate ledger: every read pair, input to molecules

Requested by the data owner. Built from the finished artifacts — STAR's per-cell logs and the
pre-split chimeric CRAMs — and reconciled against STAR's own counts at every level.

Units matter more than anything else here. STAR counts read pairs. A BAM holds alignment
records (~2 per mapped pair, but not exactly — a pair whose mate failed to align contributes one).
Mixing them silently is how you get a 3-point error and call it biology; §4 of the original post is
one example already. So every row below says what it counts.

One cell, all the way down: day1_N2_1

Level 1 — mappability (read pairs, chimeric arm)

fate	read pairs
input	9,039,826
uniquely mapped	5,467,796
mapped to multiple loci	2,932,567
mapped to too many loci	17,128
unmapped: too short	611,164
unmapped: too many mismatches	0
unmapped: other	11,171
sum	9,039,826 — exact

Level 2 — which organism (uniquely-mapped only)

	ce11	ecHT115
alignment records	10,904,111	21,765
fragments	5,456,405	11,391
of which singletons (one mate only)	8,699	1,017

5,456,405 + 11,391 = 5,467,796 — exactly STAR's uniquely-mapped pair count. Fragments with one
mate on each organism: 0.

Level 3 — the multimapping bucket

	ce11	ecHT115
fragments	2,930,490	2,077

Sum = 2,932,567 — exactly STAR's multi-loci count. But read what this is: it says where the
one emitted record landed, not what the read is ambiguous between. Under --outSAMmultNmax 1
STAR writes only the top-scoring alignment while NH still counts every locus, so the other NH−1
placements exist nowhere on disk.

Level 4 — molecules

	ce11	ecHT115
kept records carrying UB	2,352,046	3,807
kept records without UB	8,552,065	17,958
distinct UB (of 65,536 possible)	61,287	1,974
distinct (UB, ref, pos) — proxy	1,671,160	2,918

The plain arm's real combined.h5ad for this same cell: 376,940 molecules, 11,708 genes
(verified directly against the file). So the positional proxy over-counts real per-gene molecules by
4.43x here, and by 2.56–5.09x (median 3.56x) across the 16 cells.

That gap is not error, it is the annotation step: the real counter collapses per gene and drops
no_feature and ambiguous, while the proxy counts every distinct (UMI, position) including introns
and intergenic space. Read the chimeric arm's molecule numbers as "roughly 2.5–5x the truth", never
as molecule counts. The chimeric arm has no real matrix to calibrate against, because the split
refused — that is the same bug, showing up as a measurement limit.

Level 5 — saturation, defined as 1 − distinct(UB, ref, pos) / kept records carrying UB:
ce11 0.289, ecHT115 0.234. The plain arm scores 0.290 on the same reads, as it must.

All 16 cells

	value
input	222,267,418 read pairs
uniquely-mapped fragments	134,047,626 — 131,867,306 ce11 + 2,180,320 ecHT115
multimapping fragments	67,360,622 — 66,617,645 ce11 + 742,977 ecHT115
Level 2 + Level 3 residual	0 on every one of the 16 cells
cross-organism fragments	0 on every cell
saturation, ce11	median 0.287 (0.238–0.634)
proxy / real UMI ratio	median 3.56x (2.56–5.09x)

ecHT115 fragments rise ~8x from day1 (238,958) to day17 (1,941,362) — the same direction and
rough magnitude as the 14x share effect in the original post, measured a different way.

One outlier worth a name. day17_N2_7 has saturation 0.63/0.69, roughly double every other cell,
and the fewest genes detected of any day-17 cell (10,924). It is also the cell with the lowest
multimapper share on the plate, 6.1% against a 40.5% maximum — see #427. A low-complexity library
sequenced deep, not a pipeline artifact.

A methodological note, because it nearly went in wrong

An intermediate version of this analysis asserted that mates in these BAMs do not share a query
name. That is false — checked directly:

V350394876L4C006R06201145106  163  V__ce11  1604  =  1671
V350394876L4C006R06201145106   83  V__ce11  1671  =  1604

Same name, both mates, standard SAM convention; 4M records carry 2,024,071 distinct names. The
fragment counts above were reconstructed positionally (FLAG plus RNAME/POS ⟷ RNEXT/PNEXT) rather
than by name, which is a valid but roundabout route — and it is validated by the thing that matters:
it reconciles to STAR's own totals with zero residual on all 16 cells, which a broken pairing
could not do. Recorded here because the earlier diagnosis in this thread does rely on query-name
singletons, and that reasoning stands.

What the ledger cannot say

1. Whether a multimapper's other loci were on the other organism. --outSAMmultNmax 1 never
   writes them. Not answerable by re-analysis at any effort; it needs a re-map with the flag raised.
   No proxy was built, deliberately.
2. True per-gene molecules for the chimeric arm — blocked by the split bug, hence the calibrated
   proxy above.
3. dropped_secondary from the CRAM is definitionally 0, since io cram strips 0x100 before
   the file exists. It happens to also be a true zero under --outSAMmultNmax 1, but the CRAM cannot
   be the evidence for that.

Full tables: pilot16/script/ledger/read_fate.tsv (16 rows x 60 columns, units in every column name)
and summary.json.

[lhqing]

Correction: "cross-organism fragments: 0" is tautological, and I published it as a finding

The ledger comment above reports "Fragments with one mate on each organism: 0" and "cross-organism
fragments: 0 on every cell". The counts are accurate. The implication — that this is a measurement
of anything — is wrong, and the data owner caught it.

Tested directly on day1_N2_1's chimeric BAM, all 16,784,546 records:

test	count
records whose mate is on a different reference — any pair of references, not just cross-organism	0
records with RNEXT == *	16,180
records with mate-unmapped flag 0x8	16,180

Not one cross-chromosome pair of any kind. And --chimSegmentMin appears in none of the four mapping
modules, so STAR's chimeric detection is off: its paired-end model only reports a pair whose mates
align concordantly to one reference, and anything else is never emitted as a pair. A worm/bacterium
molecule therefore cannot appear as a cross-component fragment. The zero restates a flag setting.

This weakens one leg of the diagnosis above

The main post lists, among split.py's stated assumptions, that mates "sit on the same chromosome
[component]" — HELD, 0 templates with kept mates on different components. That is true and it is
vacuous. It is guaranteed by the aligner's configuration rather than observed on real data, which
is precisely the distinction split.py's docstring draws when it says these facts "were read off the
aligner's source and nobody has yet watched them hold on a real chimera." One of the two has now been
watched (the NH premise, genuinely). The other has not, and cannot be by this experiment.

It also names a latent hazard: enable --chimSegmentMin for any reason and the premise stops holding
silently, with the parity check wrong in a second and different way. Whichever fix #426 lands on
should not treat same-component mates as an invariant it can rely on.

Where cross-organism evidence would actually live

Three populations, none of them the one that was counted:

1. The singletons — a genuine cross-organism molecule fails to pair, so STAR emits only the mate
   that aligned, 0x8 set. That is the same 16,180 records as the split-parity bug, and nobody has
   looked at what their missing mates were. Recovering the unaligned mate sequence (re-map with
   --outSAMunmapped Within KeepPairs, or go back to the uBAM) and asking where it lands is a real
   experiment that this pilot did not run.
2. The multimapper bucket — reads whose sequence is ambiguous between the organisms, still
   unmeasurable under --outSAMmultNmax 1.
3. unmapped: too short — pairs where neither mate cleared the 0.66 filter.

The bacterial share stated in this thread stays a lower bound for the reason already given
(cross-organism ambiguity is dropped as multimapping). This correction does not move that number; it
removes a line that appeared to bound the problem from the other side and does not.