feat(observations): configurable observe_molecules; fix the secondary/supplementary filters (#82) - #83
Merged
Merged
Conversation
…uses The only way to change a setting was to build a whole GbcmsDnaConfig, whose four required fields -- variant_file, bam_files, reference_fasta, output -- have ZERO overlap with the seven this entry point reads. Two of them must name files that exist on disk, bam_files is a dict, and output is an OutputConfig. So adjusting one read filter meant fabricating an output directory and a variant path that are never touched: paths that read as load-bearing and are not, which is how the next person gets misled. Adds filters/quality/alignment/umi_tag/threads/apply_baq/library_type as keyword arguments. config= keeps working; an individual argument overrides the matching config field, so a pipeline config can be reused and adjusted. Found by the first real consumer (mulligan) needing a non-default read filter. Tests written inside gbcms could not have found it -- they never had to construct the config a caller does. umi_tag is the one argument where None is a choice rather than an absence: it means group by read pair, not UMI family. It uses a sentinel so that passing umi_tag=None alongside a config that sets one OVERRIDES it. Treating None as "not supplied" there would silently inherit the config's tag and change what counts as a molecule, with nothing to indicate it. Also corrects advice this module previously gave. The filters docstring told callers to set supplementary=False for cross-locus phasing. That does nothing: supplementary and secondary alignments are dropped before fragment evidence is built, so they never produce an observation. Measured across all six read filters on real data -- duplicates, improper_pair and indel work as documented; those two cannot change the result at all (#82). The new tests use `duplicates` as the discriminating filter for exactly that reason: `supplementary` would have been the intuitive choice and would have made the test pass vacuously. Tests: 420 passed (+11). ruff, black, mypy, clippy -D warnings, cargo test (215) all clean. Docs build warnings unchanged (1, pre-existing). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…n turned off Closes #82. Both flags previously could not change ANY output. An unconditional skip in the counting loop discarded those records before fragment evidence, "regardless of the user filter flags" -- redundant when the flags were on (such records never reach the cache) and defeating them when off. Measured across all six read filters on real data, these were the only two that did nothing; duplicates, improper_pair and indel all worked. Default runs are byte-identical. Both still default to on, so the records never reach the cache: total depth unchanged at 1202 across 40 real loci, and binned<->legacy parity holds. Turning a filter off now admits those records to FRAGMENT-level evidence (dpf/rdf/adf) and to the observation export, while read-level dp/rd/ad keeps its promise and excludes them. That split is not a compromise, it is the correct line, and the existing tests say so. A supplementary shares a QNAME with its primary and is the same physical read, so counting both toward dp reports depth 2 where one read exists -- test_supplementary_shared_qname_not_double_counted pins exactly that case, and two more tests assert the same read-level contract. Honoring the flag naively broke all three. Fragments are immune either way: hash_molecule keys on QNAME, so a primary and its supplementary collapse into one FragmentEvidence, which is why admitting them there cannot double-count. What this fixes concretely: a locus reached ONLY by a supplementary segment reported dpf=0 -- not a filtered read but a wrong answer, since a molecule demonstrably covers it. That is what made a molecule spanning a large deletion invisible to cross-locus phasing. A new test covers that arrangement; it is the other half of the double-count test's contract. New consequence, documented in the read-filters page, the observations page and the CLI help: with the filter off an admitted supplementary raises dpf while leaving dp unchanged, so the two stop moving together. The observation reconciliation invariant (rows == dpf) is verified under both settings. Also documents the improper-pair foot-gun found while auditing: MSK-ACCESS v1 alignments are 100% paired but 0% PROPER_PAIR, so enabling that filter there discards every read. Tests: 421 passed (+1). clippy -D warnings, cargo test (215), mypy, ruff, black, mkdocs --strict all clean. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Two commits, reviewable separately. Closes #82.
1.
feat:observe_moleculestakes the settings it actually usesobserve_moleculesreads seven settings.GbcmsDnaConfigrequires four. They do not overlap at all:GbcmsDnaConfigvariant_file,bam_files,reference_fasta,outputobserve_moleculesfilters,quality,alignment,threads,apply_baq,umi_tag,library_typeAnd it is not merely four extra keyword arguments —
variant_fileandreference_fastamust name files that exist on disk,bam_filesis adict, andoutputis anOutputConfig. Changing one read filter meant fabricating paths that read as load-bearing and are not.config=keeps working; an individual argument overrides the matching field.umi_tagneeds a sentinel. It is the one argument whereNoneis a choice — "group by read pair, not UMI family". Under the uniform "None = not supplied" rule, passingumi_tag=Nonealongside a config that sets one would silently inherit it and change what counts as a molecule. Tested.2.
fix: the secondary/supplementary filters actually do something nowBoth flags previously could not change any output — an unconditional skip discarded those records before fragment evidence, "regardless of the user filter flags". I audited all six read filters on real MSK-ACCESS data:
filter_duplicatesfilter_improper_pairfilter_indelfilter_secondaryfilter_supplementaryfilter_qc_failedThe line this draws
Simply honoring the flags breaks three existing tests, and they are right to break.
test_supplementary_shared_qname_not_double_countedbuilds a primary and its supplementary over the same locus and assertsdp == 1: they are one physical read, and counting both reports depth 2. Two more tests assert the same read-level contract, matching the CLI help.So the correct split is read-level vs fragment-level, not on/off:
dp/rd/addpf/rdf/adfFragments are immune to double-counting either way:
hash_moleculekeys on QNAME, so a primary and its supplementary collapse into oneFragmentEvidence.What it fixes
A locus reached only by a supplementary segment reported
dpf=0— not a filtered read but a wrong answer, since a molecule demonstrably covers it. That is what made a molecule spanning a large deletion invisible to cross-locus phasing. New test covers that arrangement; it is the other half of the double-count test's contract.Verification
Default runs are byte-identical: total depth unchanged at 1202 across 40 real loci, and binned↔legacy parity passes (18 tests).
New consequence, documented
With the filter off, an admitted supplementary raises
dpfwhile leavingdpunchanged — the two stop moving together. Recorded in the read-filters page, the observations page, and the CLI help.Also documented the foot-gun found while auditing: MSK-ACCESS v1 alignments are 100% paired but 0% PROPER_PAIR, so
--filter-improper-pairdiscards every read there.Gate
pytest— 421 passed; parity suite greencargo test— 215 passed;clippy --all-targets -D warningscleanruff check .,black --check,mypy— cleanmkdocs build --strict— 1 warning, unchanged fromdevelop(pre-existing)🤖 Generated with Claude Code