diff --git a/docs/research/chimera-offline-detection.md b/docs/research/chimera-offline-detection.md new file mode 100644 index 0000000..11d6cbb --- /dev/null +++ b/docs/research/chimera-offline-detection.md @@ -0,0 +1,331 @@ +# Can compose tell a chimera from its name alone, offline? + +Read 2026-08-14 out of the `liulab-genome` source for +[#407](https://github.com/liuhlab/seqforge/issues/407), under map +[#406](https://github.com/liuhlab/seqforge/issues/406). **Yes for the two questions the map's +constraint 5 actually asks — a chimeric name is decidable, and its component list recoverable, from +the name plus the wheel's own shipped table, with no read of the machine's reference store and no +network. But the offline answer is a different proposition from the on-machine one, and the two +disagree in both directions: `liulab-genome` itself states that "the record, never the metadata row" +decides whether a built assembly is a chimera, and it ships a test proving that `ce11_ecHT115` +registered from somebody's own FASTA is not one. Offline detection is therefore a claim about the +NAME, and it is sound only if every component is a row in the shipped +`assembly_metadata.tsv` — which today is 7 assemblies, so exactly 120 chimeras of any arity are +detectable and every other one needs a `liulab-genome` edit and release first. The map has not priced +that. Per-component `ncbi_taxid` is reachable offline; per-component annotation NAME is not — +the merged name is not injective, by that function's own docstring.** + +This is a measurement, not a decision. It reports what the code on disk does today. What to do about +it — whether compose selects the module, whether the `.smk` decides at parse time, whether a +processing field declares it — belongs to the compose-selection ticket and to whatever record it +writes. + +Two source trees are involved and they are not at the same version, which is itself a finding (§2). +Line citations are against the sibling checkout `/Users/hanqingliu/src/liulab-genome` at +`bd0d4d5` on branch `refactor/94-surfaces`; seqforge citations are against this worktree. + +--- + +## 1. The four checks, in order, for `ce11_ecHT115` + +The four checks ADR-0008 names are implemented in one function, +`genome/io/source.py::resolve_source`, and are enumerated in its own docstring at +`src/genome/io/source.py:198-219`. What each one does for the name `ce11_ecHT115`: + +| # | check | code | fires? | reads disk? | +|---|---|---|---|---| +| 1 | a completion record here | `source.py:256-261` | **yes**, always | **yes** — `AssemblyDir.read_record()` | +| 2 | a source the caller named | never reaches this function | no | n/a | +| 3 | the name | `source.py:262-277` | **yes**, and it decides | **partly** — see below | +| 4 | today's fetch path | `source.py:260`, `265`, `267` | no, not reached | no | + +**Check 1 is unconditional and it is a disk read.** `record = assembly_dir.read_record()` +(`source.py:256`). If a record is there it is believed outright and the name is never consulted: +`ChimeraDetails.from_record(record)` returns the components for a chimera's record, or `None` for +any other record, in which case the assembly is resolved as an ordinary fetch (`source.py:257-261`). + +**Check 2 never arrives.** `--source` / `Genome(path_or_url=...)` builds a `SeededSource` handed +straight to the registration; the docstring says so at `source.py:205-210`. There is no +compose-time equivalent, so it is not a check compose would have to reproduce. + +**Check 3 is the whole of the offline question, and it is two halves.** The syntactic half is pure: + +```python +# src/genome/io/source.py:262-267 + try: + candidates = split_name(assembly_dir.assembly) + except ChimeraNamingError: + return fetched_source(metadata, golden_path_url) + if not all(_could_be_a_component(name) for name in candidates): + return fetched_source(metadata, golden_path_url) +``` + +`split_name` (`src/genome/chimera.py:154-197`) is a `str.split` and a regex over the parts and +nothing else — `chimera.py:189-190`: + +```python + parts = tuple(name.split(_NAME_JOIN)) + if len(parts) < _MIN_COMPONENTS or not all(_COMPONENT_RE.fullmatch(part) for part in parts): +``` + +The module's own docstring states the guarantee: *"Pure means names in, names out: nothing here +opens a file, and nothing here imports from `genome.io` or `genome.genome`"* (`chimera.py:10-11`). +Confirmed by the import list at `chimera.py:51-55`: `re`, `collections`, `collections.abc`. + +The semantic half is `_could_be_a_component`, and this is the one line that touches the machine — +`src/genome/io/source.py:183`: + +```python + return is_prepared(assembly) or lookup_assembly(assembly) is not None +``` + +`is_prepared` (`source.py:150-172`) is `read_record(assembly_data_dir(assembly)) is not None`; the +directory comes from `LIULAB_DATA` or a well-known lab root +(`src/genome/io/registration.py:79-135`). It never raises when the root is missing — it falls back +to `~/liulab_data` and simply stats a path that is not there — so a caller can be offline and still +get an answer. + +**The disk read in check 3 is a disjunct, and for `ce11_ecHT115` it is the losing one.** `ce11` and +`ecHT115` are both rows of the shipped `src/genome/data/assembly_metadata.tsv`, so +`lookup_assembly` answers each of them from a file inside the installed wheel — `metadata.py:397-401` +reads it through `importlib.resources.files("genome")`, cached (`metadata.py:355-373`). Same bytes +on every machine, same answer with the reference store deleted. + +Then the canonical-order check, pure again (`source.py:268-277`): `derive_name(candidates)` is +compared to the given name, and a mis-ordered spelling raises `FileNotFoundError` naming the +canonical one. + +**So: for `ce11_ecHT115` on a machine holding nothing, check 1 fires and reads the disk, check 3 +fires and decides, check 4 is not reached — and only check 3 is reproducible offline.** Dropping +check 1 is exactly what makes an offline answer a different proposition, which is §5. + +Run against the sibling checkout's source with `LIULAB_DATA` unset and no chimera anywhere on this +laptop: + +```text +name split_name all parts listed? canonical +ce11_ecHT115 ('ce11', 'ecHT115') True ce11_ecHT115 +ecHT115_ce11 ('ecHT115', 'ce11') True ce11_ecHT115 <-- MIS-ORDERED +hg38 ChimeraNamingError - - +my_ref ('my', 'ref') False my_ref +hg38_mm10 ('hg38', 'mm10') True hg38_mm10 +tinyCe_tinyEc ('tinyCe', 'tinyEc') False tinyCe_tinyEc +sacCer3_hg38 ('sacCer3', 'hg38') True hg38_sacCer3 <-- MIS-ORDERED +test-star ChimeraNamingError - - +hg38_mm10_sacCer3 ('hg38', 'mm10', 'sacCer3') True hg38_mm10_sacCer3 +``` + +--- + +## 2. Is a table row required for every chimera? Not for the chimera — for every component + +**The chimera's own row plays no part in detection.** Check 3 splits the name and looks up the +*parts*; the whole name is never looked up. The `metadata` argument `resolve_source` receives is +used only by `fetched_source` (`source.py:146-147`), which is check 4, and its only call site +(`src/genome/io/download.py:228-232`) passes `lookup_assembly(assembly)` for that purpose alone. Delete +the `ce11_ecHT115` row and `resolve_source` still returns `ComponentSource(("ce11", "ecHT115"))`. + +`liulab-genome`'s own test says the same thing in the same words: `tests/test_metadata.py:45-55` +defines what makes a shipped row a chimera's — *"It splits into two or more parts and the table lists +every one of them"* — and `tests/test_metadata.py:237-251` asserts the `ce11_ecHT115` row is name +and nothing else, existing *"so that a machine holding neither component can still tell this name +from a free-form local key, by splitting it into components the table lists"*. Confirmed by reading +it back: `AssemblyMetadata(assembly_name='ce11_ecHT115', species=None, ucsc_name=None, +ncbi_name=None, ncbi_assembly_id=None, ncbi_taxid=None, source_url=None, sha256=None)`. + +**Every COMPONENT, though, must be a row — and this is the constraint the map has not priced.** +Offline, `is_prepared` is unavailable to compose by construction (R7: no machine fact may decide a +compose-time output; the genome deliberately resolves at run time in `rule genome_index`, +`src/seqforge/workflows/map/star-umi.smk:251-253`). That leaves `lookup_assembly` alone, so the set +of chimeras compose can detect is exactly the set whose every component is one of the seven rows the +wheel ships: + +```text +hg38 hg19 mm39 mm10 sacCer3 ce11 ecHT115 +``` + +That is 21 pairs, 35 triples, … — 120 chimeras of arity ≥ 2 in total, and **no others, ever, until +`liulab-genome` gains a row and cuts a release.** Say that plainly: + +- A chimera involving any organism not on that list — a PDX with a new mouse strain, a co-culture + with a second bacterium, a spike-in genome, *any* locally seeded reference — is **undetectable by + compose** until someone edits `src/genome/data/assembly_metadata.tsv`, releases `liulab-genome`, + and re-pins seqforge's dependency. On-machine resolution has no such constraint: `is_prepared` + covers a locally registered component the moment it is built. +- The fixtures the map leans on are already outside the set. `tinyCe`, `tinyEc`, `tinySc` and + `tinyEcDub` are not table rows, so `tinyCe_tinyEc` is **not** offline-detectable — see the run in + §1. `liulab-genome`'s own `tests/test_source.py:139-157` proves this is the intended behaviour and + that only `is_prepared` rescues those names. **The map's cheap bar (constraint 6) is a synthetic + round-trip on exactly those fixtures.** If compose's selection rule is "the table lists every + part", a fixture round-trip cannot exercise it without a monkeypatched table — which is testing + something other than what ships. +- **The pin is already behind.** This worktree's `pixi.lock` pins `liulab-genome` at + `ab3272312eb3f9ef25e3ac8320aac28187af42f5` (2026-07-23, "add a Chromap aligner"). The installed + package at `.pixi/envs/default/.../site-packages/genome/` has **no `chimera.py`, no + `io/source.py`, no `io/components.py`, no `ecHT115` row and no `ce11_ecHT115` row** — its table + stops at `ce11`. So today seqforge cannot detect a chimera offline at all, and the very first + chimera ticket owes a dependency re-pin before a line of detection code can run. This is the + "edit and release before composing" cost, already binding, not hypothetical. + +--- + +## 3. The component list comes from the name, and only from the name + +`split_name` returns the candidates (`chimera.py:154-197`); `derive_name` re-derives the canonical +spelling from them (`chimera.py:106-151`). Both pure. `ComponentSource(candidates)` — the value +check 3 returns (`source.py:277`) — carries nothing the name did not. + +The record is **not** needed for the list. Where the record is load-bearing is everything *else* a +chimera knows about itself: `ChimeraDetails` (`src/genome/io/components.py:186-217`) carries the +**separator** its chromosome names were written with, and each component's `sha256` and contributed +annotation. `Genome.components` reads it, never the name — `src/genome/genome.py:397-418`: + +> *"The single test of whether an assembly is a chimera, and the completion record is what answers +> it — never the metadata table, which lists a chimera as a cross-reference and would answer the +> same question differently on a machine where the row is stale or absent (ADR-0008)."* + +That sentence is the boundary of this whole ticket. It is a direct statement, in `liulab-genome`'s +public API, that name-and-table is *not* the is-a-chimera test. An offline compose-time answer is a +second, weaker predicate, and it needs its own name so nothing confuses the two. + +Note what this costs the splitter: **the separator is not offline-recoverable.** `__` is only the +default `split_suffixed` assumes when no component carries a doubled underscore +(`chimera.py:295-362`), and `tinyEcDub` is the shipped fixture that breaks it. The map's constraint 1 +already puts the split on the `liulab-genome` side and reads the separator off the record at run +time, so this is consistent — but it does mean a compose-time artifact may not bake a separator in. + +--- + +## 4. Taxid: yes, offline. Annotation name: no + +**`ncbi_taxid` per component is reachable offline**, straight out of the shipped table: + +| component | `ncbi_taxid` | `species` | default annotation | +|---|---|---|---| +| `ce11` | `6239` | `Caenorhabditis elegans` | `wormbase_ws298` | +| `ecHT115` | `634469` | `Escherichia coli HT115` | `refseq_rs_2025_06_26` | + +The chimera's own row carries `ncbi_taxid=None` (§2), which is correct — a chimera has no single +taxid — and it is worth flagging against seqforge's current shape: `GenomeRef.ncbi_taxid` +(`src/seqforge/models/processing.py:42-44`) is a **single** optional taxid, and +`src/seqforge/manifest/policy.py:418-422` fills it from `dataset.experiment.organism.value` — the +*dataset's* asserted organism, never the assembly's. For a chimeric run that field describes one +component at best and is silently wrong at worst. Nothing today reads it into a command line, so +this is a shape note for the identity ticket, not a live bug. + +**The per-component annotation NAME is not offline-recoverable in general, and this is stated by the +code rather than inferred.** `merged_annotation_name` (`src/genome/io/components.py:84-116`) joins the +contributing annotations with `+` in sorted-component order, and its docstring says at +`components.py:94-99`: + +> *"It needs no parse-back: what a merged annotation is made of is recovered from the components, and +> written down in its own record besides. And it is not asked to carry *which* components +> contributed — a chimera with a component that contributes nothing spells the same name a different +> subset would."* + +So `wormbase_ws298+refseq_rs_2025_06_26` splits on `+` into two names, but nothing in the string says +which component each belongs to, and a 3-component chimera where one component contributed no +annotation yields a 2-element list against a 3-element component list. Positional alignment is a +guess. Confirming the second half: `annotation_metadata.tsv` has **no row naming `ce11_ecHT115`** — +ADR-0008's *"its merged annotation gets no row at all"*, verified. + +What *is* offline is each component's **default** annotation (the table above, `default: yes` for +both). For today's `ce11_ecHT115` that happens to be the truth — the built chimera's merged +annotation is `wormbase_ws298+refseq_rs_2025_06_26`, exactly the two defaults in sorted-component +order. It is a coincidence of the build, not a guarantee: a chimera built against a non-default +component annotation spells a merged name the table cannot reproduce. **Per-component annotation is +a run-time read of the record, full stop.** Map constraint 3 ("each counted against its own +annotation, which the chimera's record names per component") already says record; this confirms +compose cannot pre-compute it. + +--- + +## 5. False positives and false negatives — both exist, both have a shipped test + +A local key containing `_` splits happily; whether it then reads as a chimera turns entirely on +whether the table lists every part. `my_ref` → `('my', 'ref')` → neither listed → ordinary +assembly, which is the documented separation (`source.py:211-214`). `hg38` and `test-star` never +reach the split at all — one has no `_`, the other's parts fail `[A-Za-z0-9]+`. + +**False positive — a real, non-chimeric assembly whose name reads as a chimera's.** Concrete and +already tested in `liulab-genome`: `tests/test_source.py:173-186`, +`test_a_plain_record_keeps_a_chimera_shaped_name_whatever_it_was_registered_as`, whose comment is + +> *"The reason the record comes first: `ce11_ecHT115` seeded years ago from somebody's own FASTA is +> not a chimera, and no amount of the name looking like one may change what a finished registration +> already is."* + +The assembly whose name the map is built on is the test's own example. `hg38_mm10` is the same case +with two very ordinary parts — a name a lab would plausibly give a hand-concatenated reference built +years before this package existed, or a liftover scratch build. Offline, compose cannot tell it from +a real chimera, because the only thing that can is check 1, the record. Cost of getting it wrong: +compose selects the chimeric `.smk`, the run maps fine (the index is whatever is there), and the +split rule finds no `__` suffixes — a failure at split time on a reference that was never chimeric. + +**False negative — a real chimera that reads as an ordinary local key.** Also tested: +`tests/test_source.py:139-157`, `test_only_a_prepared_component_counts_when_the_table_lists_neither` +— *"Neither half is in the shipped table, so nothing but a record of its own can make +`tinyCe_tinySc` read as two assemblies rather than as one name somebody chose."* Every locally +registered component is in this class, and so is every one of the map's fixtures. Cost: compose +selects the plain module, the chimeric BAM is produced and nothing downstream knows — precisely the +gap the map says it exists to close, silently reinstated for exactly the references the lab builds +by hand. + +**A third case, neither of the above: the mis-ordered name.** `ecHT115_ce11` and `sacCer3_hg38` +split to listed parts but are not the canonical spelling, and `resolve_source` *raises* +`FileNotFoundError` naming the right one (`source.py:268-276`). A compose-time rule that mirrors +check 3 inherits that raise, which turns a legal free-form local key into a hard compose refusal for +a user who never mentioned a chimera. Offline compose has no record to fall back on, so it cannot +resolve the ambiguity the way `resolve_source` does. + +**Neither direction is closable offline**, and that is structural rather than a gap in the +implementation: check 1 exists precisely because the name is not authoritative, and check 1 is the +machine. + +--- + +## 6. Verdict + +**LEGAL-WITH-CONDITIONS.** Compose *can* decide, offline and with no network, that a name is spelled +like a chimera's and that every part is a shipped table row, and recover the ordered component list — +all of it from `chimera.split_name` plus `metadata.lookup_assembly`, both pure of the machine's +reference store. Constraint 5 is not illegal as written. The conditions are: + +1. **What compose decides is a property of the NAME, not of the reference.** It needs its own term — + *chimera-spelled*, say — kept distinct from `Genome.components`, which is the real test and is a + run-time record read. A compose-time refusal must say which of the two it is talking about. +2. **Detection is capped by the shipped table.** Seven components today, 120 chimeras; anything else + needs a `liulab-genome` row, a release and a seqforge re-pin *before* the dataset can be composed. + The map should price this, and should decide what compose does when a user names a chimera it + cannot see — silently plain, or a refusal naming the missing row. +3. **The pin is behind the feature.** The `liulab-genome` this repo installs has no `chimera.py` + at all. Re-pinning is the first task of the compose-selection ticket, not a footnote. +4. **Both error directions survive**, and each hits a real reference class — a pre-existing + `hg38_mm10`-shaped assembly, and every locally registered component including the map's own test + fixtures. Whatever compose emits must be overridable by a declared processing field, so a user can + say "this is not a chimera" or "this is one" and be believed. That override is also the only way + the constraint-6 fixture round-trip exercises the shipped selection rule rather than a + monkeypatched table. +5. **Nothing about the chimera's INTERNALS may be pre-computed at compose time** — not the separator + (`tinyEcDub` forces `___`), not the per-component annotation (the merged name is not injective). + Compose may select the module; the module reads the record. That is consistent with map + constraints 1 and 3, and it is the line that keeps R7 intact. + +--- + +## Method and caveats + +- Everything above is read from source. Nothing was built, no chimera was constructed, no reference + store was consulted beyond `~/liulab_data` existing on this laptop. +- Citations are against `/Users/hanqingliu/src/liulab-genome` at commit `bd0d4d5`, branch + `refactor/94-surfaces` — **not** a tagged release and **not** what seqforge installs. Line numbers + will move. +- The run in §1 and the tables in §2 and §4 were produced by importing that checkout's `src/` ahead + of the installed package (`PYTHONPATH`, `PYTHONDONTWRITEBYTECODE=1`), inside seqforge's `default` + pixi env, with `LIULAB_DATA` unset. `liulab-genome` was not modified. +- "120 chimeras" is `2^7 - 7 - 1` over the seven shipped rows — every subset of size ≥ 2. It counts + what the table *permits*, not what would build; only a prepared set builds one (ADR-0008). +- The claim that check 3 never consults the chimera's own row is by reading, not by deleting the row + and re-running. The reading is unambiguous: `resolve_source` looks up only `candidates`, and its + `metadata` argument reaches `fetched_source` and nothing else. +- Written up on branch `research/chimera-offline-detection` (not opened as a PR). diff --git a/docs/research/star-on-a-chimera.md b/docs/research/star-on-a-chimera.md new file mode 100644 index 0000000..22814eb --- /dev/null +++ b/docs/research/star-on-a-chimera.md @@ -0,0 +1,587 @@ +# STAR on a chimeric reference: what the BAM actually says about which component a read came from + +Measured 2026-08-14 for [#408](https://github.com/liuhlab/seqforge/issues/408), a research ticket +under map [#406](https://github.com/liuhlab/seqforge/issues/406). This is a **measurement of +behaviour**, not a decision: it says what STAR 2.7.11b does when the reference is one FASTA +concatenated from two or more assemblies with chromosome names suffixed `__`. What +that behaviour *decides* about the split contract belongs to the map and to whatever record the map +produces. + +Every claim below is labelled **VERIFIED** (read off STAR's own source, its `parametersDefault`, or +its manual, with a locator) or **INFERRED** (a consequence I reasoned to and did not observe). An +inferred claim marked as such is useful; an inferred claim dressed as verified is a defect, so the +two are kept apart on purpose. Nothing here was run: no STAR invocation, no cluster job. + +## The answer + +**The load-bearing claim survives in substance and fails in two named places, and one of those two +is fatal to the split as the map currently describes it.** In detail: + +- **A window is confined to one chromosome, so no non-chimeric alignment — and therefore no proper + pair — can span two components.** A genuinely half-host/half-contaminant fragment does not become + a cross-component pair; at default filters it becomes `unmapped: too short` and, because + `--outSAMunmapped` is never passed by any seqforge module, it does not appear in the BAM at all. + STAR's chimeric detection is OFF by default (`--chimSegmentMin 0`), which is the only machinery + that could ever have written a cross-component record pair. +- **`NH` and `MAPQ` are computed over the whole chimera, not over a component.** They are both + functions of `nTrOutSAM`, the count of loci reported anywhere in the concatenated genome. So the + map's phrasing is right for the reads it describes and wrong at the boundary: STAR reports every + locus within `--outFilterMultimapScoreRange` (default **1**) of the best, not only the best. A + read that beats its cross-species hit **by one point** is still `NH:i:2`, `MAPQ:3`, and a + component-local splitter would hand that record to the host BAM carrying a MAPQ a single-assembly + run would have written as 255. +- **`--outFilterMultimapNmax` is counted across the whole chimera.** 9 host loci + 2 contaminant + loci = 11 > 10 and the read is dropped as "mapped to too many loci". Confirmed; the flag that + changes it is `--outFilterMultimapNmax` itself, and no seqforge module passes it today. +- **The fatal one: `map/star-umi` passes `--outSAMmultNmax 1`, so its BAM contains exactly one + alignment record per template.** A splitter reading that BAM cannot observe "this template's + alignments span more than one component", because N−1 of them were never written. The map's + three-way routing rule is not implementable against that artifact as it stands. The ticket brief + states the opposite — that `star-umi` deliberately does *not* set the flag — and the brief is + wrong: `src/seqforge/workflows/map/star-umi.smk:453` passes it, and it is a settled module literal + (ADR-0022, via #256 decision 7). +- **"Unmapped" is not component-attributable, by construction**: an unmapped record carries + `refID = -1`, `POS = -1`, `MAPQ = 0`, `NH:i:0`, `HI:i:0`. The contract should say this outright. + The one wrinkle: the *unmapped mate of a mapped read* carries `RNEXT`/`PNEXT` pointing at a + chimeric chromosome, so an unmapped record can still name a component in a field a naive rewrite + would miss. +- **Beyond `@SQ`, the reference is named in `@PG CL:` and in `@CO`, both of which embed + `--genomeDir`** — the chimera's index path — verbatim. `@HD` and `@RG` do not name it. Any + "restore the header a single-assembly run would have produced" clause has to reach those two lines + or record an explicit exception. + +## What was read + +STAR **2.7.11b** — the version the `align-rna` image pins, per `src/seqforge/workflows/map/star-umi.smk:212` +("2.7.11b in the `align-rna` image (2026-08-05)"). Source tarball +`https://codeload.github.com/alexdobin/STAR/tar.gz/refs/tags/2.7.11b`, unpacked and read locally; +`source/…:line` locators below are into that tree, and are stable for that tag. `doc/STARmanual.pdf` +ships in the same tarball and is cited by section and page. `source/parametersDefault` is STAR's own +authoritative default-and-description table, from which the manual's parameter appendix is generated. + +Two secondary sources appear, both for corroboration only and never to establish a fact: a reply +from the STAR author on a combined-genome question, +[alexdobin/STAR#748](https://github.com/alexdobin/STAR/issues/748), and one peer-reviewed +measurement of combined-reference error, quoted in §2. + +## 0. The flags this pipeline actually passes + +Read off the two modules a chimera would run through today. + +`map/star-umi` (`src/seqforge/workflows/map/star-umi.smk:445-453`): + +```text +--runMode alignReads --genomeDir --runThreadN N --genomeLoad LoadAndKeep +--readFilesIn --readFilesType SAM {PE|SE} --readFilesCommand samtools view +--readFilesSAMattrKeep All [--clip3pAdapterSeq …] --outFileNamePrefix

