Changelog
All notable changes to this project will be documented in this file.
[0.9.2] - 2026-08-25
Added
-
A WPS panel-anchor figure.
.WPS.panel.parquetwas written on every panel run and plotted by nothing, so the report showed only the genome-wide anchors. On a targeted assay those are two different measurements: one XS2 sample gave amplitudes of 1.6 (CTCF) and 14.7 (TSS) over 34,851 and 23,671 genome-wide anchors, against 262 and 199 over the 105 and 61 the panel targets. Both are now drawn, from one shared builder, and each legend states the anchor count behind its mean. -
scripts/lint.sh— runs the lint suite at the versions CI runs it.
CI installs .[dev], so it uses the == pins in pyproject.toml; a
developer's shell uses whatever it happens to have. On one machine that was
black 26.5.1 against a 26.1.0 pin and mypy 2.3.1 against 1.19.1 — and 1.19.1
reports [import-untyped] where 2.3.1 does not, so #71 passed locally and
failed CI on exactly that. The script reads the pins from pyproject.toml
and runs each tool through uvx at that version, keeping one source of
truth. Without uvx it still runs but prints the version it used, since a
silently different version is the whole problem.
tests/unit/test_lint_script_matches_ci.py asserts it stays faithful — the
same tools and the same flags as the workflow, versions never hardcoded, no
pinned tool left unrun. A script trusted to reproduce CI while quietly
checking less is worse than no script.
pip install 'krewlyzer[all]'— every optional capability in one target:[report](HTML rendering and charts) pluspsutil, which sharpens available-memory detection. It excludes[docs],[test]and[dev]on purpose; installing the tool should not deliver mkdocs and a linter.
tests/unit/test_optional_dependencies.py comes with it, and is the more
useful half. It walks the source for imports guarded by except ImportError
and asserts each is declared in some extra — which is the check that was
missing when markdown was imported and declared nowhere. It found psutil
in the same state while being written. It also pins that [all] covers every
runtime extra, and that it covers no contributor tooling.
Fixed
capture_maskwas computed against one chromosome's baits, for every chromosome. Neither call site inrust/src/wps.rssupplied the region's chromosome: one hardcoded0u32under a comment reading "Simplified - need proper chrom_id lookup", the other used*self.trees.keys().find(|_| true), an arbitrary key from HashMap iteration order and so not even deterministic. Every position on every chromosome was tested against whichever chromosome that resolved to.
Measured on an XS2 sample: of 33 genome-wide anchors flagged as captured, all 33 overlapped a chr1 bait while only the 9 that were themselves on chr1 overlapped a bait on their own chromosome — of the 24 elsewhere, zero did. In the panel table all 5 flagged anchors were on chr1, and all 5 chr1 anchors were flagged. The mask produced both false negatives (real on-bait positions elsewhere marked unreliable) and false positives (off-bait positions colliding with a chr1 bait).
BaitMask is now keyed by normalised chromosome name instead of a numeric
id, which deletes the lookup both callers were getting wrong rather than
fixing it in two places. Four Rust tests cover it, with baits at
deliberately colliding coordinates on different chromosomes; two fail
against the old behaviour.
Outputs produced before this fix carry the wrong mask. No feature value is affected — nothing else read the column — but any bait-coverage figure derived from it understates real coverage for anchors off chr1. Re-run WPS to correct it.
The bait-coverage line described in the next entry is computed from this column, so on pre-fix outputs it reports the wrong figure.
-
WPS figures plotted bin index but labelled the axis "bp". Every WPS profile is 200 bins of 10 bp over a 2000 bp window (
rust/src/wps.rs), and all three panels used the bin number directly as the x value. Distances read ten times too small: a TSS dip at bin +15 appeared at "+15 bp" when it is +150, and the ±1000 bp window appeared as ±100. Nothing was miscomputed — the curves match the Parquet exactly — which is why it survived: the figure was wrong only in a way a reader takes at face value.WPS_BIN_BPis pinned invalidate/claims.pyagainstNUM_BINSandWINDOW_SIZE, so changing the window fails a test rather than silently rescaling every WPS axis again. -
WPS figures did not say how much of the profile the assay captured.
capture_maskmarks bins inside a bait and away from its edges, and nothing in the Python package read it — not the plots, not validation. On a targeted panel the genome-wide anchors are almost entirely off-bait: one XS2 sample had 506 captured bins of 11,704,400. The profile is real, and off-target coverage genuinely carries fragmentomic signal, but a clean chromatin curve invites a reading it cannot support.
Both anchor panels now carry a ±1 SEM band and report bait coverage per anchor as well as per bin. The band is the honest way to convey what a mean is worth: a profile over 34,851 anchors and one over 105 are otherwise drawn identically, and on one sample the genome-wide CTCF band is 26% of the curve's own amplitude while TSS's is 6%. That says "noise-dominated" without the caption having to, which matters because deciding what a measurement supports is the downstream consumer's call, not this toolkit's. An earlier draft ended the caption "read this as a fragmentomic profile, not a targeted measurement" — a usability verdict — and a test now pins that such wording stays out.
Both panels also report bait coverage per anchor as well as per bin. The bin fraction alone misleads: a WPS window is 2000 bp against baits of order 100 bp, so even a perfectly bait-centred anchor captures a small minority of its bins — the best anchor on that sample reached 103 of 200 and the median among capturing anchors was 7. A low percentage is the expected geometry, not a badly chosen anchor set. What the anchor count does separate is real: 5 of 166 panel anchors touch a bait against 33 of 58,522 genome-wide. Counts rather than a bare percentage throughout, because every fixed-decimal format renders 0.00432% as "0.00%", which reads as rounding rather than as a fact.
describe-output -o page.htmlwrote unrendered Markdown. The file was a valid HTML document — doctype, head, styles — whose body was the Markdown source escaped into a<pre>, so it displayed##headings and|tables as literal text. Nothing said why, and the CLI reportedwrote <path>.
render_html_page converts the Markdown only when the markdown package is
importable, and markdown was declared in no dependency group at all —
not core, not report, not dev. No install of krewlyzer could satisfy it,
so every user took the fallback; it worked only where something unrelated
(mkdocs) had pulled the package in.
markdown>=3.4 now sits in the report extra, alongside plotly — that
group is human-readable output, which is what both commands produce. The
fallback also announces itself now, in the page and on the console, matching
krewlyzer report, whose charts each state when plotly is absent rather than
quietly rendering nothing. Asking for .html and getting Markdown is
confusing precisely because the file is HTML; nothing looks wrong except
the contents.
KREWLYZER_DATA_DIRis now honoured by every bundled-data lookup, not just one. The README prescribes it for pip installs — the wheel excludesdata/to stay under PyPI's 100 MB limit, so those users keep a separate clone and point at it.AssetManagerhonoured it;get_bundled_gene_bedandbuild-gc-referencedid not, and each resolved against the installed package instead.
For a pip user, that meant panel assets loaded correctly while
get_bundled_gene_bed returned None — and None is not an error here, so
gene-level FSC and region-MDS were simply absent from the output with only a
debug-level log to say so. build-gc-reference degraded more loudly,
warning that regions "may include problematic areas", but degraded all the
same.
The resolution rule now lives once, in assets.bundled_data_dir(). Where the
variable is set it is the whole answer — resolve there or return nothing,
never fall back to the bundled copy, matching what AssetManager already
did. A fallback would demote the override to a suggestion and mask a
misconfigured data clone behind whatever the package happened to ship.
Nothing caught this because the tests that would have were skipped in CI for
precisely the same reason: conftest.DATA_AVAILABLE asked the dataless
installed package too. It now resolves the way the runtime does, which puts
16 previously-skipped tests into every CI run — and they found this on the
first one.
[0.9.1] - 2026-08-15
Added
validate-output --expect— reconcile the samples you meant to run against the ones that exist. Everything else in the gate validates the samples it finds, and two failure modes leave nothing to find:
| before | |
|---|---|
| No output directory at all | invisible |
| Directory exists, no Parquet written | invisible — discover_samples skips it |
| Parquet present, no completion marker | ERROR ✓ |
A job killed between mkdir and its first write was indistinguishable from
one never submitted. At a handful of samples that is obvious; across a
16,000-sample cohort it is not, and the consumer reads a short cohort without
knowing it was short.
Takes a Nextflow samplesheet (using its sample column) or a plain list. The
two absences are reported separately, because an empty directory means the
job ran and its logs exist while a missing one usually means it was never
submitted. An unexpected sample is a WARN, not an ERROR — it corrupts
nothing, but it usually means an identifier that does not round-trip, which
makes every other count unreliable. A CSV without a sample column is an
error rather than a guess at the first column: reconciling against the wrong
field reports every sample missing, which reads as a catastrophic run failure
instead of a malformed input.
The summary table shows counts; the --json-report names every affected
sample, which is the split that keeps it readable at cohort scale.
- A check that a release actually published. Every other gate verifies the repository; nothing asked the registries whether the release landed. 0.9.0 passed every workflow while the docs site served 0.8.3 and no GitHub Release existed — both invisible, because green CI only means the configured steps ran, not that the steps you need exist.
scripts/check_release_artifacts.py compares the newest tag against PyPI,
GHCR, the docs stable alias, and the Release object. It runs weekly, and
after the Release workflow finishes via workflow_run — not a tag trigger
with a sleep, which would idle a runner for 40 minutes and cry wolf
whenever publishing ran long. The weekly schedule is the load-bearing one: a
post-release check alone would not have caught the docs, because that deploy
had been deleted three weeks before 0.9.0.
Three exit codes, because "could not reach the registry" is not "in sync": 0 matches, 1 stale, 2 unreachable. A network failure reported as success is the silent pass the whole script exists to prevent.
-
check_phi_guard.shassertion 6 — every tracked hook exists, matches the staged version, andpre-pushstill delegates to git-lfs. Compared against the index rather thanHEAD, so a deliberate edit passes while an out-of-band overwrite fails. -
A stale-extension warning at pytest session start. Python imports whatever
maturin developbuilt last, so after arust/src/change the suite tests the previous binary — passing, and proving nothing. Inconftest.pyrather than a per-tool hook so it fires for a human, for CI, and for any agent. A warning and never a failure: the comparison is mtime-based andgit checkoutrewrites mtimes, so a branch switch can trip it. -
release.ymlnow creates the GitHub Release, using the CHANGELOG section for the tag as the body. Until now the Releases page was hand-maintained, which is why 0.9.0's tag went green with nothing on it. A missing section fails the step; a section over 60,000 characters is cut at a heading boundary and linked (0.9.0's was 103,376 against GitHub's 125,000 limit). -
MIN_PON_VERSIONandKREWLYZER_ALLOW_OLD_PONpinned in the claims registry, now that three documents quote them.
Changed
- FILTER_MAF no longer runs once per sample when the MAF is shared. It
scanned the entire cohort MAF, line by line in Python, to extract one
sample's rows — so a 16,552-sample run performed 16,552 full scans of the
same file, each as its own SLURM job whose scheduling overhead exceeded the
work it did.
maxForks queue_size.intdiv(2)then reserved 250 of 500 slots for it, upstream of RUNALL, so the opening wave of a run was nothing but MAF filtering and no real work could start. This was the observed queue clog on the 0.8.3 cohort.
The new SPLIT_MAF does one pass per distinct MAF: 1 job instead of
16,552, and the whole queue freed for RUNALL.
Batching keys on (maf, single_sample), not the MAF alone — grouping a pass-through sample with a filtered one would apply a single mode to both and silently corrupt one. A sheet with per-sample MAFs yields one task each and behaves exactly as before; samples with no MAF keep FILTER_MAF and its header-only placeholder.
Every requested sample gets a file, including those with zero matching
variants, and SPLIT_MAF exits non-zero if it wrote a short set. A sample
with no file would be dropped by the caller's re-pairing and RUNALL would
never run for it — silent cohort loss reintroduced one stage earlier than
--expect can see it. The metas ride through the process so re-pairing is a
map, preserving the join-free streaming this workflow is built around.
- Reverted the RUNALL retry ladder added moments earlier in this cycle. It
was built on the claim that the pipeline had no retry policy. That claim was
wrong: queue, errorStrategy and maxRetries come from the institutional
config -profile iris loads over HTTPS at parse time, so they appear nowhere
in this repository and a grep finds nothing. iris sets maxRetries = 3,
retries any failure, falls back to ignore, and selects the partition with a
closure honouring --isolated and --partition.
Because withName outranks a generic process block, the ladder silently
replaced all of that — cutting retries from 3 to 1, ignoring --partition
and --isolated, and pinning attempt 1 to the smaller of two partitions. The
task.attempt ramps it claimed to revive had been live all along.
The pre-existing behaviour is better and needs no queue names: a growing time request makes SLURM pick the partition, since attempt 2's 4h no longer fits a 3h queue.
tests/unit/test_runall_resources.py now pins the absence of queue,
errorStrategy and maxRetries from that block, so the next reader who
greps, finds nothing and "fixes" it lands on the explanation instead.
- check_phi_guard.sh runs in 2.3s instead of ~100s. git archive applies
the smudge filter, so assertion 1 was unpacking and scanning 660 MB of
parquet, gzipped BED and FASTA (75s), and assertion 5 grepped the same
binaries for /Users/ (13s). GIT_LFS_SKIP_SMUDGE exports the 132-byte
pointers instead, and assertion 5 skips LFS paths. This gives up scanning
inside LFS payloads: accepted, because they are machine-built compressed
reference data where a plaintext regex could never have matched.
-
Release guide Phase 5 merges
mainintodevelop, not the release branch. The old form left each branch holding a merge commit the other lacked, somainsat permanently one ahead; mergingmainmakes it an ancestor andgit merge-base --is-ancestor main developbecomes a valid post-release check. -
The gates table in
AGENTS.mdhas a "Read first" column..agents/rules/was described as read "on demand" with nothing creating demand, so those files went unopened for a whole release.
Fixed
-
A 16,552-sample run produced zero samples. Three faults compounded:
-
SPLIT_MAF interpolated the sample ids into its script. Nextflow strips a script block's common indentation after interpolating, so 16,552 zero-indent id lines left no common prefix, nothing was stripped, and python saw
import oson line 2..command.shwas 16,648 lines: 96 of code, 16,552 of ids. Exit 1. The ids are now a staged file — a filename is one token and cannot defeat the dedent at any cohort size. - Retries could not submit. The scheduler runs
EnforcePartLimits=ALL, where a partition list enforces the intersection of limits and offers the union of nodes.2.h * task.attemptasked 4h on attempt 2 against a list containing a 3h partition and was refused — "Requested time limit is invalid" — not rerouted. New--long_partitionsends retries to a partition whose limit covers the grown request. ignoreon a fan-in step lost everything. SPLIT_MAF is one task for the whole cohort, so ignoring it left RUNALL with no input and the run completed having done nothing. It is nowterminate: a fan-in failure stops the run in ninety seconds instead of after the queue drains.
The unit tests could not have caught (1): the harness dedents before
substituting, Nextflow does the reverse, so it tested well-formed Python while
production shipped broken. The stub run could not either — its ids are joined
with spaces, one line, nothing to defeat. test_module_interpolations.py now
refuses any newline-joined value in a .nf.
- Nextflow 26.04 could not parse
nextflow.configat all. The nf-coretry { includeConfig } catchidiom is rejected outright — "Try-catch blocks cannot be mixed with config statements" — so the pipeline refused to start. 25.10.3 accepted the same file, which is how it went unnoticed: the only signal was a version nobody had run. Replaced with the ternary nf-core moved to, which is why their own configs already parse under 26.
The behaviour change is the better half. The old catch swallowed any
failure, so a transient network blip left the run going with no institutional
config — no process_single label, no queue closure honouring
--isolated/--partition, no maxRetries — and three processes silently
falling back to 8 CPU / 32 GB / 24h, a request a short partition refuses at
submission. A whole run's resource policy changed, reported as one line of
stderr. Skipping now requires the operator to set NXF_OFFLINE.
tests/unit/test_nextflow_config_parses_on_26.py pins the construct. It is
not a parser and cannot claim the file is 26-clean — the parser stops at the
first error, so more may sit behind any given one. Only a real
NXF_VER=26.04.6 ... -stub-run can establish that.
- SPLIT_MAF's stub block still referenced the pre-rename input. Its input
became
metasso the metas could ride through and the caller could re-pair with amapinstead of ajoin; the script block and tag were updated and thestub:block was not. Groovy resolves${...}at task runtime, so the file parsed, every unit test passed, andnextflow -stub-runfailed immediately withNo such variable: sample_ids.
tests/unit/test_module_interpolations.py now asserts, for every module,
that the stub block references no name absent from its inputs and script —
which is exactly the shape of a rename applied to one block and not the
other. Written as a subset comparison between the two blocks rather than a
declaration check: parsing Nextflow's input syntax precisely (path(x), bare
path x, multi-element tuples) produced false positives on seven modules,
and a test that cries wolf is worse than no test.
- The Nextflow pipeline reported itself as 0.8.3.
manifest.versionsat at the previous release through all of 0.9.0. Nothing functional depended on it, which is why it drifted — but Nextflow writes it into every execution report, trace file and Tower entry, so a 16,000-sample cohort would have been labelled 0.8.3 in its own provenance while running 0.9.0 containers.
The Phase 2 version bump is a sed over container tags in the modules and
never matched this line. The release guide compounded it by listing
nextflow.config under "container tag" — there is no container pin in that
file at all. Both corrected, and tests/unit/test_nextflow_version.py now
pins manifest.version and every module container tag to __version__.
The container half matters more than the manifest: a stale pin there does not mislabel a run, it executes the previous release's defects under this version's name.
- The documentation site stopped publishing on 2026-07-31 and nobody noticed
for a whole release.
ba13fd4(#32) replaced the docs workflow with the aggregate check alone, dropping themikedeploys and the strict build with it. Nothing has published since:stablestill served 0.8.3 — the documentation for the release whose defects 0.9.0 exists to fix — anddevfroze at the same date. Tagging 0.9.0 did not help, because the old workflow only deployed on pushes tomain/develop, never on tags.
Restored, with three changes. A tag publishes stable, not a push to
main: the tag is the authoritative "released" signal, and hanging it on a
main push is what let PyPI and the container ship while the site sat still.
No paths filter on push, because GitHub requires the ref filter and the
path filter to both match, so tags: ['*'] beside paths: would drop any tag
whose commits missed docs/. And workflow_dispatch takes a version, so
0.9.0 can be backfilled without re-tagging.
mkdocs build --strict returns as a PR check. Its absence is why three links
to non-existent anchors survived into the release.
- Git LFS had silently disabled the PHI guard.
git lfs installwrites its hooks intocore.hooksPath—.githooks/here — and only declines when it recognises the hook already present. Checking out a commit predating.githooks/removes the hooks, and the next LFS command installs a three-linepre-pushstub in their place. That happened during the 0.9.0 release:pre-commitwas gone andpre-pushwas the stub, so nothing scanned staged content for patient identifiers, and nothing reported it.
pre-push now reads the ref list once and delegates to git lfs pre-push
after its own checks pass, so LFS never needs to install anything and the two
cannot clobber each other. The three LFS hooks are tracked for the same
reason. A side effect: git push uploads LFS objects by itself again, rather
than needing a separate git lfs push --all.
-
The version guard hardcoded
0.9.0in one of its two messages instead of formattingMIN_PON_VERSION, so raising the floor would have left that branch naming the old one. -
Three documentation links pointed at anchors that never existed, including one whose target is an
!!! noteadmonition, which generates no anchor at all.
Documentation
-
The PON version guard is in the README. 0.9.0 refuses any PON built below the floor — the first thing an upgrading user with a self-built model meets — and neither the refusal nor its escape hatch was documented there.
-
build-pon --from-outputsis in the README andpon-models.md, not just the developer guide. It turns a ~15h build into minutes and is how all four bundled models were rebuilt for 0.9.0. Both examples previously shown for it used a--pon-variantflag thatbuild-pondoes not have. -
BreakPointMotifadded to the README feature table.
[0.9.0] - 2026-08-11
Added
- The four GRCh37 PONs re-aggregated with the σ fix and the breakpoint
blocks, and stamped 0.9.0. Rebuilt via
--from-outputsfrom the samerun-alldirectories, so all four cohort digests are unchanged —1e7158a60f3e3ca7,a0fed96df4c948ae,5c8c873ec1def225,655a5962b075cea2— and 22 blocks become 24.
The rebuild is provably a NaN-ing of residue and nothing else. Against the previous build: identical anchor sets, bit-identical means, bit-identical σ at all 24.7 M surviving positions, 1,177,647 positions newly NaN whose old σ topped out at 7.9 × 10⁻¹³, and zero positions newly finite.
Residue σ falls to zero in three models and to 2 positions of 15.8 M in
xs2.all_unique (σ 2.2 × 10⁻⁷ against a mean of 9.7 × 10⁻⁸ — the same
phenomenon two decades above the 10⁻⁹ floor, left rather than chased with a
reactive threshold change).
validate-pon passes all four with zero findings — the first clean pass.
On a real XS1 plasma sample, BreakPointMotif median |z| falls from 5.85 to
2.36, max from 159 to 8.5, and motifs beyond |z| = 10 from 70 to none,
putting it alongside EndMotif at 1.82 rather than 3× worse.
- The product now records what produced it.
{sample}.metadata.parquet— the consumer's completion marker — gains five columns:krewlyzer_version,pon_applied,pon_model,pon_cohort_digestandpon_krewlyzer_version.
The version was recorded, in {sample}.features.json, which invariant #2
says downstream never reads. So the one fact needed to tell two runs apart
lived precisely where the product ignores it. That mattered the moment this
release changed what several columns mean: given a results directory there
was no way to answer was this produced before or after the FSD bin fix?
except by re-running it.
pon_model is the basename only. The full path would put the operator's
home directory into a shipped table.
Stamped by run_features, so a single-feature CLI run records the same
provenance as run-all (invariant #6), and stamped from pon_for_zscore
rather than pon — a model loaded but withheld by --skip-pon scored
nothing, and recording it as applied would be a lie about the columns present.
validate-outputrequires the PON-derived columns when, and only when, a PON was applied. The contract declared none of them, so the seven WPS columns and FSD'spon_stabilitycould vanish entirely and the gate would pass — which is exactly what happened twice in this release. They cannot be required unconditionally, because--skip-pon, no PON, and a PON the version guard refused are all legitimately unscored;pon_appliedabove is what tells the cases apart.
A directory with no pon_applied at all was written by a build older than
0.9.0, and is reported rather than assumed either way.
-
no_collided_columnsinvalidate-output— a.1-suffixed or repeated column name means a frame collision was written to disk, and no reader can tell which copy is current. Generic and universal: it runs on every table, so the next step that appends instead of replacing is caught without anyone thinking to look for it. It would have found the FSD duplication above unprompted. -
load_pon_modelrefuses a PON older than 0.9.0. Recordingkrewlyzer_versionwas only half a guard:stamp-ponwrites it andvalidate-ponflags it missing, but nothing stopped a run scoring against a model built before the meaning changed — a fabricatedwps_background, six floored σ, a region-MDS fitted over 65–400 bp while samples are measured over 65–1000 bp. Every z-score against such a model is wrong in a way no schema check can see, because the columns are present and finite.
Refusing, not warning. A warning in a log nobody reads is not a guard, and a
plausible wrong number is the failure this release exists to remove. Set
KREWLYZER_ALLOW_OLD_PON=1 to override deliberately — without a documented
way out, someone edits the parquet instead and the version stops meaning
anything.
Writing it exposed two holes it would otherwise not have covered:
motifandregion-mdscalledPonModel.loaddirectly, so those subcommands would have scored against a modelrun-allrefuses — invariant #6. Both now use the guarded loader, and a test forbids any other caller.- The Rust scorers take a parquet path and read it themselves, and that path was set regardless of whether the Python load succeeded. A refused PON still reached FSD, WPS, OCF and region entropy: the Python half stopped and the Rust half carried on.
The floor is a tuple, (0, 9, 0), not a string — it is a compatibility
floor, not the package version, and stays put at krewlyzer 1.5.0.
- The 0.9.0 GRCh37 PONs. All four rebuilt from
run-alloutput via--from-outputs, and the first to passvalidate-ponwith zero findings: no fabricated baseline, no non-positive σ, every required block present, and a cohort digest identifying what each was built from.
test_the_gate_rejects_the_currently_shipped_models has flipped to
test_the_shipped_models_pass_their_own_gate. Its original docstring said
that flip would be the acceptance record for the rebuild; this is it. Two
further assertions join it: every shipped model records its provenance, and
no two share a wps_background — the byte-identical baseline across four
models is the defect that started this release.
plausible_z_scoresinvalidate-output— a z-score above 100 is a near-zero divisor, not biology. Every fabricated-σ defect in 0.9.0 produced plausible numbers and survived every schema check; this catches the opposite failure, where σ is so small the z explodes, which no schema notices either because the column is present, typed and finite.
Runs on every table rather than a listed few, so a new output cannot opt out
by omission, and its columns are read even where the contract does not
declare them. Covers vector z columns element-wise — wps_nuc_z carries 200
values per row, and a per-position σ is exactly where a near-zero divisor
hides. For scale: real WPS z-scores against these models reach 6.7, with
nothing above 10 across 6.9M positions.
-
scripts/check_pon_env.py— behavioural probes that answer "does this install actually have the PON fixes?".--versionreports the previous release until the bump, and a CLI flag only proves the Python side is current; neither can see a stale compiled extension, andgit pulldoes not rebuild one. Each probe drives the real backend on input whose correct answer is known, and fails with the remedy. Run it before any PON build. -
validate-ponchecks a packing list. It read the blocks present in the file and skipped the rest, so a block that vanished entirely was invisible to every check — indistinguishable from one never expected. Demonstrated on the real xs2 model: deletingregion_mds, orfsc_gene_baseline, or almost every baseline at once, produced an identical finding list each time.
That is the hole region_mds fell through when a Rust reader was handed a
format it could not parse. New findings: PON.BLOCK_MISSING for a core block
and PON.BLOCK_MISSING_BAM_ONLY for one that needs BAM input.
fsd_baseline and gc_bias also join the σ checks. Their σ lives in per-bin
columns rather than one *_std, which is why they were skipped — but both
are applied, and a degenerate gc_bias reaches every FSC log-ratio.
-
input_kindin PON provenance —"bam","bed","mixed"or"outputs".mds_baselineandregion_mdsneed a BAM, so their absence is legitimate for a fragment-BED cohort and a defect for a BAM one; without this the gate could only warn either way. Models built before it record nothing and still warn. -
scripts/build_pon_array.sh— one SLURM task per sample, feedingbuild-pon --from-outputs. Extraction across a cohort is embarrassingly parallel and the cluster already has a scheduler for it, so doing it four at a time inside one process leaves most of the machine idle:
| route | 47 samples |
|---|---|
in-process, 4 parallel (build_pon.sh) |
~15.5 h |
in-process, --parallel-samples 0 |
~6.2 h |
| this array, 12 concurrent | ~5.9 h |
| this array, unthrottled | ~1.5 h |
Nextflow would fan out the same way, but every module pins
krewlyzer:0.8.3 and the SLURM profile enables singularity — a run today
would silently execute 0.8.3 and reproduce the defects 0.9.0 removes. That
needs a 0.9.0 container, which needs the release, which needs this PON. An
sbatch array needs neither. --from-outputs will accept Nextflow's output
unchanged once a container exists; it is the same layout.
A finished sample is skipped on resubmission, keyed on the same last-written
marker --from-outputs checks, so a partial directory is never mistaken for
a complete one. SAMPLE_LIST is required rather than derived, as in
build_pon.sh, and a task index past the end of the list is refused —
a short array would otherwise build a quietly smaller cohort.
krewlyzer build-pon --from-outputs DIR— build a PON by aggregating per-samplerun-alloutput directories, reading only files. No BAM is opened.
The in-process builder re-runs feature extraction itself, four samples at a
time on one node: 55–97 minutes per sample, roughly 15 hours for a
47-sample cohort. run-all already computes every one of those features, and
its output directory is a strict superset of what the builder consumes. Every
defect found in the 0.9.0 models so far has been in aggregation, not
extraction — so each one cost a full rebuild to re-check. Re-aggregating now
costs minutes.
Both routes stay and converge on the same collectors, so neither can drift
into producing a model the other cannot read. SAMPLE_LIST and
--from-outputs are mutually exclusive: accepting both would leave it
ambiguous which cohort the digest describes.
A directory is aggregated only if it holds the last file its pipeline
writes (MDS.gene) plus the inputs no baseline can do without. An incomplete
directory is named and stops the build unless --allow-failures is given —
silently aggregating a half-written one makes the cohort smaller than its own
metadata claims.
One difference is unavoidable: EndMotif stores frequencies rounded to six
decimals, while the in-process path keeps them unrounded in memory. At a
typical k-mer frequency of 0.0039 against a cross-donor σ of 2.7e-4, that is
0.37% of σ — the two routes agree on σ to within 0.08% and on z to within
about 0.002. Since Parquet is the product and the product is what rounds,
that is the precision that actually exists.
A 0.8.3 output directory cannot seed the wps_background block. It
predates nrl_at_band_limit, and without that flag there is no way to tell a
repeat length from the edge of the window it was searched in. The refusal
says so. Every other block builds from 0.8.3 output.
krewlyzer stamp-pon— records the release a built PON ships with.
A PON is built from develop, where the version still reads the previous
release, so the model records that however new the code is. Bumping before a
four-hour build would fix it at the cost of putting a release number on
unreleased code; this is the other order — build, then stamp when cutting the
release. Added to the release guide as Phase 2.7, since the sed-based
version update cannot reach inside a Parquet file.
Afterwards krewlyzer_version means the release this model is published
with, not the code that produced it — which is the definition a
compatibility guard needs. build_date is untouched.
It refuses to stamp a model that fails validate-pon. Without that it
would be the shortest path to laundering: run it on one of the models
carrying the fabricated 167.0 / 5.0 baseline and it would claim exactly the
compatibility the guard exists to deny. PON.NO_VERSION alone does not
block, since that is the condition being fixed.
krewlyzer validate-pon— gates the reference, not the results.validate-outputchecks what a run produced; nothing checked the model those results are measured against. Every PON defect fixed in this release sat in four shipped files, visible in the file and invisible to every check.
The load-bearing assertion is the same invariant the output gate enforces one
level up: a baseline that cannot vary with its cohort was not fitted to
one. A single check that σ differs between groups would have caught the
fabricated wps_background in March.
Also checks: σ positive and finite, krewlyzer_version and cohort_digest
recorded, and every entry backed by ≥ 3 samples. Exit codes match
validate-output — 0 satisfied, 1 violation, 2 structural.
All four currently-shipped PONs fail it, which is the acceptance criterion: a gate that passes the models we know are wrong is not a gate. A test asserts that failure, and flipping it is the record that the rebuild worked.
Wired into the Nextflow pipeline as KREWLYZER_VALIDATE_PON, running once
per run rather than once per sample. Wired at the same time it was added —
an unwired module is the defect #34 fixed, and it reads as coverage while
providing none.
- PON provenance:
krewlyzer_version,cohort_digest,cohort_label. The four models in this repository recordn_samplesand nothing else, so none can be reproduced, audited, or compared against a rebuild — and the build script that produced them has a stale comment claiming 47 where the metadata says 21, so even the integer is uncorroborated.
The cohort is a salted, non-reversible digest of the sample IDs, stable
across paths and ordering. Two builds from the same cohort match; the digest
reveals nothing about who is in it. A PON ships in this repository, in the
Docker image and on PyPI, so it is the last place an identifier may appear
(invariant #4). --cohort-label adds a free-text name for humans.
Writing is inert: PonModel.load reads through meta.get(key, default), so
unknown keys are ignored and existing models still load unchanged.
- Per-motif z-scores (
EndMotif.frequency_z,BreakPointMotif.frequency_z).mds_baselinecarries 625 k-mer means and σ — every 4-mer over ACGTN — and nothing read them; the tables shippedMotif, Frequencyalone, so a shift in one motif was invisible unless it moved the whole-sample MDS enough to notice.
The join is on the motif string, never position: 625 baseline keys against 256 ACGT motifs in the output. And the sample is put on the baseline's scale first — sample frequencies sum to 1.0 across the 256, the baseline's expectations for those same 256 sum to 0.972, with the missing 2.79% in the N-containing k-mers the output never reports. Comparing them directly biased every z upward: measured on a real sample, median +0.37 naive against −0.21 corrected.
- WPS z-scoring — the largest PON baseline, previously read by nothing.
wps_baselineis ~128k anchors of 200-element mean and σ vectors, roughly 90% of every PON file, and its only consumer was a log line appending"WPS"to a list of available components.
WPS.parquet now carries wps_nuc_z / wps_tf_z (per-position z vectors)
plus three derived shape quantities with their own z-scores, from a new
wps_shape_baseline block.
Measurement drove every choice here, and ruled out the obvious ones:
| rejected | why, measured |
|---|---|
| mean of z over positions | lag-1 autocorrelation 0.986 — a fragment spans ~167 bp and touches many positions, so such a mean has nothing like σ/√200 precision |
| max of |z| over positions | expected max under pure noise is 2.97, so |z|>2 would flag nearly every anchor |
| centre-minus-flank amplitude | TSS dips at the centre (−6.8 vs −3.4 flanks), CTCF does the reverse — any fixed window is backwards for one |
| raw peak-to-trough range | correlates +0.512 with local_depth; log1p drops it to −0.036 |
| raw shape correlation | bounded at 1.0; mean 0.844 σ 0.099 means 302/400 anchors cannot reach +2 — Fisher arctanh removes the ceiling |
Displacement is measured but deliberately not z-scored.
wps_phase_shift_bp ships because it is cheap and genuinely non-redundant
(corr −0.24 and −0.28 with the two scored statistics), but it gets no
baseline: per-sample mean lag varies by 0.26 bp against a within-sample
spread of 8.43, so there is no whole-sample phasing signal, and per
anchor the intraclass correlation is 0.479 — about half of any lag is
noise, optimistically, since that estimate used a baseline containing the
samples being scored. It is also integer-valued, so on a small cohort its σ
bottoms out at std([0,0,0,0,0,1]) = 0.408 and a 1 bp shift scores z = 2.4.
wps_phase_at_search_limit marks anchors where the ±30 search ended on its
own edge (1.8% measured) — nrl_at_band_limit one level down.
- Per-exon MDS baseline (
region_mds_exon) andMDS.exon.mds_z.MDS.exonis the finest localisation krewlyzer produces — 1,006 rows on xs2, 1,725 on xs1 — and was the only feature table with no baseline at all, so it shipped a raw score with nothing to compare it against.
Keyed on (gene, name): name alone is not unique across genes, and
coordinates would break whenever the panel BED is regenerated. Aggregated in
Python like the FSC gene and region baselines rather than adding a second
Rust entry point for a groupby over ~1.7k rows.
Measured before building it, which refuted the expected shape: exon data is not sparse. Every exon appears in every sample of its assay (7/7 xs1, 19/19 xs2) with a measurable spread, and under 0.25% carry fewer than 10 fragments. So no sparsity handling was needed — only the same NaN-not-floor rule as everywhere else. Verified on the real cohort: 1,006/1,006 xs2 exons fitted, median σ 0.0079, none unmeasurable.
A PON built before 0.9.0 has no such block; mds_z is then NaN and a line
says so.
tests/real_data/— a local cohort gate. PointKREWLYZER_TEST_CORPUSat a scored cohort and it runs; leave it unset and it skips, sopytest tests/stays green and CI never sees it. No cohort data is committed, and no patient identifier may appear in a test name, parameter id, assertion message or artifact —sample_label()gives a non-reversible handle where output needs one.
It exists because every defect the PON audit found needed a cohort to see: a single fixture cannot reveal a constant, cannot match a real PON arm, and cannot show per-anchor sample support. The unit suite proves a scorer can run; only a cohort proves it did.
The four PON blocks that are built and consumed by nothing are marked
xfail(strict=True), so Phase B wiring them turns the test red and forces
the marker off.
-
read_exact_tableincore/output_utils.py— reads the given path with no sibling resolution, the counterpart to the deliberately parquet-firstread_table. -
tests/invariants/test_read_after_write.py— pins both readers' behaviour and an AST check that no function reads with the resolving reader after writing, with a reviewed allowlist that must state why each entry is safe. Also fails on a stale allowlist entry, so it cannot rot. -
tests/invariants/test_pon_format_parity.py— asserts every PON-derived column is present, populated and equal acrosstsv/parquet/both. One test that would have caught FSD, FSC gene, FSC region and the entropy processor together.bothis covered explicitly: it leaves a TSV and a parquet on disk, which is what let the parquet-first reader pick the wrong one. -
FSD was never PON-normalised under
--output-format parquet. Undertsvthe bins came out as log2 ratios; underparquet, raw counts — same column names, no warning either way. 0.9.0 makes parquet the Nextflow default, which would have turned this from affecting nobody into affecting every pipeline run.
Two independent bugs, each sufficient on its own:
- The caller guarded on
outputs.fsd.exists(), whereoutputs.fsdnames the.tsv. The Rust writer honours--output-format, so under parquet no.tsvwas ever written and the whole post-processing block — PON normalisation included — was skipped. Same class asb77909e, missed site; now usesresolve_table_path. _write_fsd_outputread back throughread_table, which is parquet-first: givenx.FSD.tsvit prefersx.FSD.parquet— the raw table the single-pass writer emitted before normalisation. So the log-ratios were computed, logged as41 arms normalized, and then overwritten with the raw counts.
FSD PON: no arms matched is now a warning rather than a debug line: zero
arms means the output is raw counts wearing the same column names as
log-ratios, and nothing downstream can tell the difference.
-
The WPS background PON baseline was never fitted.
pon/build.pylooked for columns namednrl/nucleosome_repeat_lengthandperiodicity/period_score;WPS_backgroundwritesnrl_bpandperiodicity_score. Neither ever matched, so the fallback fired on every build and all four shipped PONs carry an identical hardcoded167.0 / 5.0 / 0.0 / 1.0across all 28 groups — from cohorts of 21 and 47 samples, logging28 groupsas though it had fitted them. Combined with the 0.8.xnrl_bp ≡ 150.0degeneracy, every sample ever produced scorednrl_z = -3.4: a constant, moderately extreme z presented as a measurement. A missing source column is now fatal. -
Six σ floors replaced with NaN.
max(std, 0.001),max(nrl_std, 0.1),.max(0.01)and three siblings do not make a z-score conservative — they make it arbitrarily large, because the divisor is a number nothing measured. That isnrl_at_band_limitone level down: a boundary value reported as a result. NaN propagates to an absent z, which is what "we could not measure this" should look like. A single-sample WPS anchor previously got σ = 0.1 at every position. -
The WPS baseline had no minimum-sample requirement. FSC gene and FSC region have required ≥3 since they were written; WPS — 141k anchors, the largest block by 100× — required nothing, so an anchor seen in one sample still produced a baseline. Measured on the shipped models that was 1.6% of anchors for all_unique/xs1 and 28.8% for duplex/xs1. Now ≥3, and the count skipped is logged.
-
np.std→ddof=1in the PON baselines. These are a sample of healthy donors, not the population. The population form understates the spread and inflates every z built from it — by 2.5% at n=21. -
build-pon --keep-sample-outputs DIR. Per-sample feature outputs were extracted to a temp directory and deleted on completion, so every rebuild re-ran extraction over every BAM from scratch (~4 h for 47 samples) and leave-one-out calibration would have cost n of those. Kept outputs also survive a failed build, so the rerun can skip what finished. -
_log_baseline_qualityreports how many entries in a block carry a measured spread, warns on those that do not, and warns when every entry shares one σ — the signature of a fabricated baseline. -
docs/cli/index.mdcovers every command, and two tests keep it that way: one cross-checks the registered commands against the reference and the README, the other resolves every README documentation link against the pages mkdocs actually publishes. -
krewlyzer report SAMPLE_DIR -o report.html— a single-sample report with a cross-axis verdict, 16 charts and the tables behind them. Internal use: it contains one patient's measurements, so it is generated on demand rather than committed or published.describe-outputremains the shareable, structural view.
Verdict. Four independent axes — fragment size, nuclease cutting, tissue shedding, chromatin accessibility — each showing its value, source column and direction, summarised as "N of M assessable axes agree". Never a composite score: one would hide exactly the disagreement worth seeing. An axis with no PON z-score is not assessable, never counted as disagreement, because that would make a thinner run read as a healthier one.
Organised by table, not by chart. Each section carries one output's chart, its why / how / what, and its columns together. A chart three screens from its data is a chart nobody connects to the data.
PON state is first-class. Without one, every z-score and log-ratio is absent — most of the interpretable surface — so the report leads with a banner and marks each affected column.
Sixteen charts across the five-act structure from the output EDA notebooks, including three taken from the existing notebook's specialised panels: the per-region FSD heatmap, which shows whether the size shift is uniform across chromosome arms or concentrated on a few — something the summed density curve cannot; the WPS foreground profile around TSS and CTCF anchors; and the GC correction curve, since everything GC-corrected downstream is multiplied by it.
Plotly, self-contained (~5.6 MB, no network). The theme follows the docs site,
with an Auto/Light/Dark toggle that the figures respect. plotly is an
optional [report] extra; without it the command still runs and each chart
states that it is missing.
Three rules the charts hold to, each learned from a real defect here: a
constant column is reported as constant rather than drawn as a flat line; no
threshold is drawn as a cut-off, only labelled literature anchors; and a
value that means "nothing was observed" is never placed on the measurement
scale — an mFSD variant with no ALT fragment has ALT_MeanSize = 0, and
plotting 0 − REF would put a fabricated point at the most tumour-like end
of the very axis being read. Those variants appear at the origin with their
own marker, because their absence is itself a finding.
krewlyzer describe-output SAMPLE_DIR— says what a sample's output files contain, which is the question people ask before "is it correct?".
Row and column counts, dtypes, value ranges, distinct counts and null counts
per column, plus total size and which tables are read downstream. Markdown by
default; -o report.html renders a self-contained page.
Everything is either measured from the file or read from contract.py.
Nothing about the biology is restated — that stays in
docs/reference/output-files.md and is linked — so a table that gains a
column or stops being consumed changes the report automatically and the two
cannot disagree.
Each table carries a one-line what it measures and which way cancer
moves it, from a new validate/meaning.py. Direction is the point: it
differs per axis, and getting one backwards is the commonest misreading of
this output — MDS was documented the opposite of its own threshold table
for a year precisely because no test and no report ever stated one. A test
requires a direction for every table except metadata, which is provenance.
No thresholds are recorded. Every numeric band examined turned out to be a display default or refuted outright: the documented ATAC/TFBS entropy range flags a perfectly healthy N(167,30) distribution as abnormal. Directions are robust, magnitudes are cohort-specific, and a number in the registry would acquire an authority it has not earned — a test enforces their absence.
Identifier columns are redacted. Sample directories here are named for
the patient and several tables carry the sample id as a column value, so
the first version of this report leaked a real identifier through Sample
and sample_id examples even though it was generated from deliberately
renamed files — the filenames were clean and the contents were not. Knowing a
column holds an identifier is the useful fact; which one is not, and this
report is meant to be hosted.
scripts/build_gene_bed.py, and gene BED assets rebuilt from GENCODE. The bundled gene BEDs could not answer "which exon is first". The panel assets (xs1/xs2) had five columns and no strand at all; the WGS asset numbered exons in coordinate order, so MTOR — minus strand — carriedexon_num 0at its lowest coordinate while its real exon 1 is at its highest. Anything deriving E1 from those files got the wrong end of every minus-strand gene.
GENCODE's exon_number is transcription-ordered and MANE_Select names the
agreed transcript, so the question is resolved once at build time and written
into the asset as an explicit is_e1 column. Canonical-transcript policy is
--transcript-overrides → MANE Select → GENCODE basic protein-coding →
Ensembl canonical → longest CDS; the build fails if the GTF carries no
MANE tags at all, which is what silently happens with Ensembl's frozen
GRCh37 release 87.
New columns are additive — the first five are unchanged and the existing
parser produces identical genes and coordinates: transcript_id,
exon_number, strand, is_e1, is_alt_e1, is_first_captured.
Three panel symbols (H3F3A, HIST1H3B, PAK7) were renamed by HGNC and
match nothing in a current GTF. They are carried by an explicit alias table
and asserted, because the failure mode is a gene silently losing its
annotation rather than an error.
"First exon" is not one thing. Genes carry a median of 13 distinct annotated first exons, because alternative promoters are the norm, so a single boolean cannot describe what a tile is. Three columns:
| tile overlaps… | column | xs1 / 128 | xs2 / 146 |
|---|---|---|---|
| canonical transcript's exon 1 | is_e1 |
25 | 33 |
| another basic protein-coding transcript's exon 1 | is_alt_e1 |
+15 | +15 |
| most 5′ captured tile (always exists) | is_first_captured |
128 | 146 |
MSK-ACCESS tiles coding hotspot exons, so AKT1's canonical exon 1 sits 15 kb
past the panel's most 5′ tile — but many genes are captured at an
alternative promoter instead, which a MANE-only view misses entirely. 40
xs1 genes have a tile on some basic protein-coding first exon; 88 have
none. is_first_captured exists for every gene but is frequently an internal
exon, which is not a promoter proxy and must not be read as E1.
Which transcript is canonical is configurable per gene via
--transcript-overrides (a gene<TAB>transcript_id TSV), because a panel
designed around specific clinical transcripts should not have MANE imposed on
it. A transcript that is absent from the GTF, or belongs to a different gene,
is a hard error — a silent fall back to MANE would produce an asset that
disagrees with the file the operator wrote.
This also bounds the existing filter_fsc_to_e1 defect: of the 25 xs1
genes with a canonical E1 tile, its lowest-start pick is correct for 18 and
wrong for 7 (6 minus-strand). The larger issue is the 88 with no captured
first exon at all, for which it emits a row labelled E1 that is an arbitrary
internal exon.
krewlyzer validate-output— checks a finished output directory against the contract its consumers rely on. Three layers, in increasing order of what a schema alone can catch: every consumed table present and shaped correctly; domain invariants (frequencies sum to 1, chromosomeschr-prefixed, the six FSC channels partitiontotal, FSD carries only size-bin columns); and anti-degeneracy — a metric that is identical across every sample is an error, not a pass.
That last layer is the point. Run against a real 0.8.3 cohort the schema
passes completely, while WPS_background.nrl_bp ≡ 150.0,
nrl_deviation_bp ≡ 40.0, periodicity_score ≡ 0.3333 and
adjusted_score ≡ 0.0 in every sample: four of the five columns that table
contributes carry no information. A schema-only gate certifies that directory
as good.
Exit codes are 0 satisfied, 1 contract violation, 2 structural (missing
directory, unreadable Parquet) so a workflow can retry on 2 and escalate on 1.
Cross-sample degeneracy below two samples is reported SKIP, never PASS.
Declaring a column legitimately constant requires a written
constant_reason, so silencing a finding costs a justification.
--json-report emits stable finding ids for trend tracking.
scripts/validate_output.py runs it from a checkout.
krewlyzer validate-cohortplus theKREWLYZER_VALIDATE_SAMPLEandKREWLYZER_VALIDATE_COHORTNextflow modules — the gate split into a scatter and a gather so it scales to a real cohort.
A sample directory is ~1.5 GB (WPS.parquet alone is ~120 MB), so re-reading
tens of thousands of them is not viable. validate-output --fingerprint-out
reduces each sample to a ~20 KB fingerprint — a hash and two counts per
column — and validate-cohort compares those. Reading is projected to the
declared columns and bounded by TableRule.scan_rows, which roughly halves
per-sample time and, more importantly, makes memory independent of cohort
size.
The split is not just an optimisation: degeneracy is inherently cross-sample.
On a real cohort every sample passes validate-output individually while
validate-cohort fails, because no single sample can distinguish "this
metric is a constant" from "this is its value here".
run-allnow writes{sample}.validation.jsonand{sample}.fingerprint.jsonon Parquet runs, and--strict-validationmakes a contract violation fail the run.
Emitting is on by default because the fingerprint is a cheap byproduct here —
the tables are already written — and it is what makes validate-cohort
affordable later; leaving it opt-in would in practice disable the only check
that catches a constant metric. Failing is opt-in because a contract rule
that turns out too strict should not take down a cohort. Skipped entirely for
tsv-only runs, since the contract describes the Parquet surface downstream
reads. A checker that throws is caught and logged: it must never lose a
completed run.
WPS_background.nrl_at_band_limit— marks a right-censored NRL estimate, where the spectral peak sat on the edge of the 140-250bp search band sonrl_bpis that bound rather than a measurement.
On real plasma this is 21% of Alu groups for XS1 and 43% for XS2, and
those groups are indistinguishable from interior ones by periodicity_score
or fragment support — so a consumer reading nrl_bp = 250.0 had no way to
tell "the repeat length is 250bp" from "no nucleosomal peak was found". The
same present, plausible, wrong failure mode the NRL fix itself addressed,
one level down.
Not triggered by a long period: 400bp and 2000bp synthetic signals both resolve to 225-235bp inside the band. Only the absence of periodicity pins the edge, so the flag means "no peak found", not "period too long".
Additive and optional, so existing consumers are unaffected; the contract gate accepts directories written before it existed.
- Every feature tool now accepts
--output-formatand--compress. Eight of eleven had neither, sokrewlyzer fsccould only write TSV whilerun-all --output-format parquetwrote Parquet -- and Parquet is all the downstream consumer reads. Anyone running a single tool got output nothing could load, with no error.
The computation was never the problem: standalone wps and run-all produce
byte-identical WPS_background on the same input, because both route through
the same single-pass run_features(). Only the CLI layer hardcoded the
format defaults, and the underlying processors already accepted both
parameters.
scripts/check_output_format.py could not catch this -- it verifies that
internal call sites forward the parameters, not that the CLI ever lets
anyone set them. A new test asserts every feature command exposes both.
The eleven per-tool Nextflow modules now pass params.output_format and
params.compress_tsv through as well; previously only runall did.
Changed
- WPS PON scoring moved from Python to Rust (
rust/src/wps.rs::apply_pon_zscore). It was the last PON z-score still computed in Python, and it was the largest: 89k anchors, a ±30 lag search each, about 5.4M correlations..agents/rules/architecture.mdputs "PON z-score", "loops over >1000 rows" and "row-level computation" on the Rust side; this was all three.
Measured on a real 76,595-anchor output against a shipped PON: 9.5 s against ~17 min, roughly 107×.
Written in Python first on purpose, and that is not an apology — three of its
decisions changed under measurement and would not have survived being written
in the faster language first: the log1p amplitude (raw range correlates
+0.512 with depth), the Fisher transform (bounded r left 302 of 400 anchors
unable to reach +2), and the deliberate absence of a phase-shift baseline
(intraclass correlation 0.479). The port is bug-for-bug, including a tie
in the lag search resolving to the most negative lag — measured incidence 0
in 317 anchors, and correcting it would have destroyed the oracle's ability
to tell a porting slip from an intended change.
The Python reference is frozen as that oracle in
tests/unit/test_rust_python_equivalence.py. Nothing imports it outside
tests; the shipping module is now a ~40-line call. Tolerance was fixed at
1e-6 relative before the first comparison; measured worst case 2.2e-13
on synthetic anchors and 9.0e-14 against real data.
apply_wps_pon now takes the PON path rather than a loaded model, and
baseline_attr is renamed baseline_table — the argument names a table in
the parquet, not an attribute on an object.
- Two silent-read traps found while porting, both of the "present, plausible, returns nothing" class that this release is mostly about:
Every shipped PON stores table and region_id as Arrow large_string,
never string. A reader that downcasts to StringArray gets None and
yields an empty baseline — a legitimate state that degrades silently
instead of raising. pon_model.rs::PonModel::load does exactly that and is
blind to every PON this project has ever shipped. It survived because it has
no callers: each reader loads what it needs, and pub items are exempt
from dead-code warnings. Recorded in architecture.md rather than deleted,
pending an explicit decision on scope.
Vector columns are list<double> from the Python builder and list<float>
from the Rust one. The new reader accepts both and logs the type it could not
read, rather than reporting no anchors.
-
Four
eprintln!("DEBUG: …")calls inwps.rsnow go through the logger. They wrote to stderr on every WPS run, below no level and past every filter. -
The two PON blocks that were built, stored, and read by nothing are now applied. Found by tracing all 21 baseline blocks from build → parquet → read → apply; 19 were reaching scoring code and two were not.
wps_baseline_panel (326 anchors in the xs2 model) — apply_wps_pon ran
only on the genome-wide WPS, so {sample}.WPS.panel.parquet shipped raw.
The anchors closest to the targeted regions were the one WPS output with no
comparison to a healthy cohort. They now carry the same derived and z-scored
columns as the genome-wide file. The shape statistics borrow the genome-wide
wps_shape_baseline: a few hundred panel anchors is too few to fit a second
one, and they overlap the genome-wide set.
gc_bias_ontarget (25 bins) — PonModel.get_mean and get_variance read
self.gc_bias only, so on-target FSC and FSR normalised against the
genome-wide GC curves. Capture enrichment shifts the GC profile, which is the
entire reason panel mode fits a second block. {sample}.FSC.ontarget
*_log2 and *_reliability values change as a result. A PON without the
block falls back to the genome-wide curves rather than dropping the column.
- A scored column is NaN when its baseline is absent, not 0.0. Six
fabrications in the apply path, after
zscore_or_nanremoved nine from the baseline classes:
| where | when | was | now |
|---|---|---|---|
fsc_processor |
no gc_bias |
{ch}_log2 = 0.0 |
NaN |
fsc_processor |
no variance | {ch}_reliability = 1.0 |
NaN |
fsr_processor |
no long fragments | ratio = short_count |
NaN |
fsr_processor |
ratio ≤ 0 | short_long_log2 = 0.0 |
NaN |
region_entropy.rs |
σ unmeasurable | z = 0.0 |
NaN |
region_entropy.rs |
σ column absent | σ ← 1.0, mean ← 0.0 | NaN |
Zero is never the cautious choice. A log2 ratio or z-score of zero says
"this sample sits exactly at the healthy baseline" — the most confident
statement either column can make, asserted precisely when nothing could be
compared. Measured: with gc_bias absent, three windows whose raw values
differed fourfold all read 0.0.
The Rust pair defeated the Python fix until both were done. entropy_std
defaulted to 1.0 and entropy_mean to 0.0 — a standard normal — so a
builder storing NaN for a single-donor label produced the raw difference
rather than no score. And position(...).unwrap_or(0) read column zero
when a column was not found, which in a PON parquet is table.
region_entropy_processor was the fifth site of the sample_std_or_nan
defect and the only one outside pon/build.py.
scripts/build_pon_unfiltered.sh→scripts/build_pon.sh, taking assay and variant as arguments. The old script was hardcoded to xs2 and carried a header claiming "47 samples" inherited from the xs1 copy it was made from — while the xs2 model it produced records 21. Four models built from four edited copies of one script is how that happens.
It also passes --keep-sample-outputs and --cohort-label, and ends by
running validate-pon on what it built. A build that produces a model the
gate rejects has not succeeded, whatever build-pon's exit code said.
docs/guides/building-pon.md gains the rebuild runbook: what to watch in the
log, why the sample outputs are kept, and the LFS ordering that has to be
right.
-
MIN_SAMPLES_PER_KEYnames the ≥3 floor once inpon/build.py. FSC gene and FSC region have required it since they were written and WPS acquired it in4cd634b; the four now agree by construction rather than by three separate literals happening to match. -
OCF PON logging distinguishes "absent from the baseline" from "matched but unscoreable". Both correctly write NaN, but collapsing them reported
0/1 regions matchedfor a region that matched perfectly well and simply has no variance in the cohort. That is the difference between "rebuild the PON" and "expected", and only the log can say which. -
tests/integration/test_pon.py's PON fixture used a schema that has never existed —version/genome/sample_countwhere real PONs writeschema_version/assay/n_samples. SincePonModel.loadreads throughmeta.get(key, default), it loaded as a completely empty model and the tests asserted onlyis not None.test_pon_model_loadinghad been markedskip("PON parquet schema requires production format")— true of the fixture, not the loader. Fixture corrected against the bundled PON, test unskipped, and both now assert real field values. -
PON scoring was skipped under
--output-format parquetin three more places. Same shape as the FSD defect below: callers name a hardcoded.tsv, the Rust writer honours--output-format, and a bare.exists()is then False. -
apply_fsc_gene_ponreturnedNonewith "file not found" and added nodepth_zscore apply_fsc_region_ponand the e1-only filter, identicallyprocess_region_entropyskipped the whole step
All four now resolve with resolve_table_path. Measured before and after on
the bundled PON: depth_zscore absent → present with 5/5 populated.
- TFBS and ATAC emitted
z_score = 0.0when there was no baseline. Zero is not a neutral placeholder — it is the most confident claim the column can make ("this sample sits exactly at the healthy baseline"), asserted on the strength of having no baseline at all, and indistinguishable from a measured zero. NowNaN, with a WARNING naming how many labels are affected. Applies both when no--pon-modelwas given and when the PON lacks the table.
Removed
- 271 lines of dead Rust.
_core.wps.apply_pon_zscore(173 lines) plusload_wps_baseline_from_parquet,phase_shift,PHASE_MAX_LAGand the parquet imports that existed only to serve them. The crate now builds with zero warnings.
_core.wps.apply_pon_zscore in detail: It emitted one scalar z
per anchor from wps_long_std/wps_short_std — v1.0 baseline field names,
while every shipped PON is v2.0 vector format — and pushed 0.0 where σ was
not positive. It was also unreachable: the call was gated on
pon._source_path, an attribute nothing in the codebase sets. Dead,
schema-obsolete, fabricating, and implementing the statistic measurement
refuted.
PonModel.save— a second serializer that wrote only the metadata block, producing a PON with no baselines at all, whilebuild-ponused_save_pon_model. Nothing in production called it. Removed rather than completed: a second writer is a second thing to keep in step with every new block.
Fixed
pon_stabilitywas an unnormalised inverse variance and collapsed to zero.1 / (mean(σ²) + 0.01)only ever looked sane because it was fed the wrong σ —get_statsreturned the first size bin's spread for every size. Once each bin got its own, the true magnitudes showed: FSD σ is in fragment counts, 111–4,573 on xs1 and 947–35,961 on xs2, somean(σ²)reached 1.6×10⁶ and 8.3×10⁷ and the reciprocal wrote0.000000at six decimals for every arm. The0.01floor assumed σ ≈ O(1); it never was. Caught by the anti-degeneracy check on a real XS2 sample — invariant #1 working.
Now 1 / (1 + mean(CV²)) with CV = σ/expected, bounded in (0, 1].
Dividing by expected is justified, not convenient: measured on both models,
d log σ / d log expected is 0.83 and 0.95 — close to 1, so the
spread is between-donor variation rather than Poisson counting noise (which
would give 0.5), and the ratio is genuinely scale-free.
On real output: 41 distinct values per assay at six decimals, against 5 and 1 before. xs1 spans 0.714–0.758, xs2 0.865–0.926 — and that difference is interpretable, the 21-donor xs2 cohort agreeing more tightly than the 47-donor xs1 one. The between-arm spread is only ~6% but it is real rather than estimation noise: split-half over odd/even bins gives r = 0.834 and 0.997, a full-length reliability of 0.91 and 0.999.
Old and new values are not comparable — the scale changed from an unbounded inverse variance to a bounded agreement score.
-
pon_stabilitywas declaredvary=BOTH, which it cannot satisfy. It is computed from the PON's σ alone, with no sample value entering it, so every sample scored against one model gets the same 41 numbers. Nowvary=WITHIN, which still requires it to differ down the rows — the check that caught the collapse above. -
The gate could not see four of the six tables Parquet mode dropped.
.OCF.parquetand all three.synctables were absent from the output contract; only the two on/off-target summaries were declared, and they are what reported the loss. All four are declared now, withocf_zmarkedrequires_pon— the.synctables carry the per-position profiles the summaries are reduced from, which is the large half of OCF. -
FSC.regions.is_e1/is_alt_e1were declaredstringand writtenint64.gene_bed.pyparsesfields[8] == "1"into a Python bool, which lands as int64 0/1 — so the contract was wrong, not the writer. The synthetic cohort wrote strings too, so fixture and contract agreed with each other and neither with the code, which is how it survived. These columns are new in 0.9.0, so nothing older disagrees. -
fsc_gene_ratios_sum_to_oneflagged regions that have no fragments. All six channels are zero, so the sum is zero and the row reads as a maximal partition failure — "worst deviation 1.000000" — when nothing is wrong with it. 3 of 1725 regions on a real XS1 plasma sample, 13 of 1725 on a shallow one, making the count a function of sequencing depth and burying any real break among them. Zero-fragment regions are now excluded; a region with even one counted fragment must still partition it. -
BreakPointMotifwas scored against theEndMotifbaseline — and so were the two on-target motif tables.apply_motif_ponusedpon.mds_baselinefor all four, an end-motif, genome-wide block.
An end motif is the 4-mer at the fragment's 5′ terminus, what the nuclease
left. A breakpoint motif spans the cut site and includes reference bases
not present in the fragment. The two frequency vectors correlate only
0.696 on the same sample, so frequency_z there measured the offset
between two definitions rather than any departure from healthy:
EndMotif median / >10 |
BreakPointMotif median / >10 |
|
|---|---|---|
| XS1 | 1.82 / 0 of 256 | 5.85 / 70 of 256 |
| XS2 | 4.47 / 39 of 256 | 11.25 / 136 of 256 |
A correctly fitted baseline gives a median near 0.67 — the half-normal median. Demonstrated on synthetic cohorts with the real correlation structure: median |z| 6.75 with 86 of 256 beyond 10 against the end-motif baseline, 0.75 with none beyond 10 against its own.
The PON gains breakpoint_motif_baseline and
breakpoint_motif_baseline_ontarget, built by both routes. mds_mean and
mds_std are NaN in those blocks: MDS is defined on end motifs, and a number
there would be a different statistic under the same name.
mds_baseline_ontarget needed no new work — it already existed and was
already used for the whole-sample MDS z, just not for the per-motif
frequencies. A PON without the breakpoint blocks (every one built before
0.9.0) now yields no frequency_z on those tables rather than a wrong one.
- A σ of 10⁻¹⁷ was used as a divisor. Every PON built to date carries
positions whose baseline spread is floating-point residue: 4.6 % of positive
σ in xs1.all_unique, 12.0 % in xs2.all_unique, 47.2 % and 55.4 % in the
duplex pair. On a real XS2 plasma sample that produced 728,007
wps_nuc_zabove 100 and 354,260 above 10⁶, peaking at 6.1 × 10¹⁸.
WPS is (fragments spanning a position) − (fragments ending near it), and each
fragment carries a fractional GC-correction weight. Where no donor had
coverage, or the two terms cancel, the true value is zero — but it computes
as ~10⁻¹⁷ rather than 0.0, so the exact min == max test added earlier in
this release never fired. wps_tf corroborates it: 0 % residue and 185,079
honest NaN, because TF-sized fragments are rare enough that those positions
sum to exact zero.
SIGMA_FLOOR = 1e-9 now applies in the builder (element_wise_std,
mean_and_sd) and in the scorers, so the four models already in the wild
are safe without a rebuild. Not a taste threshold: across the 24.7 M positive
σ of the shipped xs1.all_unique baseline, 1,177,647 sit below 10⁻¹², none
sit between 10⁻¹² and 10⁻⁶, and the rest are measurements — 10⁻⁹ is the
log-space midpoint of six empty decades. Mirrored in
pon/model.py::SIGMA_FLOOR and asserted equal in validate/claims.py.
Measured on a real XS1 plasma sample against the unchanged shipped PON: max |z| falls from 2.8 × 10¹⁷ to 3.7 × 10⁴ and values above 10⁶ from 18,398 to zero.
scripts/check_pon_env.py gains a matching probe. The existing one asserts
that identical donor vectors give NaN and passed throughout, because the
real case is near-identical; the new one was confirmed to fail against the
pre-fix binary before it passed against the fixed one.
--output-format parquetsilently discarded every OCF table. All six —OCF,OCF.sync,OCF.{on,off}target,OCF.{on,off}target.sync. Measured on a real XS1 plasma BAM: 44 logical tables underboth, 38 underparquet, exactly the OCF six missing.
The OCF temp directory is an internal intermediate — Python moves the files
to their final names and converts them, because _core.ocf.apply_pon_zscore
reads TSV line by line. pipeline.rs nonetheless dispatched that
intermediate write on the user's --output-format, so Parquet mode wrote
all.ocf.parquet while the mover looked for all.ocf.tsv. Nothing matched,
nothing moved, the temp directory was deleted. No warning, because an absent
file is indistinguishable from "OCF was not requested".
The FSC counts three lines above already carried the right rule — "an
internal intermediate; keep as TSV". OCF follows it now, and the mover
reports an empty intermediate as an error. The .sync tables are the large
ones (483 KB TSV vs 56 KB Parquet), so this removed orientation-aware
fragmentation from precisely the format the release ships in (invariant #2).
-
On-target FSD was never PON-scored under Parquet. The same defect one branch below its own fix:
.exists()on a hardcoded.FSD.ontarget.tsv, fifteen lines under the comment explaining why the genome-wide branch had stopped doing that. Confirmed on real XS1 and XS2 plasma samples — 69 columns and zero_logRagainst the genome-wide table's 137. Resolved withresolve_table_path; an absent table is now warned about. -
On/off-target OCF left its intermediate
.tsvbeside the Parquet.apply_ocf_python_ponwrote the final format but never calledcleanup_intermediate_tsv, which the genome-wide path does — so a Parquet run shipped two files for one table with no way to tell which was current. -
fsd_only_size_binsfired on every PON-scored FSD table. The check's reserved column set was written for the raw table and never updated when PON scoring began writing into the same file, so it reported all 67{bin}_logRcolumns pluspon_stabilityas stray — 68 spurious errors on every scored sample. A gate that fires on correct output is a gate nobody reads. It went unnoticed because the synthetic cohort carried no PON columns until now, so the check had never once seen the shape it runs against in production. -
The WPS column reference documented five columns that do not exist, missed eight that do, and typed the profiles as scalars. Found by diffing the reference against a real scored table instead of reading it.
wps_nuc_smooth, wps_tf_smooth, wps_nuc_mean, wps_tf_mean and
wps_tf_z were never written. capture_mask, local_depth and the six
derived PON columns were. And wps_nuc, wps_tf, prot_frac_nuc,
prot_frac_tf and wps_nuc_z are 200-element lists, documented as
float — so the three worked examples indexed columns that do not exist and
raised KeyError. All three are rewritten and were run against a real
76,595-anchor table before being committed.
A sweep now checks every df["col"] reference in every Python block of
output-files.md against real schemas from a full run: 0 of 40 table
sections reference a column that does not exist.
-
FSD's histogram range was documented as the filter range.
meaning.pysaid "5 bp bins over 65–1000 bp"; the histogram is 67 bins over[65, 400). Fragments of 400 bp and above are excluded from the bins and fromtotal, so a long-fragment fraction computed againsttotalhad a denominator that silently excluded them. Corrected, with the range pinned invalidate/claims.pyso the two cannot drift again (invariant #5). -
The FSD PON columns were documented nowhere.
{bin}_logRandpon_stability— the entire product of PON scoring for FSD — appear in neitheroutput-files.mdnormeaning.py. Both are now described, including whatpon_stabilitymeans and why it is NaN rather than 0.990 when no σ could be measured. -
WPS z-scores were written to a file nothing reads, then destroyed their own input when that was fixed. Two defects stacked, both introduced by the Rust port earlier in this release and both caught by reviewing the product rather than the code.
run-all passes {sample}.WPS as the output base — a stem carrying a
compound extension. Path("s.WPS").with_suffix(".parquet") reads .WPS as
the suffix and replaces it, so the scored 18-column table landed on
{sample}.parquet while {sample}.WPS.parquet stayed the raw 11-column
profile. Downstream reads Parquet only (invariant #2) — the same failure the
--output-format fix removed earlier in this release, in a second disguise,
under a comment that asserted the correct behaviour while the code did the
opposite.
Fixing the path exposed the second: the normal case is in place, and the
scorer streams its input in batches while File::create truncated that same
file out from under the reader. It now writes a sibling temp and renames
atomically, so a crash leaves the original intact rather than a half-written
product, and a return of 0 still writes nothing.
Every earlier test used a simple stem (tmp_path / "o"), which has no
compound extension and so could not see either. The new tests use the shape
the caller actually passes.
-
The four shipped PONs are stamped
0.9.0. They were built fromdevelop, where the version still read0.8.3, so every model recorded the previous release however new the code was.krewlyzer_versionnow means the release a model is published with — which is what the compatibility guard needs — whilebuild_datestill records when it was built. Only the metadata row changed: all four cohort digests, block counts and row counts are unchanged, andvalidate-ponpasses on the stamped files, which are the files that ship. -
validate-ponreported0 sample(s)after checking four models. The summary line is shared withvalidate-output, which counts samples; the PON gate fills no samples, so it printed a line saying nothing was inspected directly above one saying the contract was satisfied. It now reports what it actually checked:4 model(s). -
Every FSD log-ratio was scored against the wrong size bin. The largest output defect found in this release, and it reached every sample ever normalised against a PON.
size_bin is an integer in every sense that matters, and it is stored as a
double: a PON is one long-format table, so every row belonging to another
block carries a null there and the union column comes out float64. The
reader called row.get_int(), which errors on a Double rather than
coercing, and the error fell through to unwrap_or(0).
So every size_bin became 0. An arm's size_bins was 67 zeros,
size >= *size_bins.last() held for every size, and get_expected returned
the last row's expectation for all 67 bins.
Measured on a shipped PON and a real sample: 41/41 arms matched the last-bin baseline exactly and the bin-matched baseline not at all. After the fix, 41/41 match bin-matched and 0/41 match the old behaviour. The correction is large — median log2 shift −1.05, maximum |Δ| 5.10.
Nothing could have caught it downstream. The log-ratios still varied across bins, because the sample numerator varies, so they were present, finite, non-degenerate and wrong — invariant #1's failure mode one level subtler than the constant it was written for.
The PONs are not affected and do not need rebuilding; the reader was broken, not the model. Samples normalised against a PON do need re-running.
-
pon_stabilityused one size bin's sigma for the whole arm.FsdBaseline::get_statsreturnedstd.first()regardless of the size asked for. Measured on the shipped xs1 all-unique PON: sigma varies 41.6× (median) across an arm's 67 bins, up to 56.4×, sopon_stabilitywas wrong by a median of 4709% on all 41 arms. Sigma is now interpolated at the requested size, exactly asexpectedalready was. -
The same column-lookup defect closed in
ocf.rs, the third and last copy.position(...).unwrap_or(0)made an absent column read column zero, andocf_mean/ocf_stddefaulted to0.0/1.0— together making z the raw difference from zero. OCF's values parse correctly today (verified against a shipped PON: an observation 2σ above the mean scores exactly 2.0), so this closes a landmine rather than fixing a wrong number. It would have fired the moment a column was renamed. Thetablelookups infsd.rsandregion_entropy.rsare by name now too.
Two unwrap_or(0) lookups remain, both on TSV headers rather than PON
columns, where column zero genuinely is the label in most files. Left alone
deliberately.
-
Three more fabricated defaults on the FSD baseline path, the same family removed from nine baseline classes as
zscore_or_nanthis release: an unreadable sigma defaulted to1.0(which with the 0.01 floor reads as apon_stabilityof 0.990, indistinguishable from a measurement), an unreadable expectation to0.0, andposition(...).unwrap_or(0)made a column that could not be found read column zero —table— the same defect already removed fromregion_entropy.rsand missed here. Rows with an unreadablesize_binorarmare now skipped and counted in a warning rather than silently collapsing onto bin 0 or an empty key. -
WPS z-scores were written to a file nobody reads.
apply_wps_ponhonoured--output-format, whose default istsv, while the Rust step writes.WPS.parquetunconditionally. A default run therefore produced two files: a.WPS.tsvcarrying all seven PON columns and a.WPS.parquetcarrying none. Downstream reads Parquet only (invariant #2), so WPS scoring was present on disk and absent from the product.
Measured on a real 89,034-anchor run: the TSV had 18 columns at 928 MB, the parquet 11 at 144 MB, and the parquet was the older file. WPS is Parquet-only by contract precisely because of that size ratio — 200-element vectors are ~6× larger as text.
The writer no longer takes an output_format at all, so the knob cannot be
set wrong; when a run requests another format, the WPS exception is logged
rather than silently applied. Introduced in 8265ad3, and missed because
every unit test passed output_format="parquet" explicitly.
- FSD log-ratios accumulated on every re-run.
apply_pon_logratiofound its size bins by taking any header whose text before the first-parsed as an integer — so65-69_logR, a column it had written itself, came back as bin 65 on the next pass. Each run then appended a fresh set of log-ratios and log-ratios of the previous ones.
Measured on a real 67-bin sample: 69 columns raw, 137 after one pass, 273
after two, every _logR name repeated. _write_fsd_output reads that back
with pandas, which renames the collisions to _logR.1 / _logR.2 and writes
those to disk. The correct answer was in the file three times with no way to
say which was current — and read_csv returns the first, which is the oldest.
Nothing marked the work as done, so a rerun into the same directory, a
Nextflow retry or a resumed job all triggered it. The normaliser now drops the
columns it previously produced and rebuilds them, so running it twice gives
byte-identical output; a size bin is recognised as exactly {int}-{int}; and
it logs when it is replacing rather than appending.
The PONs are unaffected — verified, not assumed. Their cohorts were built
with --skip-pon, so no _logR column ever existed for the builder to
mistake for a bin, and all four models carry exactly 41 arms × 67 bins with no
duplicates.
-
Two more fabrications in the same function.
log_ratio = 0.0when the arm is absent from the baseline, andpon_stability = 1.0when no σ was measured. Zero says "exactly at the healthy baseline" and 1.0 is what an average variance of 0.99 would give — both specific readings asserted where nothing was measured. Now NaN, as inzscore_or_nan, the FSC log2, the FSR ratio and the entropy z. -
The vector σ had the same residue defect, in a different function. The previous fix went to
mean_and_sd, which serves the scalar shape statistics.wps_baseline's 200-elementwps_nuc_stdcomes fromelement_wise_std, a separate implementation carrying its ownsd > 0.0guard — so the fix left every vector σ untouched, and a re-aggregated cohort came back byte-identical.
This is the one that matters. Measured in the shipped models: residue σ
reaches 4.5% of usable positions in xs1.all_unique, 11.8% in
xs2.all_unique, 47.2% in xs1.duplex and 55.4% in xs2.duplex, touching
38–74% of anchors. A typical real σ there is 0.4–1.2, so a one-unit deviation
scores z ≈ 10¹¹.
- The Rust
mean_and_sdreported floating-point residue as a spread.if sd > 0.0 { sd } else { NAN }is the same holesample_std_or_nanclosed on the Python side, and it misses the case it exists for: summation error grows with n, so identical donors give residue rather than a clean zero. 21 copies of0.1— the xs2 cohort size — yield7.6e-9; 8 copies of-0.3333333yield3.2e-8. Both pass> 0.0and become divisors. The existing test used[3.0, 3.0, 3.0], which cancels exactly, so it passed throughout.
Measured in the freshly rebuilt models, where these reach compute_z_vector
through the same std > 0 guard:
| model | positions affected | anchors touched |
|---|---|---|
xs1.all_unique |
4.5% | 38% |
xs2.all_unique |
12.2% | 57% |
xs2.duplex |
55.4% | 74% |
A typical real σ there is 0.4–1.2, so a one-unit deviation against a residue σ scores z ≈ 10¹¹. Those are positions where every donor looked the same, which is exactly where a sample doing something unusual should have been reported as absent rather than astronomically significant.
Identity is now tested on the values (min == max), which is exact and
invents no tolerance — values one ULP apart keep their spread, pinned by a
test.
Found by inspecting wps_baseline, the block that is 99% of the file and the
one validate-pon cannot check, because its σ lives inside 200-element
vectors rather than a column.
build-ponran region-MDS at a different fragment cap thanrun-all. It hardcodedmax_len=400;run-alluses 1000, andrun-allis how every sample is measured at scoring time. MDS is normalised entropy, so the wider window sees more fragments and reads higher — a baseline fitted over 65–400 bp and a sample measured over 65–1000 bp are not the same quantity.
Measured on the xs2 duplex cohort, the same 21 donors aggregated both ways: all 146 genes shifted up by 0.0043 ± 0.0016. Against the baseline's own σ of 0.0043 that is a median z bias of +1.15 — 86 genes past one σ and 11 past two — in every healthy sample, from the cap alone.
Both call sites now use sample_params.maxlen, which defaults to the same
1000. Invariant #6, caught by diffing a rebuilt PON against its predecessor
rather than by any test, because nothing crashed.
The four models rebuilt via --from-outputs are unaffected: they aggregate
run-all output, so they already carry the 1000 bp measurement.
- A kept
build-poncache could not be re-aggregated.--keep-sample-outputsexists so a rebuild does not re-read every BAM — 55–97 minutes per sample, ~15 hours for 47 — and every defect found in the 0.9.0 models has been in aggregation, so re-checking one should cost minutes.
It could not, because process_sample never called write_motif_outputs:
run-all did it at its own call site and build-pon did not. A kept cache
therefore had no EndMotif or MDS table, even though the same k-mer counts
reached the model through memory. That is the one input --from-outputs
cannot reconstruct, and krewlyzer motif needs a BAM, so it could not be
backfilled either. Measured against the real 0.9.0 cluster cache: all 47 xs1
duplex directories were refused, each with missing .EndMotif.
process_sample now takes write_motif_files, which build-pon passes.
It defaults to False so run-all keeps writing them itself and nothing is
written twice. The call has to live there because write_motif_outputs takes
an ExtractionResult, which is live inside process_sample and never
returned.
The resume story is deliberately "re-aggregate the cache" rather than
"skip finished samples in the extraction loop". The in-process path
collects from live SampleOutputs objects, not from the directory, so
skipping extraction would skip collection too and drop the sample from every
baseline without saying so. On success the build now prints the exact
--from-outputs command instead of an unactionable "reuse them".
-
--from-outputsrequired extraction assets it never reads. No BAM, no reference, no bin file — butbuild-ponrefused when the bundled bin file was missing, before reaching the branch that would not have opened it. The bin file is LFS-backed, so every--from-outputsbuild failed in CI while passing locally. -
Two Rust baseline readers parse plain TSV and were being handed parquet.
compute_fsd_baselineandcompute_region_mds_baselineread withBufReader::lines()— no gzip, no parquet.File::opensucceeds on both anyway, the header parse then finds no usable columns, and every sample is skipped: three samples in, zero out.
Not an edge case. A run-all directory writes .parquet and .tsv.gz and
no plain .tsv, so it is the normal case for the layout --from-outputs
reads. Found against the real 0.8.3 corpus, where region_mds came back
empty from the same files region_mds_exon read fine.
The two failed differently and the quieter one was worse. FSD raised, but
blamed the backend — "no data returned from Rust", the same misleading
wording already fixed once in _compute_wps_baseline. region_mds logged
one warning and returned None, so the block was simply absent from the
model; validate-pon iterates the blocks present in the file and skips empty
ones, so a block that vanished entirely cannot be checked. Build, gate and
downstream all reported success.
Inputs are now normalised to plain TSV where the constraint lives, in the two callers that know about it. Unreadable inputs are reported rather than dropped. Both empty-result paths raise and name an input file.
-
region_mdsfell back to a standard normal.data.get("mds_mean", 0.0), data.get("mds_std", 1.0)makes z equal the raw MDS — about 0.95, an ordinary-looking number for a statistic that was never computed. Third instance after_compute_mds_baselineandget_periodicity_stats. -
A zero spread was reported as a measurement in three of the four PON builders.
pandas.std(ddof=1)returns0.0for identical values, not NaN — NaN only whencount < 2. Only_compute_fsc_gene_baselineconverted that; the other three did not, while carrying comments claiming they did ("NaN where pandas could not measure a spread", "an unmeasurable spread yields NaN"). A comment asserting behaviour the code lacks is worse than no comment: it stops the next reader checking. The Rust twinmean_and_sdwas correct, which is why the WPS blocks passedvalidate-ponand the Python ones did not.
One sample_std_or_nan helper now serves all six sites. Identity is tested
on the values rather than the result: std([0.95, 0.95, 0.95]) is
1.36e-16, not 0.0 — cancellation in the variance sum leaves residue, so
a sd <= 0 guard misses the exact case it exists for, and a z divided by
1.36e-16 is 1e16, worse than the zero it replaced.
- The PON averaged the edge of the NRL search band in with real
measurements.
nrl_bp = 250is the top of the FFT window, andnrl_at_band_limitsays so on every affected row — the builder ignored the flag. Across the xs2 duplex cohort all 174 band-limited rows carry exactly 250.0, one unique value with zero variance, against 194.9 ± 24.0 for the other 414. Twelve of 28 groups are partly limited, four entirely, so those groups reported the edge of the search as the healthy expectation.
The NRL is now fitted from the rows that measured one, and the group keeps
its row either way — absent, not dropped, because a group xs1 can measure
and xs2 cannot is information when the two models are compared. New
n_at_band_limit and n_nrl_fitted columns record why a baseline is
missing. Below MIN_SAMPLES_PER_KEY measured rows the mean goes too, not
just the spread: a cohort baseline averaged over one donor is the same
fabrication as a hardcoded one.
Periodicity is deliberately still fitted from all rows. Those same 174
band-limited rows hold 174 distinct periodicity values spanning 0.37–0.86 —
only the peak's position hit the edge, its strength was measured. nrl_bp,
periodicity_score and now nrl_at_band_limit are all required.
WpsBackgroundBaseline.get_periodicity_statsinvented a standard normal.row.get("periodicity_mean", 0), row.get("periodicity_std", 1)made a missing column produce a z equal to the raw periodicity — about 0.47, an unremarkable number nobody would have questioned. The same fabricated baseline the block itself was rebuilt to remove, one level down in the reader. Now indexed, so a missing column raises.
Verified against all three models already built on the cluster by rebuilding
only the affected blocks from cached per-sample outputs — no BAM re-read.
xs1.duplex went from three PON.NONPOSITIVE_SIGMA findings to none;
xs2.duplex and xs2.all_unique from one each to none. All three now exit
0 under validate-pon.
- A
region-mdsfailure during a PON build was swallowed byexcept: pass.region_mdsandregion_mds_exonare both built by globbing the files that step writes, so a failure means two baselines silently missing from the model — and the build reported success. Found while inspecting a live rebuild: 4 of 4 completed samples had noMDS.gene.tsv, andgrep -i mdsover 1,803 log lines returned nothing, because there was nothing to find.
Now warns with the exception type and message, and separately says when the gate did not open at all — no gene BED resolved, or the input is a fragment BED rather than a BAM. At aggregation, a BAM cohort that produced no MDS files at all is an error rather than an empty block.
scripts/build_pon.shrefuses to run on a pre-0.9.0 krewlyzer. A 0.8.xbuild-ponrebuilds every defect this release fixes, exits 0, and produces logs that look clean — the failure is invisible until someone opens the model. Caught in practice: the first 0.9.0 rebuild attempt ran 18 minutes onv0.8.3before the banner was noticed.
Also corrects the sample-list default from {assay}_all_unique_pon.txt to
{assay}_allUniq_pon.txt, which is what the lists are actually called — the
script had invented a convention rather than using the existing one.
nrl_zandperiodicity_znever reached a single output file — and the machinery was fully plumbed.compute_nrl_zscoredefaulted togroup_id="all", and no PON has ever held a group by that name: they areGlobal_All,Chr1_All…Family_AluY. Every lookup missed, the function returnedNone, andNonewas indistinguishable from "this PON has no baseline". The default is now a namedGENOME_WIDE_GROUPconstant, defined once, and a group that matches nothing logs what it did find.
Two defects had to stack for the column to be absent: the fabricated
wps_background baseline (4cd634b) and this. With the old PON the value
is −3.4 for every sample, which is what the fabricated 167.0 / 5.0 gives
for an nrl_bp pinned at 150.
-
27 of 28 NRL baselines were built and never used. Only
Global_Allwas scored, though the output carries annrl_bpper chromosome and per Alu family and the baseline has one for each. All 28 groups are now scored; per-chromosome NRL drift is the reason those baselines exist. -
An unmeasurable spread produced
z = 0.0, in twelve places. Ten baseline classes inpon/model.pyindependently ended withif std > 0: return (x - mean) / stdfollowed byreturn 0.0— nine copies of the same three lines, and the same defect nine times. Zero is not a cautious answer: it is the most confident claim the column can make ("sits exactly at the healthy baseline"), asserted precisely when the baseline measured no spread, and indistinguishable from a genuine zero.
Replaced by one zscore_or_nan helper, so it can only be wrong once. Also:
FsdBaseline.get_expected/get_std returned 0.0 for an unknown arm (a
caller dividing by that gets infinity, not an absence), and
compute_shape_score returned 0.0 for an undefined correlation, which is
a real claim about shape agreement.
-
WpsBaseline.compute_z_vectorsubstituted1.0for an unusable σ. Since4cd634bthe builder writes NaN where it could not measure spread;np.where(std > 0, std, 1.0)is False for NaN, so every one of those became a finite z — undoing the builder's honesty at the read side, which is the harder place to notice. This is the function WPS z-scoring will use. -
Reading back a file just written could return an older run's data.
read_tableis parquet-first: givenx.tsvit prefersx.parquet. That is right when asking "what was produced for this output" and wrong immediately after writing a specific file. The Rust backends write plain TSV and Python reads it back to honour--output-format, so a.parquetsibling left by an earlier run into the same directory won — which is what a Nextflow retry or-resumeproduces.
Reproduced end to end for region entropy: with a stale parquet present, the
output carried the previous run's z_score and entropy rather than
this run's. Fixed at every read-after-write site — region entropy, mFSD, UXM,
region-MDS gene and exon, on-target MDS — via a new
read_exact_table helper.
Also reverts an over-correction from earlier in this branch: the post-Rust
existence check in process_region_entropy must be exact, not resolving, or
a stale sibling satisfies it and the degradation branch never fires.
-
A whole-genome run lost an entire verdict axis. A panel run splits OCF into on- and off-target; a WGS run writes neither and emits a plain
.OCF.parquet.verdict.py,plots.pyand the report's run-facts strip all read only the panel pair, so every WGS sample reported OCF absent with the file sitting in the directory. Losing an axis makes the verdict read stronger, not weaker — it shrinks the denominator — which is the dangerous direction. The preference order now lives once, incontract.py, and all three read it. -
Six optional outputs rendered with no interpretation.
mFSD,UXM,OCF,FSC.regions.e1only,fsc_countsandcorrection_factorsgot a section with a measured shape and nothing else — no lede, no direction, no why/how/what — becausetest_meaning_registry.pykeyed onCONTRACTalone and those six areNOT_CONSUMED. mFSD in particular had a chart with no explanation beside it. The registry now covers everything either command can render. -
Numeric identifier columns were not redacted.
describe-outputprefers a min…max range over the example value whenever one exists, so an integer patient key or accession was printed in full by the branch meant to hide it — a range over one distinct value is the value. -
describe-outputandreportdisagreed about the same fact. The report says gated / inventory;describe-outputstill said read downstream, a claim about another repository that nothing here keeps in sync. Aligned on the report's wording, with the distinction spelled out in the output. -
README linked to six pages that do not exist.
/introduction/,/installation/,/usage/,/features/extract/,/pipeline/and/citation/all 404 — they point at a docs layout reorganised out of existence. Absolute URLs, so neither mkdocs nor the internal link checker ever resolved them. -
tests/unit/test_cli_documented.pychecked 8 of 19 commands. Its regex matched onlyapp.command(name="x")and silently skipped the eleven registered as bareapp.command()(fn). It is now cross-checked against typer itself, which immediately surfaced a fifth undocumented command,validate. -
The Nextflow pipeline discarded its own validation artifacts, and never ran the cohort check. Three compounding gaps, so the scatter/gather validation added earlier was inert end to end:
-
run-allwrites{sample}.validation.jsonand{sample}.fingerprint.jsonon Parquet runs, but therunallmodule declared neither as an output — both were produced and left to die in the work directory. KREWLYZER_VALIDATE_COHORTexisted as a module and was included by nothing and called by nothing. The gather half of the check had no inputs and never ran.--strict-validationwas not exposed as a pipeline parameter at all.
The workflow now collects the fingerprints and runs the cohort step, and
emits cohort_report. Tool-level mode (--use_runall false) produces no
fingerprints, so the cohort step is skipped rather than run on a partial set,
which would report degeneracy that is an artefact of the missing samples.
output_format now defaults to parquet (was tsv). kreview reads
Parquet only, and a tsv-only run additionally skips per-sample validation
because the contract describes the Parquet surface — so the old default
produced a cohort that was both invisible downstream and unchecked, silently,
since every reader there swallows exceptions and yields an empty feature
dict. The pipeline warns at start whenever the selection still excludes
Parquet.
Breaking for anyone relying on the pipeline emitting TSV by default; pass
--output_format tsv or both to restore it.
Three parameters were also read by modules but declared nowhere:
validate_min_samples, silent and assay. Nextflow returns null for an
undeclared parameter rather than failing, so the ?: default swallowed it —
the parameters were invisible to --help and the docs, and would become hard
errors the moment nextflow.enable.strict is enabled. All three are now
declared and documented, and a test cross-checks every params.* read in a
.nf file against the config.
docs/nextflow/outputs.md gains the validation artifacts and a table of the
six output families whose values change in 0.9.0; parameters.md gains
--strict_validation, --validate_min_samples, --gc_correct,
--queue_size and the Parquet caveat.
FSC.regions.e1onlyselected by coordinate, and emitted a row for every gene whether or not it had one.filter_fsc_to_e1sorted bystartand took the first row per gene. It had no strand to work with — the temp BED handed to the Rust aggregator was five columns, so the asset's strand and E1 flags were discarded at that boundary regardless of what the asset contained.
On a real xs2 sample, 24 of 146 rows were the canonical exon 1; the rest were internal exons labelled E1.
FSC.regions now carries strand, is_e1 and is_alt_e1 end to end
(Python Region → 8-column temp BED → Rust GeneRegion → output), and
selection is on the flags: a region qualifies if it overlaps the canonical
transcript's exon 1 or another basic protein-coding transcript's first
exon. Both are real transcription starts, and requiring the canonical one
discards most of a panel — 25 of 128 xs1 genes against 40.
Genes with neither flag are now omitted rather than back-filled. The table drops from 146 rows to 48 on that sample, and every remaining row is on an annotated first exon.
A legacy or custom gene BED still works: the flags are absent, and selection
falls back to lowest-start with a warning that it is not strand-aware. .
rather than 0 marks an unknown flag, so "the source could not say" stays
distinguishable from "no".
Output-contract impact. FSC.regions gains three columns;
FSC.regions.e1only changes row count and membership. Neither is read by
kreview today (e1only is in NOT_CONSUMED), so nothing downstream breaks,
but any local analysis keyed on the old 146-row shape will.
-
MDS.generow order was non-deterministic. The writer iterated aHashMap, and Rust randomises hash iteration per process, so two runs on identical input produced byte-different files. A comment above the loop claimed "stable, reproducible output ordering" — it sorted regions within a gene, not the genes themselves.fsc.rsalready sorted its genes; this did not. Verified on a real sample: identical SHA-256 across two runs after the fix. -
region-mdsE1 was never strand-aware on panel data, and the new assets would have silently disabled it for WGS too. Two compounding problems.
The strand fix in identify_e1_regions() reads RegionInfo.strand, but
parse_gene_bed() hard-coded '+' for the panel format — the 5-column panel
BEDs have no strand column to read. So mds_e1 reported the last exon
for every minus-strand xs1/xs2 gene, which the docs claimed was fixed.
Worse, format detection keys on column count (5 → panel, 8 → WGS,
anything else → warn and treat as panel). Regenerating the assets took them
to 11 columns, so the WGS BED would have fallen into the panel branch and
lost its strand as well — turning a panel-only defect into a universal one.
The existing E1 tests build RegionInfo by hand and never call
parse_gene_bed, so none of them would have failed.
region-mds now recognises the annotated format, reads strand from it, and
consumes the precomputed is_e1 instead of re-deriving it — the build-time
flag comes from a GENCODE exon_number and is simply better than the
coordinate heuristic. Legacy formats still parse; the 5-column one now logs
a warning that E1 will not be strand-aware rather than quietly producing a
plausible number.
mds_e1 now distinguishes three states instead of two: a value, 0.0 for
"E1 exists but had no fragments", and NaN for "this gene has no E1 at
all". The last previously collapsed into 0.0 — the worst available
choice, since MDS lives in [0, 1] and lower means more abnormal, so a
fabricated 0.0 read as maximal tumour signal. It affects 88 of 128 xs1
genes, and every gene when a strandless legacy BED is supplied.
A legacy 5-column gene BED still parses and still produces per-exon MDS, but
region-mds now refuses to derive E1 from it — without strand the
heuristic returns the last exon for every minus-strand gene. mds_e1 is
NaN for that input, with a warning naming the fix.
Output-contract impact. MDS.gene.mds_e1 and mds_e1_z change for every
minus-strand gene on panel data, and for panel genes whose canonical exon 1
is not the most 5' captured region. Not comparable across the 0.9.0 boundary.
region-mds.mddocumented the wrong MDS scale, and contradicted itself. The Formulas section showed an unnormalised Shannon entropy and a "~6.0 to ~8.0" range — raw bits, which the tool has never emitted — while the clinical table further down the same page quoted ~0.95–1.0.motif_utils.rsdivides bylog2(256) = 8, so MDS has always been in[0, 1].
Anyone who built a threshold from the formula section is out by a factor of 8.
The divisor is now pinned in validate/claims.py, so the doc and the code
cannot drift apart again silently.
- Fragment coordinates were wrong whenever R1 was the rightmost mate. The
BED writer computed
pos() + |tlen|, correct only when R1 is leftmost. For the other orientation the interval was shifted right by roughlytlen - read_length— measured at 105bp on a real read.
Output change:
{sample}.bed.gzmoves for every affected fragment, and with it every positional feature — OCF end phasing, WPS fragment centres, TFBS/ATAC and gene/exon overlap, and the GC value in BED column 4. Fragment lengths are unchanged, so FSD/FSR/FSC size distributions are unaffected. Values from before this fix are not comparable; PON baselines built on uncollapsed input should be rebuilt.
Measured incidence: ~48% of R1 reads on an uncollapsed BAM, versus 0.4% after collapsing, 0.8% simplex and 0.0% duplex. Consensus callers normalise R1 to forward — Marianas (XS1) and fgbio (XS2) both — which is why this survived: the primary MSK-ACCESS path barely triggers it, while WGS and uncollapsed input run at the full rate.
The pre-Rust implementation was strand-aware; the 2025-12 rewrite dropped the
reverse branch while writing a correct sign-aware version forty lines below,
in the same function. No commit message, doc, comment or test ever mentioned
it. The same pattern was later copied into region_mds.rs, fixed here too.
Branches on the TLEN sign rather than is_reverse(): htslib gives the
rightmost segment a negative TLEN, and real data contains records flagged
forward whose mate lies to their left, where the two disagree.
-
The fragment BED could be emitted unsorted, which makes tabix indexing fail and takes every downstream feature with it. A consequence of the fix above: a fragment whose R1 is rightmost begins before the read that produced it, so read order stopped being coordinate order. The writer now sorts each chromosome before writing. Caught by the invariant suite, not by review.
-
Gene- and region-level FSC used different size bands than the genome-level table, so identically-named columns meant different things. The genome-bin counter split at
<=100 / <=149 / <=220 / <=260 / <=400with a sixthultra_longchannel. The gene/region aggregator split at<100 / <150 / <260 / <400with no sixth channel.mono_nucltherefore meant 150–220 bp inFSC.parquetand 150–259 bp inFSC.gene.parquet; a column labelleddi_nuclheld the genome table'slongrange; and genelongheld what the genome table calledultra_long. A consumer reading both tables into one feature matrix was combining channels that do not describe the same thing.
The bands were never a deliberate gene-level choice. Three places already
documented the correct (genome) bands for these tables: the column table in
docs/features/core/fsc.md, the column table in
docs/reference/output-files.md, and a line in the latter describing
"146 genes × 6 channels" against an implementation that emitted five. Only the
code disagreed, so this is a correction, not a redefinition.
Both paths now call a single SizeChannel::of, and the bounds exist in
exactly one place. Duplicating them is what allowed the drift, so the fix is
the deduplication rather than a second edit to the same numbers. Gene and
region tables gain ultra_long and ultra_long_ratio.
A fragment that passes the length filter but matches no channel is now counted and reported once at write time. It is not logged where it is detected: that point is inside the rayon loop, and routing a per-fragment warning through pyo3-log into Python deadlocks under the GIL — observed while testing this change, where it hung the run outright rather than reporting anything.
Output-contract impact. FSC.gene.* and FSC.regions.* change value for
every fragment in a moved band, and gain two columns. The six ratios sum to 1;
the five that existed before sum to 1 - ultra_long_ratio, and the contract
gate now asserts both so a downstream renormalisation cannot quietly use the
wrong base. Values from before this fix are not comparable. FSC.parquet and
FSC.ontarget.parquet are unchanged — the genome path was always correct.
--generate-jsonsilently dropped most features for compressed and Parquet runs. Every probe inFeatureSerializer.from_outputs()wasPath(f"{sample}.FOO.tsv").exists(), and that gate ran beforeread_table()(which does understand.tsv.gz/.parquet). Reconstructing a real MSK-ACCESS xs2 output directory (written with both--output-format bothand--compress, so every table exists as.parquetand.tsv.gz, but never as bare.tsv) shows the payload going from 5 to 16 feature families:
| features captured | |
|---|---|
| before | mfsd, ocf, wps, wps_background, wps_panel |
| after | the above plus fsd, fsr, fsc, fsc_gene, fsc_region, fsc_counts, motif, tfbs, atac, gc_factors, region_mds |
The three WPS entries survived only because they are probed as .parquet;
mfsd and ocf survived only because a separate cleanup defect leaves stray
uncompressed .tsv copies of those two outputs behind. Every fragmentomics
feature — the entire size, motif and regulatory signal — was absent from the
ML payload.
Note: this means tidying up those stray
.tsvleftovers without this fix would have madefeatures.jsonstrictly worse.
-
E1-only FSC was never generated for compressed or Parquet runs, and was misnamed when it was. Two compounding bugs:
unified_processorhard-codedoutputs.fsc_region = ...FSC.regions.tsv, so theoutputs.fsc_region.exists()guard was False whenever the real file was.tsv.gz/.parquetand E1 generation was skipped entirely; andfilter_fsc_to_e1()derived its output name viaPath.stem, which strips only the last dot-segment, so a.tsv.gzinput producedS1.FSC.regions.tsv.e1only.tsv.gzinstead ofS1.FSC.regions.e1only.tsv.gz. Addedstrip_table_extension()for compound suffixes. -
EndMotif1mer was unreadable whenever
--compresswas used.write_end_motif_1mer()gzipped the table viawrite_table()and then appended the# c_fraction/# entropy/# c_bias/# samplefooter with a plainopen(path, "a"). The result is a valid gzip member followed by raw bytes, whichgzip.decompress()and pandas both reject withBadGzipFile: Not a gzipped file (b'# ')— so the entire table was lost, not merely polluted with footer rows. Every run using--compress(i.e. any pipeline emitting.tsv.gz) is affected. The footer is now written through the gzip codec, andread_table()gained a recovery path that decompresses the first gzip member and skips trailing bytes, so files already produced by <= 0.8.3 remain readable. -
WPS: nucleosome repeat length (NRL) was a constant, not a measurement. The Alu background profile covered only the ~300bp Alu body in 30 bins, so the DFT period grid was
300 / iand the only index falling inside the nucleosomal search band wasi = 2(150bp). Every sample therefore producednrl_bp = 150.0,periodicity_score = 0.3333andadjusted_score = 0.0regardless of its data, and the documented "healthy NRL ~190bp, quality > 0.7" was unreachable. The profile is now 2000bp (850bp flank + 300bp body + 850bp flank) binned at 200 x 10bp, the DFT is zero-padded 8x, the search band is 140-250bp, and the deviation tolerance is the 50bp already documented (the code used 20bp, which forcedadjusted_scoreto 0 given the pinned 150bp NRL).
Output change:
{sample}.WPS_background.parquetnrl_bp,nrl_deviation_bp,periodicity_scoreandadjusted_scorebecome data-dependent. Earlier values are constants and must be discarded, not compared.
- UXM: the
X(mixed-methylation) class was unreachable. The CLI passedmethy_threshold = unmethy_threshold = 0.5. Because the backend evaluatesratio >= methy_thresholdfirst, every fragment collapsed intoMorUand the publishedXcolumn was identically0.0for every region of every sample. Thresholds are nowMETHY_THRESHOLD = 0.75/UNMETHY_THRESHOLD = 0.25, matching the documented Loyfer et al. (2022) definition.
Output change:
{sample}.UXM.tsvXbecomes non-zero andU/Mshrink correspondingly. Models trained on the previous columns must be refit.
- region-MDS: E1 (first exon) selection ignored strand.
identify_e1_regionsalways chose the lowest start coordinate, so for every minus-strand gene the reportedmds_e1was the gene's last exon — roughly half of a typical panel. E1 is now selected by transcription order (lowest start on+, highest on-).write_gene_outputno longer re-derives E1 by coordinate; it reads the strand-awareis_e1flag, somds_e1also stops silently falling through to the next covered exon when E1 has no fragments.
Output change:
{sample}.MDS.gene.tsvmds_e1/mds_e1_zchange for all minus-strand genes. Earlier values are not comparable.
-
region-entropy crashed when the PON lacked a matching baseline. The Rust
apply_pon_zscorereturns early without writing an output file if the PON has notfbs_baseline/atac_baselinetable. The Python caller then ranload_entropy_tsv()on the missing file, tripping itsassert df is not None, and the raw Rust output was unlinked afterwards — losing the entropy results entirely. It now degrades to the documented no-PON behaviour (z_score = 0) with a warning. -
Metadata footers parsed as data (
read_table).{sample}.EndMotif1mer.tsvappends# c_fraction,# entropy,# c_biasand# sampleafter the data rows, butread_table()calledpd.read_csvwithoutcomment="#". Those lines became data rows andFeatureSerializerpropagated them intofeatures.motif.edm_1meras junk keys with NaN values.read_table()now defaults tocomment="#"(callers may override).
Output change:
features.motif.edm_1merloses the four"# ..."keys.
-
The tool-level Nextflow path could not run at all.
modules/local/krewlyzer/extract/main.nfdeclaredpath("*.json")as a required output, but metadata moved from.metadata.jsonto.metadata.tsv/.parquetin 0.7.0. The tool has not written a.jsonsince, so the process failed on a missing required output. Unnoticed becauseuse_runall = trueis the default. Now emits the three real metadata variants as optional outputs, and the stub writes TSV to match. -
Rust test suite did not compile.
src/output_utils.rsusedtempfile::tempdirwith no[dev-dependencies]section inrust/Cargo.tomland referencedArcwithout importing it, socargo testfailed to build and no Rust unit test was running in CI. Also corrected the stalegc_referencebin-index test, which asserted a 68-bin ceiling althoughget_length_bin_indexspans 60..999 across 188 bins. -
--generate-jsoncrashed on any output directory without a metadata table. The metadata probe was the one call site that handed its resolver result straight toread_table()without anis not Noneguard, raisingAttributeError: 'NoneType' object has no attribute 'suffix'. The format-matrix test always wrote a metadata table, so it never reached the branch; it now covers the absent case across tsv/tsv.gz/parquet. -
E1-only FSC and FSC PON z-scoring were still skipped on compressed and Parquet runs.
aggregate_by_gene()returned the pre-conversion.tsvpath even aftercleanup_intermediate_tsv()had deleted it, so every caller guarding on.exists()saw a missing file. It now returns the path it actually wrote, which fixes the E1 guard,apply_fsc_region_ponandapply_fsc_gene_pontogether;outputs.fsc_genewas additionally never resolved at all. The two "✓ … (with PON z-scores)" log lines no longer claim success when z-scoring bailed. -
Removed
src/krewlyzer.libs, a broken symlink into a build sandbox (/usr/local/lib/python3.11/dist-packages/) committed by mistake. -
Restored
blackandruffcleanliness. Both pass ondevelopbut were broken by the changes above this line, and CI's lint job gates on them.
Documentation
-
Corrected feature documentation against the implementation: the inverted "higher MDS = tumor" summary in
motif.md; mFSD LLR sigmas (human()is 167/30 and 145/35, not 167/35 and 145/25), the canine/ssdna presets, their unreachability, and theMIN_FOR_KS = 2floor; the non-existent 20-500bp region-entropy window and its unusable absolute entropy bands; OCF's 1bp resolution and per-10,000 per-label normalization; the FSR channel table (FSR and FSC share one counter) and the fact thattotalspans 65-1000bp while only five channels are returned, sochannel / totalratios do not sum to 1; and that WPSz_vector/shape_score/z_amplitudeare documented but never emitted, whileapply_pon_zscoresilently yieldsz = (0 - mean) / std. -
Corrected three doc statements the pass above missed:
fsr.py's module docstring still listed the retired 65-99/100-149/150-259/260-399/400+ channel table (FSR reads FSC's per-bin counts, so it uses FSC's channels), and now also records thatshort_frac/long_fracdivide by atotalspanning 65-1000bp;estimate_periodicity's doc comment still claimed a 140-200bp search band after the NRL fix widened it to 140-250bp; andmfsd.mdstill asserted a KS floor of ≥5 in its tip and worked example, contradictingMIN_FOR_KS = 2and its own column reference.
Note: those retired FSR boundaries are not fictional —
GeneResultinfsc.rsstill uses exactly that split, somono_nuclmeans 150-259bp inFSC.gene/FSC.regionsbut 150-220bp inFSC. Tracked separately.
[0.8.3] - 2026-04-28
Fixed
- GIL Deadlock in Parallel Modules: Wrapped all Rayon
par_iter()calls inpy.allow_threads()acrossmfsd.rs,uxm.rs,extract_motif.rs, andregion_mds.rs. Previously,pyo3-logattempted to acquire the GIL from Rayon worker threads while the main thread held it, causing indefinite hangs on multi-threaded runs (observed as 16-hour wall time on IRIS HPC). - mFSD Silent Error Swallowing: Replaced
Err(_) => continuein BAM record iterator with a logged error breaker (max 1000 consecutive errors). Previously, corrupt BAM regions could cause infinite silent loops. - FILTER_MAF Substring Match: Changed
sample_id in tsbto exact match (sample_id == tsb) preventingT01from matchingT01-XS1. - FILTER_MAF Comment Lines: Stripped
#comment lines from filtered MAF output in both multi-sample and single-sample modes. - mFSD 0-Variant Guard: Added early exit at both Python (
wrapper.py) and Rust layers when MAF has 0 data lines. Produces header-only TSV instead of attempting BAM access. - mFSD GC Correction Fallback: When reference FASTA is unavailable or GC lookup fails for a region, GC correction is now skipped (weight=1.0) instead of silently using a hardcoded 50% GC.
Added
- mFSD BAM I/O Diagnostics: Per-variant timing, BAM open/fetch latency logging,
record counts, and slow-variant warnings (>30s). All diagnostics use structured
logmacros (visible with--verbose). - mFSD Header Constant: Extracted 46-column TSV header into
MFSD_HEADERmodule constant, eliminating duplication between normal output and 0-variant early-exit.
Changed
- Minimum Python:
requires-pythonraised from>=3.8to>=3.10(Python 3.8 and 3.9 are EOL). - Minimum Rust: Added
rust-version = "1.87"MSRV toCargo.toml. - Cargo Dependencies: Updated 81 semver-compatible transitive dependencies (rayon 1.11→1.12, anyhow 1.0.100→1.0.102, flate2 1.1.5→1.1.9, etc.).
- Diagnostic Logging: Converted all
eprintln!breadcrumbs inside parallel closures to structureddebug!/warn!macros for proper timestamp/level formatting through Python's logging framework.
Documentation
- Developer Guide: Added "PyO3 + Rayon GIL deadlock" Known Gotcha with correct and incorrect code examples.
- Developer Guide: Updated contributing checklist test count (244→357).
Tests
- Added
test_mfsd_zero_variants— verifies 0-variant input produces header-only TSV. - Added
test_mfsd_maf_with_comment_lines— verifies MAFs with#comment headers are parsed correctly. - Total: 357 passed, 4 skipped.
[0.8.2] - 2026-03-26
Fixed
- Region MDS Output Collision: Fixed
Path.with_suffix()compound-extension bug that causedMDS.exon.tsvandMDS.gene.tsvto both resolve toMDS.tsv, silently destroying per-exon and per-gene data. Affects standaloneregion-mds,motif, andrun-allcommands. - Motif Tracking Paths: Fixed 5 additional
with_suffix()bugs in motif output tracking (EndMotif, BreakPointMotif, MDS, MDS.ontarget, EndMotif1mer) that caused path collisions and silent PON z-score normalization failures. - Gene Format Conversion: Region MDS gene output now correctly undergoes format conversion (Parquet/gzip) when no PON is provided.
- MDS Z-score Logging: Z-score append failures upgraded from
debugtowarninglevel to surface issues in production runs.
Added
- Compound Extension Tests: 57 new parametrized tests (
test_compound_extension.py) covering all 13 compound base names across TSV, Parquet, both formats, and roundtrip. - Developer Guide: "Known Gotchas" section documenting the
Path.with_suffix()anti- pattern with safe alternative and complete table of affected names.
Changed
- Nextflow Outputs: Added missing Parquet emit declarations for TFBS sync (2), ATAC sync (2), ATAC/TFBS ontarget sync (2), and fsc_counts (4 — TSV+Parquet, off/on-target).
- SLURM Script: Tuned head process memory (16G→32G), queue size (100→200), added
--output_format both --compress_tsv true --verbose truefor 14K+ sample production runs.
[0.8.1] - 2026-03-24
Fixed
- WPS Panel GC Correction: Panel WPS now uses on-target GC correction factors instead of off-target. Panel anchors overlap capture regions, making on-target correction more accurate. Falls back to off-target factors when on-target unavailable.
- Rust WPS Background: Fixed
coitreesmetadata type mismatch in WPS background consumer (wps.rs:1684) — uses.clone()for cross-platform compatibility (returns&usizeon macOS butusizeon CI Linux).
Changed
- Dead Code Cleanup: Removed 11 unused Python z-score functions from
pon_integration.py(all replaced by Rust equivalents). Module reduced from 448 to 98 lines, 14 to 3 functions. Remaining:load_pon_model,compute_nrl_zscore,compute_periodicity_zscore. - GC Factor Resolution:
gc_str/gc_ontarget_strpath resolution hoisted to a shared section inunified_processor.py, eliminating duplication between panel WPS and TFBS/ATAC.
Added
- build-pon Logging: OCF on-target/off-target baseline status now logged during PON build.
- sbatch Script:
scripts/build_pon_unfiltered.shfor building PON from high-coverage unfiltered BAMs on SLURM clusters.
Data
- PON Models Updated: Rebuilt xs1/xs2 PON models for both
all_uniqueandduplexvariants with krewlyzer 0.8.x.
[0.8.0] - 2026-03-17
Added
- On-target PON Z-scores: Panel mode now computes on-target/off-target PON baselines for MDS, OCF, and FSD features, providing clinically-scoped z-scores in panel assays.
- FSR On-target Output: FSR now emits a separate
.FSR.ontarget.tsvfile in panel mode. - FSR Real Genomic Coordinates: Region labels now reflect true genomic window coordinates instead of internal indices.
Fixed
- Output Format / Compress: All output files now correctly respect
--output-formatand--compressflags (gzip compression, path handling, GC correction factors loading). - Rust GC Correction: Added missing
PathBufimport; fixed path handling for correction factor files. - Rust Clippy: Replaced manual string strip with
strip_suffix(manual_striplint).
Documentation
- Corrected stale FSR column names in
concepts.mdandjson-output.md. - Updated FSR on-target output docs and window/panel-mode descriptions.
- Documented
_core.pyistub maintenance requirements.
CI
- GitHub Actions Node.js 24: Bumped
actions/checkoutv4→v5,actions/setup-pythonv5→v6,actions/cachev4→v5 to address Node.js 20 deprecation (enforced June 2, 2026).
[0.7.0] - 2026-03-02
Added
- Configurable Output Formats:
--output-format tsv|parquet|both(default:tsv) controls all tabular feature outputs.--compressgzip-compresses TSV outputs (.tsv.gz). WPS outputs remain always-Parquet regardless of setting.
Fixed
- build-pon Intermediate Files: Explicitly force
output_format="tsv"andcompress=Falsein allprocess_sample()calls withinbuild-pon(both parallel and sequential paths). Prevents silent failure if default output format changes — intermediate files are internal scratch consumed bypd.read_csv(sep="\t"). - Feature Serializer: Include
mds_zin JSON output for thefrom_outputs()code path. - OCF Base File:
OCF.tsv= all reads (on + off combined), not off-target.OCF.offtarget.tsvis the true panel off-target score. Corrected in docs and code comments. - Rust wps.rs: Remove erroneous
*dereference onnode.metadata(E0614 —usizeis Copy). - Rust gc_correction.rs: Prefix unused
valid_regions_pathparam with_to silence compiler warning; parameter retained for API symmetry. - MkDocs Snippets: Fix
--8<-- "CHANGELOG.md"/--8<-- "CONTRIBUTING.md"broken includes by changingpymdownx.snippets.base_pathfrom scalardocsto list['.', 'docs']so repo-root files resolve without path traversal blocked by CI deploy sandbox.
Documentation
- Post-0.6.0 Docs Audit (12 files, 7 issue categories):
- Fixed broken
--8<-- "../CHANGELOG.md"/--8<-- "../CONTRIBUTING.md"snippet includes - Corrected outdated "no global
--output-formatflag" note incli/run-all.md - Updated
metadata.json→metadata.tsvacross 7 files (8 references total) - Added WPS always-Parquet exception note to
reference/output-files.md - Added build-pon intermediate TSV format note to
guides/building-pon.md - Updated test count (248 → 244 + 4 skipped) and added CI lint steps to
developer-guide.md - Added
--output_formatand--compressparameters tonextflow/parameters.md - Output Format Options section: New section in
reference/output-files.mddocumenting--output-format,--compress, WPS always-Parquet exception, and--generate-json - OCF Variant Clarification: Added 3-variant table and note block in
reference/output-files.mdexplainingOCF.tsv(all reads) vsOCF.ontarget.tsvvsOCF.offtarget.tsv - docs/index.md: Replaced
:latestDocker tag with explicit:0.7.0per release policy
CI
- Lint Job: Added parallel
lintjob (ruff · black · mypy · cargo clippy -- -D warnings) running alongside tests on all push/PR events
[0.6.0] - 2026-02-28
Added
- mFSD Base Quality Filtering:
--min-baseq/-Q(default 20) gates variant evidence by base quality - mFSD GC Correction: Rust-native LOESS GC bias correction for variant fragment size distributions
- mFSD Duplex Weighting: Proper consensus fragment handling via
--duplex - Region MDS
--sample-name: Consistent output naming without post-hoc rename - Feature Serializer: Auto-load
fsc_counts,region_mds,uxminfrom_outputs() - IRIS Batch Submission:
scripts/run_krewlyzer_iris.shfor SLURM/IRIS cluster runs with--generate_json - nf-core Institutional Configs:
custom_config_baseparam and IRIS profile support - Versioned Documentation: Implemented
mikefor dev/stable doc versions - Nextflow mfsd Module: Full standalone params (
--reference,--correction-factors,--mapq,--minlen,--maxlen,--min-baseq,--duplex,--no-skip-duplicates) - Nextflow runall:
fsc_counts.tsvoutput declaration,--min-baseqwired - mFSD Integration Tests: 161 lines of new test coverage
Fixed
- mFSD MAF Parsing: Header-based column lookup (fixes column-index mismatch with different MAF flavors)
- Nextflow FILTER_MAF: Complete overhaul — eliminated join operator blocking, replaced regex with substring match, fixed SyntaxError in versions.yml, dynamic maxForks for SLURM
- Nextflow Workflow Streaming: Fixed RUNALL blocking from
remainder:true,failOnMismatch, channel round-robin; usedmultiMapfor proper routing - Nextflow RUNALL Outputs: Added 14 output declarations, fixed BreakPointMotif casing, explicit publishDir
- Region MDS Nextflow:
--sample-namereplacesmvworkaround - Nextflow Config: Executor queueSize placement,
-qsCLI flag, global publishDir removal - WPS CLI Tests: Fixed
--inputflag (was positional arg) — recovered 2 skipped tests - Pandas FutureWarning: Fixed
pd.concatwith all-NA columns in PON test fixture
Changed
- Code Quality: Black formatted 71 files, ruff fixed 129 lint errors, cargo clippy applied
- Ruff Config: Added
[tool.ruff]topyproject.tomlwith documented E402/F821 ignores - Agent Config:
.agent/→.agent/rules/withalwaysApplyfrontmatter
Documentation
- 45-item Audit: Corrections across 25 doc files including
.csv→.tsv(7 files),.WPS.tsv.gz→.parquet(3 files), phantom--output-formatremoved, Docker versions→X.Y.Z, parameters.md 12→28, outputs.md 14→41, JSON schema corrected, developer guide Rust table 10→19, architecture pipeline signature updated - PDF Embedding: Fixed rendering with mkdocs-pdf plugin
[0.5.3] - 2026-02-06
Documentation
- Admonition Formatting: Converted blockquotes to proper MkDocs admonitions across 10+ docs files for consistent styling
- Table Rendering: Fixed tables not rendering by adding required blank lines before table headers
- LaTeX Formulas: Fixed underscore escaping in math blocks (mFSD, FSR, WPS formulas)
- Abbreviations: Added glossary with auto-append for consistent acronym tooltips (cfDNA, WPS, etc.)
- Docs CI: Added PR validation trigger with
mkdocs build --strictto catch issues before deployment
[0.5.2] - 2026-02-05
Added
- Dual BAM Support (mFSD): New
--mfsd-bamparameter inrun-allto use a dedicated duplex BAM for mFSD while other features use the main all_unique BAM - Auto-duplex detection from filename patterns (
-duplex,_duplex,.duplex) - Startup banner shows mFSD BAM path when specified
- Nextflow: Added
mfsd_bamcolumn to samplesheet schema
Fixed
- Correction Factors Loading: Fixed delimiter detection for
.correction_factors.tsvfiles (CSV format with TSV extension)
Documentation
- Release Guide: Added version format note (use
0.5.2notv0.5.2in code) - mFSD docs: Updated with dual BAM workflow examples
- Samplesheet docs: Added
mfsd_bamcolumn documentation - Footer: Added site footer with author attribution and Antigravity acknowledgment
- Theme: Changed color scheme to blood-red (#ef5552) to match "krew" (blood) branding
[0.5.1] - 2026-02-04
Fixed
- Docker Image: Data folder now bundled in container (was missing in 0.5.0)
- CI Tests: Tests requiring bundled data now skip in PyPI installs via
@requires_datamarker
Changed
- Docker Tags: Use versioned tags only (e.g.,
:0.5.1); no:latesttag published
Documentation
- Installation: Added Singularity/Apptainer section for HPC clusters
- Installation: Clarified Docker uses release tags, not
:latest - Nextflow Examples: Added IRIS HPC cluster profile example (
-profile iris) - Testing Guide: Added Data Availability section explaining PyPI vs source
- Structure: Moved Changelog to Development section; removed duplicate files
- Mermaid Diagrams: Upgraded to
mkdocs-mermaid2plugin for theme-aware rendering - Removed: Full-site PDF export (unreliable rendering)
[0.5.0] - 2026-02-02
Added
PON Framework
build-poncommand: Generate PON models from cohort samples with FSD/WPS/OCF/MDS baselines- Bundled PON assets: Pre-computed xs1/xs2 PON models for
all_uniqueandduplexvariants --pon-variantflag: Select betweenall_unique(max coverage) orduplex(highest accuracy) PON models- Duplex tag warning: mFSD warns when
--duplexis used but no cD tags found in BAM
Panel Mode (MSK-ACCESS)
- Assay-aware asset resolution: Auto-load targets, anchors, PON via
-A/--assay xs1|xs2 - On/Off-target splitting: All tools produce separate
.ontarget.tsvoutputs - Gene-level FSC: New
{sample}.FSC.gene.tsvand{sample}.FSC.regions.e1only.tsv - Bait edge masking: WPS
--bait-paddingto remove capture artifacts --skip-target-regions: Force WGS mode for panel assays
Duplex Sequencing
--duplexflag (mFSD): Enable family size (cD tag) weighting for duplex BAMs- LLR scoring: Log-likelihood ratio classification for cross-species support
- GC-weighted mFSD: 5 new GC-corrected mean fragment size columns
Region-Based Analysis
- Region Entropy: TFBS/ATAC dual-output architecture with per-cluster entropy
- Region MDS: Per-gene Motif Diversity Score with E1-only filtering
- Rust backend: New
region_entropy.rsfor high-performance calculation
Feature Enhancements
- Jagged Index: 1-mer End Motif analysis with C-end fraction
- WPS v2.0: Hierarchical stacking, extended anchors, panel-specific normalization
- Output formats:
--formatflag (tsv/parquet/json) for all tools --generate-json: Unified JSON export with all features
Infrastructure
- Git LFS: Large files (.gz, .parquet, .bed) tracked via git-lfs
- BGZF compression: All BED outputs use block-gzip format
- GC references: Pre-computed Parquet format for GRCh37/GRCh38
- Unified processor: Single-pass Rust pipeline for Extract→Motif→FSC→FSD→WPS→OCF
Changed
Nextflow Pipeline
- nf-core compliance: Refactored to shared INPUT_CHECK subworkflow
- KREWLYZER_RUNALL: Unified module with 30+ output channels
pon_variantparameter: New pipeline parameter (default:all_unique)
CLI
- 9 tools with
--pon-variant: run-all, fsc, fsd, fsr, wps, ocf, motif, region-entropy, region-mds - Centralized asset resolution:
resolve_pon_model()andresolve_target_regions() - Filter flags:
--mapq,--minlen,--maxlen,--skip-duplicates,--require-proper-pair - Parallel processing:
--threadsoption across all tools
Data & Assets
- Terminology: Renamed "blacklist" to "exclude_regions" for inclusive language
- Directory structure:
data/{type}/{genome}/{variant}/for organized assets - Memory optimization: Parallel sample processing with spawn context in build-pon
Fixed
- mFSD: Filter discordant reads with extreme TLEN values
- mFSD: Fix verbose mode hanging by moving debug logging outside parallel loop
- Proper pair detection for legacy BAMs without proper pair flags
- GC correction for gene-level FSC in panel mode
- CIGAR handling improvements for INDELs and complex variants
- BAM extraction for v1 ACCESS and bgzip output format
Documentation
- 11 docs updated: PON Variant Selection across pon.md, pipeline.md, usage.md, feature docs
- Release Guide: New
docs/advanced/release-guide.mdfor version release process - Math rendering: LaTeX formulas in all feature documentation
- Glossary: New terms and definitions
Tests
- 28 new test files: Unit, integration, and e2e coverage
- 4 PON variant tests: test_asset_resolution.py
- Real data tests: test_real_data.py for end-to-end validation
[0.3.2] - 2025-12-18
Fixed
- CI Build: Removed
gfortranandscikit-misc- GC correction now fully in Rust - FSC GC Correction: Added
--gc-correct/--no-gc-correctand--verboseflags to FSC CLI - Python LOESS Removed: Removed
gc_correct()fromhelpers.pyandpostprocess.py
Changed
- Dockerfile: Removed
gfortrandependency (no longer needed) - GitHub Actions: Removed
gfortranfrom CI workflows - Nextflow Modules: Updated all container versions to
0.3.2
[0.3.1] - 2025-12-18
Added
- Rust LOESS GC Correction: New
rust/src/gc_correction.rsmodule using thelowesscrate - Per-fragment-type correction (short, intermediate, long)
- Configurable LOESS parameters (fraction, iterations, delta)
- FSR GC Correction:
--gc-correct/--no-gc-correctflag (default: True) - Uses corrected counts from Rust before calculating ratios
--verboseflag for detailed logging- WPS GC Correction:
--gc-correct/--no-gc-correctflag (default: True) --reference, -r: Reference FASTA for computing region GC content- FASTA-based GC computation using rust-htslib::faidx
- Graceful fallback if reference not provided
Changed
- FSC Rust Backend: Added
count_fragments_gc_corrected()function for integrated GC correction - WPS Rust Backend: Updated
calculate_wps()withreference_path,gc_correct,verboseparameters - lowess Dependency: Updated from 0.3 to 0.4 for improved API
Documentation
- Updated
docs/features/fsr.mdanddocs/features/wps.mdwith GC correction options
[0.3.0] - 2025-12-16
Added
- mFSD Variant Type Support: Now handles all small variant types:
- SNV (single nucleotide)
- MNV (multi-nucleotide)
- Insertion
- Deletion
- Complex (substitution + indel)
- 4-Way Fragment Classification: Fragments classified as REF, ALT, NonREF, or N (low quality)
- Comprehensive Statistics: 6 pairwise KS comparisons (ALT-REF, ALT-NonREF, etc.)
- Derived Metrics: VAF_Proxy, Error_Rate, N_Rate, Size_Ratio, Quality_Score
- Quality Flags: ALT_Confidence (HIGH/LOW/NONE), KS_Valid
- Distributions Output:
--output-distributionsflag generates per-variant size distributions - Verbose Logging:
--verboseflag for debug-level logging with variant type breakdown - MRD Support: Proper handling of low fragment counts (1-2) common in MRD settings
Changed
- BREAKING: mFSD output format changed from 11 columns to 39 columns
- Nextflow MFSD: Distributions and verbose logging enabled by default
- Fragment Counting: Now counts fragments (R1 only) instead of reads to avoid double-counting
Fixed
- CIGAR Handling: Improved sequence extraction for INDELs and complex variants
[0.2.3] - 2025-12-16
Changed
- Nextflow Pipeline: Added
FILTER_MAFmodule for multi-sample MAF filtering: - New
mafandsingle_sample_mafcolumns in samplesheet - Filters MAF by
Tumor_Sample_Barcodematching sample ID (regex:.*sample_id.*) single_sample_maf=truebypasses filtering for per-sample MAFs- Skips MFSD when filtered MAF has zero variants (with warning)
- Memory-efficient streaming for large MAF files (100s of MBs)
- Outputs filtered MAF + mFSD results
- Nextflow Modules: Updated all container versions to
0.2.3
[0.2.2] - 2025-12-15
Changed
- Project Structure: Migrated to recommended maturin "src layout" for better Python/Rust separation:
- Python code moved from
krewlyzer/tosrc/krewlyzer/ - Rust code moved from
krewlyzer-core/torust/ - Single
rust/Cargo.toml(removed duplicate rootCargo.toml) - Rust module now imports as
krewlyzer._core(private) - Dockerfile: Rewritten as multi-stage build for smaller image size (~200MB vs ~1GB); amd64 only
- OCI Labels: Added container metadata for GitHub Container Registry
Fixed
- Distribution Compatibility: Updated release workflow to build
manylinux_2_28wheels (compatible with RHEL 8+, AlmaLinux 8, etc.) - Source Builds: Included
clangandllvm-develin build environment forbindgen/hts-sys - Docker Build: Added
patchelfto maturin installation for Linux wheel building - Test Imports: Updated test files to use new
krewlyzer._coreimport path - CI OpenSSL: Use system OpenSSL (
OPENSSL_NO_VENDOR=1) instead of building from source - CI Python Versions: Build wheels for Python 3.10, 3.11, 3.12 in release workflow
- Documentation: Updated logo paths for new
src/layout
[0.2.1] - 2025-12-15
Fixed
- Rust Compilation: Resolved cross-platform build issues with
coitreesmetadata types (usizevs&usize) by using explicit.to_owned()conversion. - CI Build: Added
gfortran,clang, andlibclang-devto CI workflows to fixscikit-miscandrust-htslibcompilation failures. - Permission Errors: CI scripts now robustly handle
sudopermissions when installing system dependencies.
[0.2.0] - 2025-12-12
Added
- Unified Engine: New high-performance Rust core (
krewlyzer-core) that processes Extract, Motif, FSC, FSD, WPS, and OCF in a single parallelized pass. - Fragment Extraction: dedicated
extractcommand (via Rust) to convert BAM to BED with configurable filters. - Extract Documentation: New
docs/features/extract.mddetailing extraction logic and JSON metadata. - Calculation Details: Comprehensive formulas and interpretation guides added to all feature documentation.
- Root Cargo.toml: Added to support standard
maturinbuilds for the hybrid Python-Rust package.
Changed
- Performance: Significant speedup (3-4x) for end-to-end analysis due to multi-threaded processing and single-pass I/O.
- Build System: Migrated to
maturinbackend for robust Rust extension compilation. - CLI (
run-all): Now defaults to the Unified Engine. - CLI Filters: Added
--mapq,--minlen,--maxlen,--skip-duplicates,--require-proper-pairflags torun-all,extract, andmotif. - Motif Outputs: Renamed output files to use
.tsvextension consistently (e.g.,{sample}.EndMotif.tsv). - Data Handling:
motifnow uses the unified engine, eliminating the need forbedtoolsbinary entirely. - Documentation: Updated
README.md,usage.md, andpipeline.mdto reflect the new workflow.- Corrected
pipeline.mdsamplesheet format documentation to match Nextflow schema. - Updated
usage.mdand feature docs to correctly specify output directory arguments.
- Corrected
Fixed
- Test Suite: Cleaned up
tests/directory, removing obsolete scripts and fixing integration tests (test_science.py,test_run_all_unified_wrapper.py). - Validation: Fixed BAM header issues in tests.
Removed
- Legacy Python Backends: Removed pure Python implementations of
extract,motif,fsc,fsd, ensuring all paths use the unified Rust core. - Validation Artifacts: Deleted temporary validation scripts and data.
[0.1.7] - 2025-11-26
Fixed
- PyPI Metadata: Added
readmeandlicensefields topyproject.tomlto ensure the long description is correctly displayed on PyPI.
[0.1.6] - 2025-11-26
Fixed
- Docker Build: Removed
libatlas-base-devdependency fromDockerfileas it is not available in thepython:3.10-slim(Debian Trixie) base image.
[0.1.5] - 2025-11-26
Fixed
- Docker Publishing: Switched to
GITHUB_TOKENfor GHCR authentication to fix permission issues. - PyPI Publishing: Added
skip-existing: trueto handle existing versions gracefully. - CI/CD: Added build checks for Python package and Docker image to PR workflows.
[0.1.4] - 2025-11-26
Fixed
- Test Dependencies: Removed unused
pybedtoolsimports fromfsr.py,fsd.py,uxm.py, andfsc.pywhich were causingImportErrorin CI environments wherepybedtoolsis not installed.
[0.1.3] - 2025-11-26
Changed
- Dependency Reduction: Removed
pybedtoolsdependency. - Refactor:
motif.pynow usespandasfor blacklist filtering and sorting, removing the need forbedtoolsbinary. - CI/CD: Added
pytestandpytest-mocktotestoptional dependencies inpyproject.toml.
[0.1.2] - 2025-11-26
Added
- Mutant Fragment Size Distribution (
mfsd): New module to compare fragment size distributions of mutant vs. wild-type reads using VCF/MAF input. - Enhanced Fragment Size Ratios (
fsr): Added "Ultra-Short" (<100bp) ratio bin. - Documentation: Comprehensive MkDocs website (
docs/) with material theme. - Pipeline:
run-allcommand now supports--variant-inputformfsdanalysis. - Nextflow: Pipeline updated to support optional variant input in samplesheet.
Changed
- Updated
README.mdto point to the new documentation site. - Added
mkdocsandmkdocs-materialas optional dependencies.