Running the 784-worm chimeric plate: how it was driven, and what the driving taught us

[lhqing]

STATUS: LIVE — this is a running log, not a finished write-up. The six mapping shards are in
flight as of writing. Numbers marked provisional rest on partial data and will be replaced. This
post gets finalized when the run completes and the fan-in has produced both objects.

The scientific result belongs elsewhere: #425's deliverable is a dated measurement in
docs/research/, and the 16-worm pilot's findings are in #426. This is the logistics — how the
run was actually driven, and what the driving taught us. Engineering findings are being filed as
#448 and #449.

What is being run

The in-house SMART-seq3 aging plate — 784 worms, 9 ages x 3 strains, one animal per library — mapped
against the chimeric ce11_ecHT115 reference so the bacterial share of each library becomes a number.
Chimeric arm only; the plain arm was deliberately not re-run at 784 (see "Decisions taken" below).

Workspace aging_SS3/seqforge, pipeline ss3-784-chimera-c54b1418fa98, six shards across
cpu03-cpu06 + fat01 + fat02, WORKFLOW_VERSION 2026.8.18.

How it was actually driven

Enough detail to repeat it. Every seqforge/snakemake invocation goes through script/sfrun, the
container shim.

1. Refresh the code — on the LOGIN node, because compute nodes have no internet and git is not
even installed on them (CentOS 7 base).

2. Rebuild the metadata chain — make_metadata.py -> harvest verify -> mk_assertions.py.
Deterministic and offline: the drafts are parsed from directory names, so no model runs. Regenerated
rather than reused; the output was byte-identical to the August files, which confirms harvest verify
has not moved.

3. Fill the manifest — manifest fill <1568 files> --records --assertions --organism 6239 --offline.
Not seqforge run --no-llm: that path skips harvest and leaves assertions_path = None, so the
manifest would come out with no age and no strain and nothing would say so.