+--outSAMtype BAM SortedByCoordinate --limitBAMsortRAM --outSAMmultNmax 1 +``` + +`map/star` (`src/seqforge/workflows/map/star.smk:247-253`): the same shape minus the uBAM input and +minus `--outSAMmultNmax`, plus `--quantMode GeneCounts`. + +What matters for a chimera is what is **absent**, because every absent flag is a default: + +| flag | value in effect | consequence on a chimera | +|---|---|---| +| `--outFilterMultimapNmax` | **10** (default) | counted across all components — §2 | +| `--outFilterMultimapScoreRange` | **1** (default) | near-best cross-component hits enter `NH` — §3 | +| `--outSAMunmapped` | **None** (default) | no unmapped record reaches any BAM at all — §5 | +| `--chimSegmentMin` | **0** (default) | chimeric detection OFF; no cross-component supplementaries — §1 | +| `--outSAMattrRGline` | unset | no `@RG` line to rewrite — §6 | +| `--outSAMattributes` | **Standard** = `NH HI AS nM` | plus input tags kept by `--readFilesSAMattrKeep All` | +| `--outSAMmultNmax` | **1** in `star-umi`; **−1** in `star` | the split is blind in one module and not the other — §4 | + +That last row is the whole story of this document, and it means **the two modules do not behave the +same way under a splitter**. `map/star` writes every reported locus; `map/star-umi` writes one. + +One aside for whoever prices `map/star`'s twin later: that module gets its counts from STAR itself, +`--quantMode GeneCounts` (`src/seqforge/workflows/map/star.smk:250`), which emits a single +`ReadsPerGene.out.tab` over whatever GTF the index was built with — for a chimera, the merged +annotation, both components' genes in one table. Splitting *that* is a different operation from +splitting a BAM, and the map's constraint 3 (a matrix per component, against its own annotation) +does not currently say which it means for `map/star`. Out of scope here; flagged because it is not +in the map's "Not yet specified" list either. + +## 1. Can a proper pair span two components? + +**No. VERIFIED.** + +STAR builds alignments inside *windows*, and a window carries exactly one chromosome. A new window +records its chromosome once, `WC[iWin][WC_Chr] = mapGen.chrBin[aBin >> P.winBinChrNbits]` +(`source/ReadAlign_createExtendWindowsWithAlign.cpp:66`), and an align may only be merged into a +neighbouring window when the two agree on that chromosome — the guard appears twice, once per +direction: + +```cpp +flagMergeLeft = flagMergeLeft && (mapGen.chrBin[iBin>>P.winBinChrNbits]==mapGen.chrBin[aBin>>P.winBinChrNbits]); +flagMergeRight = flagMergeRight && (mapGen.chrBin[iBin>>P.winBinChrNbits]==mapGen.chrBin[aBin>>P.winBinChrNbits]); +``` + +(`source/ReadAlign_createExtendWindowsWithAlign.cpp:28` and `:47`; window extension in +`source/ReadAlign_stitchPieces.cpp:100` and `:107` carries the same guard.) `chrBin` maps a genome +bin to exactly one chromosome (`source/Genome.cpp:214`, `chrBin[ii]=ichr-1`) because +`--genomeChrBinNbits` guarantees "each chromosome will occupy an integer number of bins" +(`source/parametersDefault:68-69`). + +Both mates of a paired alignment live in one `Transcript`, distinguished by a mate gap in the exon +list, and `nMates` is 2 only when that gap is present (`source/ReadAlign_alignBAM.cpp:70-82`). One +`Transcript` has one `Chr`. The record shape follows mechanically: for a two-mate alignment `RNEXT` +is written as the alignment's *own* chromosome — + +```cpp +//6: next refID +if (nMates>1) { pBAM[6]=trOut.Chr; } else if (mateChrmaxScore < (intScore)(P.outFilterScoreMinOverLread*(Lread-1))) + || (trBest->nMatch < (uint)(P.outFilterMatchNminOverLread*(Lread-1))) // ReadAlign_mappedFilter.cpp:8-9 + → statsRA.unmappedShort++; unmapType=1; +``` + +with both thresholds at **0.66** (`source/parametersDefault:480-487`, "normalized to read length +(sum of mates' lengths for paired-end reads)"). For a 2×150 fragment the bar is 0.66 × 300 ≈ 198, +and a single 150 nt mate cannot clear it however cleanly it aligns. **INFERRED** (the arithmetic is +mine, the filter and the `Lread` definition are verified): a balanced-length PE fragment that is +genuinely half host and half contaminant is filed as `unmapped: too short` and, with +`--outSAMunmapped None`, is written nowhere. This is the same denominator effect +`docs/research/smartseq3-tn5-read-through.md` measured for the Tn5 read-through, arriving from a +different direction. + +An **index hop** is the same shape and the opposite outcome: an index hop delivers a *whole* +foreign fragment, both mates from the contaminant, under a host cell's barcode. Both mates align to +one component, the pair is proper, and nothing in the BAM distinguishes it from real contamination. +**INFERRED.** A splitter cannot see index hops; it routes them to the component they came from, +which is the honest answer but not a detection. + +### Chimeric detection, and what "off by default" implies + +`--chimSegmentMin 0`, "if ==0, no chimeric output" (`source/parametersDefault:690-691`). Manual §6 +p. 15: chimeric segments are precisely those that "belong to different chromosomes, or different +strands, or are far from each other" — i.e. STAR's chimeric machinery is *exactly* the machinery +that could emit a cross-component template, and it is switched off. When it is on and +`--chimOutType WithinBAM` is chosen (`source/parametersDefault:682-688`), chimeric alignments enter +the main BAM with `0x800` set on the supplementary records +(`source/ReadAlign_alignBAM.cpp:200`, `alignType` −11/−12/−13 documented at `:47-55`) and the +default `Junctions` mode writes `Chimeric.out.junction` instead. + +**The implication for the contract, stated plainly:** with chimeric detection off, the BAM contains +no record whose *own* alignment spans components. Every cross-component signal is therefore a +*multi-locus* signal — several alternative placements of the same read — never a *split-read* +signal. A splitter that reasons about "spanning" is reasoning about a set of alternative loci, not +about one fragment straddling a boundary. If a later effort ever turns chimeric detection on, the +routing rule needs a fourth branch and the `0x800` records need a home; nothing in the current +design anticipates that. + +## 2. `--outFilterMultimapNmax` across the whole chimera + +**Counted across the whole chimera. VERIFIED. The exposure is real and one-sided.** + +`multMapSelect` scans **every** window of the read, over the whole concatenated genome, and collects +every transcript within `--outFilterMultimapScoreRange` of the best: + +```cpp +for (uint iW=0; iWmaxScore + P.outFilterMultimapScoreRange) >= maxScore ) { … nTr++; } +``` + +and then + +```cpp +if (nTr > P.outFilterMultimapNmax || nTr==0) { return; } // :46 +… +} else if (nTr > P.outFilterMultimapNmax) { // ReadAlign_mappedFilter.cpp:15-17 + statsRA.unmappedMulti++; unmapType=3; +``` + +There is no chromosome, contig-group or component term anywhere in that comparison. +`source/parametersDefault:463-465` says the same in words — "maximum number of loci the read is +allowed to map to. Alignments (all of them) will be output only if the read maps to no more loci +than this value" — as does manual §4.1 p. 10. **9 host loci + 2 contaminant loci = 11 > 10, and the +read is discarded entirely, with `uT:A:3` / "mapped to too many loci" in `Log.final.out` +(`source/Stats.cpp:132-133`), where a `ce11`-only run would have kept all 9.** Confirmed. + +**The flag that changes it is `--outFilterMultimapNmax`.** Raising it is not free: the manual (§4.1 +p. 11) requires `--winAnchorMultimapNmax ≥ --outFilterMultimapNmax`, and `winAnchorMultimapNmax` +(default 50) "also controls the overall sensitivity of mapping: increasing it will change (improve) +the mapping of unique mappers as well, though at the cost of slower speed" — i.e. raising the cap to +protect the split changes the alignment of reads that have nothing to do with the split. That is a +reason to leave it alone, not a reason to raise it. + +### Characterising the exposure + +The shape of the exposure, all **INFERRED** from the mechanism above: + +- **It is one-sided and it is a loss, never a gain.** Adding a component can only add loci, so + `nTr` on a chimera is ≥ `nTr` on either component alone. A read can cross the cap that would not + have; no read can fall back under it. +- **Only reads already near the cap are exposed.** A read at 1 host locus is at no risk from any + plausible number of contaminant loci. The population at risk is reads with 9 or 10 host loci — + already the tail of the repeat distribution — that additionally hit the other component. +- **The two conditions are close to independent and both are small**, so their conjunction is very + small. For the `ce11_ecHT115` case specifically the second condition is doubly unlikely: a + *C. elegans* read multimapping at 9–10 loci is repeat-derived, and *E. coli* HT115 shares + essentially no repeat family with a nematode. **This is a judgement, not a number** — the + measurement that would replace it is §"What would settle this cheaply". +- **It is invisible in the output.** The discarded read appears only as a `Log.final.out` counter; + it is not in the BAM, so no splitter and no downstream metric can recover it. The only honest + handle is the delta in `% of reads mapped to too many loci` between a chimeric run and a + single-assembly run of the same cells. +- **It bites `map/star` and `map/star-umi` identically**, because the cap is applied before any + output flag is consulted. + +For scale from the literature, cited as a secondary source and not as our measurement: Choi et al., +*BMC Bioinformatics* 23 (2022), "Expression-based species deconvolution and realignment removes +misalignment error in multispecies single-cell data" +([PMC9063264](https://pmc.ncbi.nlm.nih.gov/articles/PMC9063264/)), report that on combined-reference +multispecies single-cell data, "all error in combined reference accounted for only 0.4–1.4% of total +reads, these reads were concentrated to few genes, leading to strong false signals", with +"13,000–300,000 fewer reads aligned to human genes in the combined reference than in the human +reference". That paper measures a different quantity — misassignment, not cap-crossing — but it is +direct evidence that the aggregate effect of a combined reference is around a percent while the +per-gene effect is not. + +## 3. `MAPQ` semantics + +**The convention is confirmed; "computed over the reported alignment count" is confirmed with one +correction; "already component-local" is FALSE at the boundary.** + +VERIFIED, `source/ReadAlign_alignBAM.cpp:276-283`: + +```cpp +MAPQ=P.outSAMmapqUnique; +if (nTrOut>=5) { MAPQ=0; } +else if (nTrOut>=3) { MAPQ=1; } +else if (nTrOut==2) { MAPQ=3; } +``` + +with `outSAMmapqUnique` defaulting to 255 (`source/parametersDefault:372-373`). The map's stated +convention — 255 unique / 3 for two loci / 1 for three or four / 0 for five or more — is exactly +this. The manual states it as a formula, §5.2.1 p. 11: "The mapping quality MAPQ (column 5) is 255 +for uniquely mapping reads, and `int(-10*log10(1-1/Nmap))` for multi-mapping reads. This scheme is +same as the one used by TopHat". The identical ladder appears in the SAM text writer +(`source/ReadAlign_outputTranscriptSAM.cpp:210-216`) and the splice-graph writer +(`source/ReadAlign_outputSpliceGraphSAM.cpp:61-67`), so there is one convention and three copies of +it, not three conventions. + +`nTrOut` is the *reported* count, and `writeSAM` passes the **full** reported count into every +writer even when the write is truncated: + +```cpp +auto nTrOutWrite=min(P.outSAMmultNmax,nTrOutSAM); // :168 — how many records get written +for (uint iTr=0;iTr 1` — placed at several loci, so no gene +owns them"). Consequences, chained: + +- A cross-component multimapper arrives at the counter as `NH > 1` and is routed to + `_multimapping`. **No component's matrix ever counts it.** This is a genuine safety property the + map does not currently claim, and it holds *today*, before any splitter exists. +- Symmetrically: a read that would have been `NH:i:1` on `ce11` alone and is `NH:i:2` on the chimera + is **lost from the host's matrix** — moved from a gene into `_multimapping`. That is the counting + face of the §3 margin-of-1 problem, and it is a systematic under-count of the host, proportional + to how much of the host's transcriptome has a near-hit in the other component. **INFERRED.** +- The frozen umite fixture gives the order of magnitude for how much traffic rides on `NH` at all: + 1 640 of 12 977 aligned read names, **12.6%**, carry `NH > 1` on a single-assembly run + (`src/seqforge/workflows/umite/count.py:70`). That is the population whose `NH` a chimera can + perturb. + +### `HI`, and what a splitter actually gets to iterate over + +`HI = iTrOut + --outSAMattrIHstart` (default 1) is an **output-order index**, not a stable locus id +(`source/ReadAlign_alignBAM.cpp:296`; `source/parametersDefault:346-347`). + +Grouping depends entirely on the output type, and the manual is explicit (§5.3 p. 13–14): + +- `BAM Unsorted`: "The paired ends of an alignment are always adjacent, and multiple alignments of a + read are adjacent as well." +- `BAM SortedByCoordinate`: **no such guarantee**, and this is what both seqforge modules write. The + coordinate sort is by genomic position (`source/BAMoutput.cpp:77-94`); the per-record tiebreak key + is `(iReadAll<<32) | (iTr<<8) | mate` (`source/ReadAlign_outputAlignments.cpp:200`), which orders + records *within* a coordinate bin, not across the file. + +**So a splitter iterating a seqforge BAM does not get bundles.** To see all of a template's +alignments it must either name-sort (a full extra pass and, per +`src/seqforge/workflows/map/star-umi.smk:399-403`, 2× peak disk that the module deliberately refuses +to pay) or hold a QNAME → components map across the whole file. + +### The fatal interaction: `--outSAMmultNmax 1` + +**VERIFIED, and it is the finding that contradicts the map.** `map/star-umi` passes +`--outSAMmultNmax 1` (`src/seqforge/workflows/map/star-umi.smk:453`); it is a settled module literal +(ADR-0022 via #256 decision 7, recorded at `src/seqforge/workflows/umite/count.py:67`). Therefore +`nTrOutWrite = min(1, nTrOutSAM) = 1` and **exactly one alignment record per template reaches the +BAM**, whatever `NH` says. + +The consequences for the split contract: + +1. **"A template whose alignments span more than one component" is not observable.** The splitter + sees one RNAME. `NH:i:3` tells it three loci exist; it tells it nothing about *where*. The map's + three-way routing — one component / ambiguous / unmapped — cannot be computed from this artifact. + A cross-component multimapper and a within-component multimapper are byte-indistinguishable + except by re-alignment. +2. **Which component the single record lands in is a tie-break, not a decision.** When + `--outSAMmultNmax != -1`, `multMapSelect` partitions `trMult` so that top-scoring alignments come + first (`source/ReadAlign_multMapSelect.cpp:62-68`) and marks `trMult[0]` primary (`:87-88`). Among + equal-best alignments the survivor is the first in *window order* under the default + `--outMultimapperOrder Old_2.4` (`source/parametersDefault:282-285`; manual §5.2.1 p. 11, "the + order of the multi-mapping alignments for each read is not truly random"). Window order is the + order STAR created windows while scanning the read's anchor seeds, "going through ordered + positions in the suffix array" (`source/ReadAlign_stitchPieces.cpp:40-49`) — **suffix-array order, + not coordinate order and not component order**. It is deterministic for a given index and read, + and it is not something a splitter, a concatenation order or a flag can steer. **So for a read + that ties exactly across components, which component's BAM it lands in is arbitrary in the precise + sense that nothing downstream can predict or control it.** The repo already records the tie-break + half of this at `src/seqforge/workflows/__init__.py:287-292`. +3. **Every written record is primary.** `trMult[0]->primaryFlag=true`, so no `0x100` survives and + `cram.py`'s `-F 0x100` (`src/seqforge/workflows/cram.py:148`) is the cheap invariant its docstring + says it is. A splitter must not expect secondaries to exist. +4. **`HI` is always `HI:i:1`** on this module's BAM, since only `iTrOut == 0` is written. It carries + no information a splitter can use. +5. **`map/star` does not share this**, because it omits the flag: its BAM carries every reported + locus (up to 10), so spanning *is* observable there — after a name sort or a QNAME pass. + +The map's constraint 4 already forces `map/star-umi`'s chimeric twin to be a standalone `.smk` copy. +**INFERRED, and this is the cheap way out:** that copy is free to drop `--outSAMmultNmax 1`, or to +set it to `--outSAMmultNmax -1`, at the cost of a larger sort and a larger BAM (the repo measured +~18% of the sort spent on records `-F 0x100` later discards, `src/seqforge/workflows/__init__.py:280-282` +— so dropping the flag re-incurs roughly that). It cannot be dropped *silently*: it changes +`workflow_version`, hence `run_id`, and the counter's `_multimapping` behaviour is defined against +`NH` rather than bundle length precisely so the flag can move without moving the counts +(`src/seqforge/workflows/umite/count.py:66-74`). That last property is what makes this fixable +rather than fatal. + +## 5. `--outSAMunmapped`, and whether "unmapped" is component-attributable + +**It is not, by construction. VERIFIED.** + +First, the state of the world: **no seqforge module passes `--outSAMunmapped`**, so the default +`None` — "no output" (`source/parametersDefault:349-353`) — is in force and **no unmapped record +reaches any BAM seqforge writes today.** The map's routing branch "an unmapped read stays as STAR +left it" is, against the current modules, a branch over an empty set. The unmapped population is +visible only in `Log.final.out`. + +If a chimeric twin turns it on, this is what it gets. An unmapped record is written with: + +```cpp +if (alignType<0) { pBAM[1]=trOut.Chr; } else { pBAM[1]=(uint32) -1; } // refID — :521-525 +if (alignType<0) { pBAM[2]=…; } else { pBAM[2]=(uint32) -1; } // POS — :527-532 +``` + +`MAPQ` left at its initialisation of 0 (`source/ReadAlign_alignBAM.cpp:116`), and + +```cpp +attrN+=bamAttrArrayWriteInt(0,"NH",…); attrN+=bamAttrArrayWriteInt(0,"HI",…); +attrN+=bamAttrArrayWriteInt(trOut.maxScore,"AS",…); attrN+=bamAttrArrayWriteInt(trOut.nMM,"nM",…); +attrN+=bamAttrArrayWrite((to_string((uint) alignType)).at(0), "uT",…); +``` + +(`source/ReadAlign_alignBAM.cpp:156-161`). So: **`RNAME = *`, `POS = 0`, `MAPQ = 0`, `NH:i:0`, +`HI:i:0`, plus a `uT` tag giving the reason** — 0 no seed/window, 1 too short, 2 too many +mismatches, 3 too many loci, 4 unmapped mate of a mapped pair (manual §5.2.2 p. 13 and §5.4 p. 14; +`source/ReadAlign_mappedFilter.cpp:4-18` sets 0–3 and +`source/ReadAlign_outputAlignments.cpp:213` sets 4). In the coordinate-sorted BAM every such record +goes to the final bin — `if (bamIn32[1] == ((uint32) -1)) { iBin=P.outBAMcoordNbins-1; }`, +`source/BAMoutput.cpp:89-90` — so they are appended at the end of the file, ordered by read number +(`source/BAMbinSortUnmapped.cpp:16-48`). + +**There is no field on such a record that names a component**, and there could not be: STAR did not +place the read. **The contract should say this explicitly**, and should say what it does with the +pile — an unmapped record is not evidence of *either* organism, and duplicating it into every +component's BAM would double-count it while dropping it would lose it. Neither the map nor this +document decides that; the map's "Not yet specified" list should grow a line for it. + +**One wrinkle a naive rewrite will miss, VERIFIED.** For `uT:A:4` — the unmapped mate of a mapped +read — STAR fills `RNEXT`/`PNEXT` from the *mapped* mate: `mateChr=trOut.Chr; mateStart=…` +(`source/ReadAlign_alignBAM.cpp:131-134`), reaching `pBAM[6]=mateChr` at `:551-552`. So the record +is `RNAME = *` but `RNEXT = __`. **An unmapped record can carry a chimeric +chromosome name in `RNEXT`.** Any `@SQ`/RNAME rewrite that walks RNAME only will leave a dangling +reference id behind and produce a BAM that fails `samtools quickcheck`-grade validation against the +rewritten header. + +## 6. What besides `@SQ` names the reference + +VERIFIED, entirely from `source/samHeaders.cpp` — one function builds the whole header. + +| line | what it emits | names the chimera? | +|---|---|---| +| `:29-31` | `@SQ SN: LN:` for every real chromosome | **yes** — this is the `__` set | +| `:35-54` | extra `@SQ` lines read verbatim from `/extraReferences.txt` | **yes**, if the index was built with on-the-fly insertions | +| `:56-62` | optional extra `@PG` from `--outSAMheaderPG` | only if the user puts it there | +| `:64` | `@PG ID:STAR PN:STAR VN: CL:` | **yes — this is the one that hurts** | +| `:66-76` | `@CO` lines slurped from `--outSAMheaderCommentFile` | only if the user puts it there | +| `:79-81` | `@RG` from `--outSAMattrRGline` | no; and seqforge passes none, so there is no `@RG` at all | +| `:84` | `@CO user command line: ` | **yes** | +| `:86` | `P.samHeaderExtra` | no | +| `:88-100` | `@HD VN:1.4`, plus `SO:coordinate` for the sorted BAM | no | + +`commandLineFull` is not the argv — it is STAR's *final effective* command line, rebuilt from every +parameter whose `inputLevel > 0`, i.e. every parameter the user set from any source +(`source/Parameters.cpp:446-456`, logged as "Final effective command line"). `commandLine` is the +verbatim argv (`source/Parameters.cpp:332-362`). **Both contain `--genomeDir `**, and `commandLineFull` additionally contains every other flag the run used. + +So a "restore the header a single-assembly run would have produced" clause has exactly three lines +to deal with, and they are not equal: + +- **`@SQ`** — must be rewritten. Already in the map. +- **`@PG … CL:`** — names the chimera index in a free-text field. Rewriting it fabricates a command + that was never run; leaving it makes the per-component BAM self-describing as chimera-derived, + which is arguably *correct* provenance. **This is a real fork in the contract and the map does not + mention it.** +- **`@CO user command line:`** — the same fork, one line down. + +**INFERRED, offered as a recommendation and not a decision:** leaving `@PG`/`@CO` alone and adding a +`@PG` of the splitter's own (with `PP:STAR`, which is what the SAM spec's `PP` chain is for) is both +honest and cheaper than rewriting, and it means the byte-for-byte "identical to a single-assembly +run's header" bar in map constraint 6 cannot be met literally. **That bar needs restating as +"identical `@SQ` block and `@HD`" or it will fail on its first synthetic round-trip** — including on +the `tinyEcDub` fixture, whose whole point is to break naive assumptions. + +Two further header facts worth having: the binary BAM reference dictionary is written from +`chrNameAll`/`chrLengthAll`, which is `chrName` plus any extra references +(`source/samHeaders.cpp:33-34`, `:106`; `source/BAMbinSortByCoordinate.cpp:63`), so the text header +and the binary dictionary must be rewritten together. And STAR ships its own *narrow* answer to +"extra sequences in one index" — `--outSAMfilter KeepOnlyAddedReferences` / +`KeepAllAddedReferences` (`source/parametersDefault:396-399`), implemented as an index-range test +`trOutSAM[itr]->Chr < mapGen.genomeInsertChrIndFirst` +(`source/ReadAlign_outputAlignments.cpp:141-164`). It is a *filter on one output*, not a split into +N, it keys on **index contiguity rather than on names**, and it applies only to sequences inserted +at mapping time with `--genomeFastaFiles`, not to a genome built chimeric. **It is not usable here** +— but it is confirmation that STAR itself has no name-based notion of a component, which is why the +map's constraint 1 (names belong to `liulab-genome`) has nothing on the STAR side to consume it. + +## 7. Verdict on the map's standing decision 2 + +The decision reads: *"with the spanning templates removed, each component's BAM needs an `@SQ` and +RNAME rewrite and no MAPQ or NH surgery at all, because STAR only reports best-scoring loci — a read +that beats its cross-species hit already carries the `NH` and `MAPQ` a single-assembly run would have +given it."* + +Graded clause by clause: + +| clause | verdict | +|---|---| +| a template's alignments all land in one component ⇒ no MAPQ/NH surgery needed | **TRUE, and for a stronger reason than stated.** MAPQ and NH are pure functions of `nTrOutSAM`. If every reported locus is in one component, `nTrOutSAM` is what a single-assembly run would have computed, so both fields are already right. Nothing to fix. | +| "STAR only reports best-scoring loci" | **FALSE as written.** It reports loci within `--outFilterMultimapScoreRange` (default 1) of the best. A read beating its cross-species hit by exactly 1 is still `NH:i:2 MAPQ:3`. | +| the routing rule is computable from the BAM | **FALSE for `map/star-umi` as it stands** — `--outSAMmultNmax 1` writes one record per template, so "spans more than one component" is unobservable. TRUE for `map/star`. | +| `@SQ` + RNAME rewrite is the whole header job | **INCOMPLETE.** `@PG CL:` and `@CO` embed `--genomeDir`. And `RNEXT` on a `uT:A:4` record names a chromosome too. | +| "an unmapped read stays as STAR left it" | **VACUOUS today** — `--outSAMunmapped` defaults to `None` and no module passes it, so there are no unmapped records in any seqforge BAM. | + +**Standing decision 2 survives as a design intent and does not survive as a statement about the +artifact `map/star-umi` produces.** The property it wants — no arithmetic on the species BAMs — is +real and rests on a mechanism (MAPQ and NH are functions of the reported-locus count) that is more +robust than the reason given for it. What fails is the *premise* that a splitter can see the loci it +needs to see. The cheapest repair is a one-flag change in the chimeric twin the map already requires +to exist, and the counter is already insulated from that flag by design. The decision should be +amended rather than reversed. + +## What is INFERRED and not verified, collected + +Kept in one place so nobody has to re-derive which is which: + +1. That a balanced-length PE fragment split across components is filed `unmapped: too short` rather + than as a one-mate alignment. The filter and `Lread` are verified; the arithmetic that a 150 nt + mate cannot clear 0.66 × 300 is mine and is unobserved. +2. That an index hop presents as an ordinary within-component proper pair and is undetectable. +3. That the `--outFilterMultimapNmax` exposure is negligible for `ce11_ecHT115` specifically. The + mechanism is verified; the magnitude is a judgement. +4. That the margin-of-1 population — reads beating a cross-species hit by exactly one point — is + what makes "no surgery" and "correct" come apart, and its size. Unmeasured in either direction. +5. That the systematic host under-count via `NH:i:1 → NH:i:2` is real. Follows from verified + mechanics; never observed. +6. That dropping `--outSAMmultNmax 1` in the chimeric twin costs roughly the ~18% sort overhead the + repo measured for the reverse change, and nothing else. +7. That leaving `@PG`/`@CO` alone and chaining a splitter `@PG` with `PP:STAR` is the better of the + two forks. A recommendation, not a measurement. +8. That `--genomeChrBinNbits` needs scaling for the 94-sequence `ce11_ecHT115` index (manual's own + guidance, `min(18, log2[max(GenomeLength/NumberOfReferences, ReadLength)])`). This is a memory + and build concern belonging to `liulab-genome`, and it does not touch correctness of the split. + +## What would settle the inferred claims cheaply + +All of it fits inside the map's existing bar (constraint 6, the synthetic round-trip on +`tinyCe`/`tinyEc`/`tinySc`/`tinyEcDub`), and none of it needs a benchmark re-run: + +- **Claims 1 and 2**: synthesise a fragment with mate 1 from `tinyCe` and mate 2 from `tinyEc`, run + the twin with `--outSAMunmapped Within`, and read the `uT` tag. One cell, one read. +- **Claims 4 and 5**: plant a read that matches `tinyCe` exactly and `tinyEc` with one mismatch, and + read `NH`. Confirms or refutes the score-range boundary in a single record. +- **Claim 3**: on the real plate, diff `% of reads mapped to too many loci` in `Log.final.out` + between the existing `ss3-ce11-ws298-ce321d3fc6d4` run and a chimeric run of the same cells. The + numbers are already on disk for one side. +- **Claim 6**: measured by the twin's first real run; no separate experiment. + +Claims 7 and 8 are decisions and a build parameter respectively, not measurements, and belong to the +map and to `liulab-genome`.