4. Fresh recipe, fresh --id — processing new --assembly ce11_ecHT115 --annotation wormbase_ws298+refseq_rs_2025_06_26 --threads 8 --mem-gb 48 --id ss3-784-chimera. Fresh because
run_id folds the recipe's fill-time workflow stamp, so re-composing an older recipe re-stages into
a directory that run already owns (#447).

5. Compose, then dry-run — 5493 jobs over 784 cells, module map/star-umi-chimera selected off
the assembly name with no second flag.

6. Materialize the genome index ONCE — snakemake --cores 1 results/index/ce11_ecHT115 — before
any shard starts. See hazard 2.

7. Partition — partition.py does longest-processing-time-first bin packing weighted by fragment
count against heterogeneous nodes, emitting one list per node. fat nodes get 234 cells (20 slots),
cpu nodes 79 (7 slots), so all six are projected to finish together.

8. Six shards, each sbatch --partition=... --nodelist=<node> --cpus-per-task=<n> --export=ALL,BATCH=<node> run-chimera-shard.sbatch, each asking for its own cells' four leaves and
never for all.

9. The fan-in, one job, gated --dependency=afterany on all six — afterany not afterok,
because --keep-going makes a non-zero exit the normal outcome of a shard that refused one cell.

Logistic observations

Things that would have silently produced a wrong or incomplete result

The obvious file glob is wrong. fastq/*/*.fq.gz yields 1597 files; the dataset has 1568.
The extras are 28 per-lane source files under day13_CF_26..32 — already concatenated into the
merged pairs, so passing both would have doubled every read for seven worms — plus one orphan mate.
It surfaced only because Singularity caps its argv JSON at 128 KiB and refused, initially with a
misleading error about the runtime engine; bisecting argument counts produced the real message. The
shell never complained. Use the file list records.yaml declares, not a glob.

Six shards race rule genome_index. It is shared upstream work and is not part of the disjoint
per-cell split. Its run: block does symlink_to(), which raises on an existing link, and snakemake
removes an output before re-running a job. Up to five of six shards can lose their entire run.
Materialize it once first. By contrast rule load_genome's touch() flag is safe to race, and the
difference is not visible from either rule.

--export=BATCH=<node> strips the environment. The bare form sets SLURM_EXPORT_ENV to that one
name, so the job starts with no PATH or HOME. Must be --export=ALL,BATCH=<node>.

--notemp is mandatory for a shard, and the reason is module-specific. split_chimera's
per-Component BAMs are temp() with two consumers: unique_to_cram per cell, and the plate-wide
umi_count. A shard's DAG sees only the first and reclaims them when the CRAMs land; the fan-in then
finds every input missing, re-runs split_chimera, whose own input is the STAR BAM — also temp(),
also reclaimed — and silently re-maps the whole plate.

--nolock is required and throws away the guard that catches a bad split. Six instances cannot
share one lock. The only thing between that and a corrupt output is a human verifying the lists are
disjoint. Do it explicitly: union must be 784 lines with zero duplicates.

<cell>.qc.json.gz is NOT a completion signal. rule qc_bundle does not depend on the CRAMs, so
it lands before them. A shard naming only the QC bundle — the natural shorthand, since it is what the
report reads a cell's row from — would finish every cell with its archives never written, and nothing
downstream consumes an archive, so nothing would say so. Name all four leaves.

Measurement traps

sacct is unavailable on this cluster (accounting DB refuses connections), so MaxRSS per job does
not exist. Any resource figure has to come from your own sampler.

The peak-RSS sampler in the run scripts measures NODE-TOTAL, not per-process. It sums RSS across the
process session, which double-counts the shared genome segment once per attached STAR job. Honest
per-rule memory needs /proc/<pid>/smaps_rollup (Pss, private RSS), because the index is one ~1.3 GB
SysV segment with nattch equal to the number of attached STAR processes — not one copy per process.

ps percentages are instantaneous. For setting a threads: declaration what you want is
(utime+stime)/elapsed integrated over a process's whole life.

Scheduling and utilisation

Nodes run well below saturation, and the cause is a thread declaration. Observed mid-run: cpu03 at
44% (load 24.6/56), fat01 at 72% (115.8/160). split_chimera holds an 8-core slot and delivers
~1.7 cores — 24 OS threads, most blocked behind a serial per-record Python routing loop — stranding ~6
cores per running job. STAR by contrast converts its 8 threads at ~94% efficiency. Filed as #449; the
structural point is that the module already solves this on the memory axis (memory.py turns one recipe
figure into three per-rule-class requests) and has no equivalent for threads.

Work is back-loaded. Each shard runs most of its mapping before any of its conversion, so a
completed-cell count stays near zero for hours and then finishes in a rush. Partly the thread
reservation above; possibly also --notemp removing snakemake's temp-reclamation scheduling incentive.
Do not read early low completion as trouble.

A canary shard was worth less than expected. One shard was submitted alone to validate before the
other five. Because conversions queue behind all of that shard's mapping, the end-to-end check it
existed to provide was hours away, while the target list was already evidenced by a finished
pilot16r2 cell under the same workflow file. Cost: the other five nodes idle. If the stage you want
to validate is at the back of the DAG, a canary does not validate it — find the evidence elsewhere.

Things that were fine, contrary to expectation

Log volume is a non-issue: 0.24-0.61 MB per shard after four hours. Worth recording because the
August run left a 7.2 GB .err for a single mapping job. The difference is #429 declaring STAR's run
files as rule outputs instead of letting them stream to stderr — a fixed regression, noted so nobody
re-opens it.

Snakemake's coordinator overhead is modest at 5493 jobs: 199 MB RSS, .snakemake metadata 6.6 MB
over 5404 files. (It does carry ~384 threads at ~47% CPU sustained, which is unexplained but harmless
here.)

Memory is nowhere near the request. The recipe asks 48 GB per mapping job and derives a 36 GiB
--limitBAMsortRAM cap; measured peak on a whole 20-slot node is a small fraction of that. The cap
permits rather than reserves. A mem_gb ticket is still owed.

Decisions taken, and what they cost

Chimeric arm only, no plain 784 arm. #425's stated methodology was both arms under the same code,
because two of its five items are arm-vs-arm comparisons. Those now rest on pilot16r2's 16-cell
comparison — defensible, since both arms there ran under #429-era code on the same cells, but the
research artifact must say so and must not present a 784-cell bacterial share beside a 16-cell
multimapper bound as one result.

--mem-gb 48 kept unchanged despite being ~20x the measured peak. In the shard model
(--cores N, no --resources mem_mb) it throttles nothing, so cutting it buys no throughput while
risking a FATAL on the deep tail; and keeping it identical to the pilot's recipe is what keeps the
16-cell result comparable.

Two fat nodes and four cpu nodes; cpu07/cpu08 left free for other lab members.

Corrections to earlier claims in this run

Recorded because they were acted on before being checked:

• The bacterial share is not a lower bound. Run the in-house SMART-seq3 plate against the built chimera, and get the bacterial fraction #425 says an organism-ambiguous read is "dropped as a
  multimapper". It is not — split.py routes it by its representative record's Component and keeps
  it, marked. multiplaced is a subset of kept; dropped holds only unmapped/secondary/supplementary.
  The point value is contaminated in both directions, and the honest bracket is wide.
• Three different statistics have all been reported as "the bacterial share": pooled depth-weighted
  share of kept records (1.445%), per-cell median of the same (0.978%), and per-cell median of
  unique/kept_total (the 0.144% -> 2.081% day1/day17 figures). They differ by up to 1.5x. Label which
  one you mean.
• Timing estimates were revised twice — first from a stale partition.py constant (~6x
  conservative), then upward again once node utilisation turned out to be 44-72% rather than saturated.

Open items, to be resolved before this is finalized

• All six shards complete; fan-in produces combined.ce11.h5ad and combined.ecHT115.h5ad
• The fan-in's own cost at 784 x 2 components — the one rule nobody has ever measured at any scale
• Real per-cell disk and the --notemp multiplier (provisional: ~3.1 GB/cell, ~3 TB peak)
• Per-rule threads and memory profile -> what each rule actually earns (in flight)
• split_chimera cost normalised to records/sec, and the profiling research task on split_chimera spends 81% of its wall time rebuilding records in Python, over half of that round-tripping aux tags #449
• Failure count and modes across 784 cells
• seqforge report behaviour at 784 samples (that code was exercised on 16)
• Total wall clock, and the gap against each successive estimate

[lhqing]

FINAL — the run completed 2026-08-21. Every open item below is resolved.

784/784 cells, run_state: finished, no missing deliverables, 0 alerts.

The eight open items

• All six shards complete; fan-in produced both objects. combined.ce11.h5ad (784 x 46,926
  genes, 315 MiB) and combined.ecHT115.h5ad (784 x 4,555, 12.2 MiB).
• The fan-in's own cost at 784 x 2 — the rule nobody had ever measured. Wall 28 min 46 s;
  peak 16.43 GB PSS per Component tree at 64 threads, 2.06 GB at 8. It scales with thread
  count, which fan_in_mem_mb (max(8 GB, mem_gb/4) = 12 GB) does not model — a ~6x
  over-request at 8 threads and a 34% under-request at 64. Filed on One threads figure serves every rule, so a plate-wide fan-in takes 16 of 160 cores while a split reserves 8 to use 2 #452.
• Real per-cell disk and the --notemp multiplier. Peak 2.90 TB; 255 GB after cleanup.
  8.7x, and see the correction below — most of it was not reclaimed automatically.
• Per-rule threads and memory. Integrated CPU against a uniform threads: 8: STAR 7.1,
  unique_to_cram 4.0 (ce11) / 1.2 (ecHT115), multiplaced_to_cram 3.2, split_chimera
  1.7, qc_bundle 0.8 of 1, fan-in wanting ~64. Filed as One threads figure serves every rule, so a plate-wide fan-in takes 16 of 160 cores while a split reserves 8 to use 2 #452.
• split_chimera normalised. 19.45 s per million records, linear (r = 0.920), no fixed
  overhead — the flat per-cell figure was the wrong unit. 0.548x STAR median. Bottleneck is the
  per-record Python rebuild at 81.4% of wall, over half of that the aux-tag round-trip; BGZF is
  12%, and lowering it is a pessimisation. A copy() path measures 3.01x with
  byte-identical output. split_chimera spends 81% of its wall time rebuilding records in Python, over half of that round-tripping aux tags #449.
• Failure count and modes. 3 of 784 (0.38%), all one cause, all split_chimera — a
  multiply-placed survivor whose unmapped partner STAR never wrote. A multiply-placed survivor whose unmapped partner STAR never wrote makes split_chimera refuse healthy cells — the mirror of ADR-0053 #451, fixed and deployed
  mid-run; the 3 cells re-ran clean and are statistically unremarkable (largest robust |z| = 1.79).
• seqforge report at 784 samples. Renders correctly, 4.1 MB self-contained HTML, and says
  finished with missing_deliverables: []. Worth noting: in the pilot this same report called a
  broken run 16/16, alerts: [] while both deliverables were missing. Finalize SS3 chimera support: unblock the split, partition the archives by mappability, bundle the QC, and say when a run failed #429 fixed that, and
  this is its first exercise on a complete plate.
• Total wall clock, and the gap against each estimate. Mapping ~6.4 h per shard; fan-in 29
  min. See below.

Estimates against outcome

estimate	when	actual
~2.7 h mapping	from partition.py's stale constant	—
~5 h	after re-deriving from pilot QC bundles	—
6-7 h	after measuring STAR throughput	~6.4 h per shard
8-10 h	after finding nodes at 44-72% utilisation	(over-corrected)
1-3 h fan-in	extrapolated from the 8-thread pilot	29 min at 64 threads

The fan-in estimate is the instructive one. It would have been correct at the declared 8 threads —
and the job would have been 8x slower than the node allowed. A prediction that matches an outcome
does not mean the outcome was right.

Corrections to this log

--notemp does not defer reclamation. It cancels it. This log said the fan-in would reclaim what
the shards retained. It reclaimed only the per-Component split BAMs (700.7 GiB, exactly and
completely) because its executed DAG was 3 jobs; every other temp() belonged to jobs already
complete and never in the DAG. ~2 TB was permanently orphaned and needed a manual sweep. Recorded
on #448 — a partitioned-execution feature cannot treat --notemp as a flag to pass.

Log.out is declared, contrary to a monitoring report during the run. It is the first member of
STAR_PROGRESS_LOGS in workflows/h5ad.py, imported and declared temp() via STAR_PROGRESS. A
grep of the .smk for the literal string misses it because the constant lives in the shared module —
which is the indirection that file exists to create. There is no undeclared-file bug.

Peak memory figures quoted mid-run were node-total. The samplers in these scripts sum RSS across
the process session, double-counting the shared genome segment once per attached STAR job. Every
honest per-rule number here comes from /proc/<pid>/smaps instead. The trap is documented above and I
still fell into it once.

What this run produced beyond the data

Four issues, one PR, one measurement, all from driving it:

• A large plate has to be sharded across sbatch jobs by hand, and every hazard is rediscovered each time #448 — partitioned execution as a first-class feature, and the container/Slurm root cause
• split_chimera spends 81% of its wall time rebuilding records in Python, over half of that round-tripping aux tags #449 — split_chimera's cost is a Python record rebuild, not compression
• A multiply-placed survivor whose unmapped partner STAR never wrote makes split_chimera refuse healthy cells — the mirror of ADR-0053 #451 — the parity bound is two-sided (fixed, PR The multiply-placed survivors bound the shortfall, and the split stops refusing it #453)
• One threads figure serves every rule, so a plate-wide fan-in takes 16 of 160 cores while a split reserves 8 to use 2 #452 — one threads figure serves every rule
• PR The plate says what fraction of a worm library is its food #454 — docs/research/worm-bacterial-fraction-aging-plate.md

The result

5.886% pooled bacterial share, 0.753% by the typical worm, working range 4.93-5.89%.
rho(share, age) = +0.699 across 784 cells, holding inside every strain. Day1 -> day17 is 42x on
the per-cell median.

The rise is real; the monotonic reading is not — day3 falls below day1, day7 returns to baseline,
and the peak is day15. Full numbers and caveats in PR #454.

For the next session

• PR The multiply-placed survivors bound the shortfall, and the split stops refusing it #453 is open and unmerged, and /share/home/lhq/src/seqforge is checked out on that
  branch, not main. Switch it back once merged.
• 3 of 784 cells carry two payload keys the other 781 lack (unanswered_survivors,
  multiplaced_singletons) — they were split under the A multiply-placed survivor whose unmapped partner STAR never wrote makes split_chimera refuse healthy cells — the mirror of ADR-0053 #451 fix. Treat them as optional in any
  plate-wide analysis.
• Artifacts: aging_SS3/script/metrics784/{final,split_timing,resource_profile,split_parity,fanin}/.

[lhqing]

Another engineering finding off this run, filed as #457.

Opening day9_N2_9.ce11.unique.cram in IGV turned up reads carrying "introns" of up to 840,537 bp
— in a genome whose longest annotated intron is 86,515 bp. They are in the .unique archive at
MAPQ 255/NH:i:1, so no uniqueness filter removes them.

Cause: none of the three map modules set any splice parameter, so STAR's derived alignIntronMax
(~590 kb) and --alignSJoverhangMin 5 apply. That default is mammal-calibrated — measured against
the GTFs installed on ircbc, the longest intron is 1.24 Mb in hg38 and 870 kb in mm39, but 86.5 kb
in ce11 and 2,483 bp in sacCer3. So the flag we never set is about right for human/mouse and
wrong by ~7x for worm, ~238x for yeast.

The ecHT115 arm turns out to be a free control: E. coli has no introns, so every N gap there is
spurious by construction — and 100% of them carry the same signature (mean 9 bp anchor, ~128 bp
soft-clip) as the worm-arm reads.

Numbers, method and what is not claimed are in #457. Not fixing it here — the right ceiling is
genome-dependent, so where that value lives is a design question rather than a one-line edit.

[lhqing]

The multimapper question this run opened is answered — #427 is closed

This log recorded that 32.593% of kept records are multiply placed, and used that to call #425's
37.5% upper bound near-vacuous on the reasoning that they are "overwhelmingly within-worm repeats".
That was right, and they now have names.

They are ribosomal RNA. The chrI 45S rDNA array is collapsed in ce11 to ~2.1 copies of a 7,197 bp
unit, so an rRNA read can never map uniquely; rrn-1.1/rrn-1.2/rrn-3.1/rrn-3.56 are 53.09%
of all placed molecules on the plate.

And the NH=3 that is 81% of the bucket is not three copies of anything — it is two copies counted
three ways. A single mate has two placements, but STAR places a pair and this recipe bounds no
mate gap, so both-in-A, both-in-B and one-in-each all score the same. A mate-gap sweep steps from
14.41% to 84.28% exactly at 7,197 bp, and the same reads run single-end collapse from 83.45% to
3.73%.

Two things from this log get sharper as a result:

• The same "which statistic do you mean" trap this log flagged for the bacterial share applies
  here too. STAR's read-level share and the counter's fragment-level share differ by a median of 3
  points and a maximum of 50, and rho between them is only 0.79 — because the counter's denominator
  is ce11 fragments while STAR's is all input reads. day17_DA_1 reads 7.90% against 47.88% for
  exactly that reason.
• The saturation lead from What ARE the multi-mapped reads? 6.1%-40.5% across one plate, and nothing says why #427's second comment is dead, and backwards. Across 784 cells the
  correlation is +0.354 — more saturated cells have more multimappers, not fewer.

The full report, with every table and the ruled-out list, is
discussion #474. The dated measurement is
docs/research/multimapper-composition-aging-plate.md, merged in #473, and the artifacts sit beside
this run's others at aging_SS3/script/metrics427/.

Three follow-ups filed: liuhlab/liulab-genome#111, #471, #472. The recommendation was no bar and no
panel — the share is a proxy for rRNA carry-over, so map/star-umi and map/star-umi-chimera stay
in MODULES_WITHOUT_CROSS_CHECKS, now on evidence rather than on absence of it.