fix: add-samples wrote to a flat Parquet the query engine never reads - #37
Merged
Conversation
The rule that a chromosome's bucket directory takes priority over its flat Parquet file was reimplemented in every module that opens a variant file: the query engine, the dump and annotate workers, the benchmark helper, the compactor and the updater. Eight copies, none of them shared, and no single place to look when a writer and a reader disagree about where a row belongs. This adds the resolver as a standalone module. Nothing is migrated onto it yet, so the change is a pure addition; the call sites move over separately. partitioned_chroms() deliberately accepts any subdirectory rather than validating names against the canonical chromosome list. A database may hold unplaced or alt contigs — normalize_chrom() preserves bodies like 'chrGL000209.1' and 'chr1_KI270706v1_random' precisely so they stay outside the canonical set — and filtering them in the read path would make their data unreadable rather than merely untidy.
_merge_chromosome_parquet only ever knew the flat layout: it read and wrote variants/<chrom>.parquet and had no branch for the bucketed variants/<chrom>/bucket_N.parquet files. But create-db has produced the bucketed layout since it became the build default, so on every real database the merge found nothing to merge and wrote the incoming samples to a brand-new flat file instead. That file is dead weight. The readers resolve a chromosome's bucket directory first and never fall through to the flat file, so the added samples were never counted as carriers. They still counted toward AN through the capture index, which made them look homozygous reference at every position and pulled allele frequencies down across the whole database. Nothing failed, and afquery check had nothing to report. The read, merge, Phase 2 recompute, sort and write now live in one helper that operates on a single output file, and the chromosome-level function dispatches it either once for a flat file or once per affected bucket. The flat path is therefore exercised by the same code as the bucketed one. filtered_bitmap needs care: it is a function of the whole cohort rather than of which bucket received rows, because enlarging a WES technology changes (tech_bm - carrier_set) at every row of the chromosome. When min_covered is active, every existing bucket of an affected chromosome is revisited, and a per-row comparison keeps buckets whose bitmaps did not move from being rewritten. A chromosome new to the database follows the layout the database already uses, resolved once per update so that all chromosomes in one batch agree.
Eight call sites resolved a chromosome's storage layout for themselves: the query engine's path and glob helpers, the dump and annotate bucket workers, dump's chromosome availability check and work-unit enumeration, the benchmark variant sampler, the compactor, remove-samples, check-database, and the build resume check. Each reimplemented "a bucket directory wins over a flat file", and dump kept its own copy of BUCKET_SIZE. All of them now call the resolver. Behaviour is unchanged: the same precedence, the same per-chromosome decision, the same treatment of an empty chromosome directory. Verified by capturing dump, annotate, check, info, query, batch and region output for a bucketed and a flat database built from the same cohort, and diffing it against the same capture before the change — 22 files, all byte-identical. BUCKET_SIZE now has one definition, re-exported from build and dump so existing imports keep working.
A database updated by the previous add-samples ends up holding both variants/chr1/ and variants/chr1.parquet, and nothing said so. check-database validated each file on its own merits, and each is individually well-formed; the damage is in the pair, because queries read only the bucket directory and every sample whose calls live in the flat file is counted as homozygous reference. check now errors once per affected chromosome and names that consequence, so the failure is legible without opening the Parquet files. add-samples refuses to run on such a database, and does so before ingest starts rather than after an operator has waited out a long VCF pass. remove-samples is deliberately left unguarded. It already clears bits from both layouts, is safe on a mixed database, and is the first step of the recovery procedure — refusing there would leave no way out.
add-samples passed the next free sample id where _update_manifest expects the sample count. The two agree only until a sample is removed, because ids are never reused: after a removal the next id runs ahead of the count. remove-samples already writes SELECT COUNT(*), so the mismatch only appears on a remove-then-add cycle — which is exactly what the recovery procedure for a split variant layout asks operators to perform. Without this, following that procedure ends with check-database warning about a sample_count mismatch the operator caused by doing as they were told.
Every existing add_samples test copies the hand-built fixture from conftest, which writes a flat variants/<chrom>.parquet directly. That is the one layout create-db does not produce, so a merge that understood only the flat layout passed the whole suite while corrupting every real database it touched. These build a database through run_preprocess and assert, through the query engine, that an added sample is counted as a carrier — for a new variant, for a variant the cohort already had, and at bit level in the bucket file, which no test checked after an add before. Also covered: a position past the last existing bucket, a chromosome absent from the database, and that the per chromosome layout is unchanged by the add and never split in two. The coverage-evidence case adds a wes_kit_a sample whose only variant lives in bucket 2 and asserts that bucket 0's filtered_bitmap picked it up, since that bitmap is a function of the cohort rather than of which bucket received rows. A partial fix that only rewrote buckets with new data would pass everything else and fail that one. Eight of the ten fail against the previous merge; the two that do not are guards against a fix that rewrites more than it should.
On chrY the eligible set included females, so they were reported as homozygous reference and counted in n_samples_eligible. A female has no chrY to genotype: she is neither a carrier nor homozygous reference there. AN has always excluded such samples, which left the two disagreeing — in the test cohort a chrY variant reported AN=1 alongside n_eligible=2, and the per-genotype counts summed to more samples than could hold an allele. Fixed where eligibility is decided rather than in the five places that derive N_HOM_REF from it, so point, batch, region, dump, annotate and variant-info all pick it up. split_ploidy already knows which samples carry alleles at a position; intersecting with it leaves autosomes, chrX and chrM untouched, since there both sexes are genotyped and only the ploidy differs. The golden case for chrY recorded the old n_eligible and is updated. Found by the query oracle added in the following commit.
compact selected the rows to keep with table.take(keep_indices) on a plain Python list. When every row of a file is dropped the list is empty, pyarrow infers a null-typed index array, and take has no kernel for (uint32, null) — so the command died with ArrowNotImplementedError partway through, leaving the database half compacted. Emptying a whole file is ordinary: remove a sample that was the only carrier of everything in a 1 Mbp bucket and that bucket has nothing left to keep. Spelling out the index type makes the empty case produce an empty table, which is what the rest of the function already handles.
The suite was green while add-samples silently dropped every added sample, because it only ever asked whether the code agreed with a fixture the same project had written by hand. Nothing compared the numbers against the VCF, BED and manifest text the database was built from. tests/oracle.py works those numbers out from the raw files and imports nothing from afquery — reusing the code under test to predict its own output would only restate its assumptions. The coverage half is the half that matters: AC and the genotype counts come from the VCFs, but AN, N_HOM_REF and N_NO_COVERAGE come from capture-BED intervals crossed with sex-dependent ploidy, which is where every silent counting bug this project has had actually lived. It is checked against a bucketed build, a flat build, a phenotype subset, and a database grown by add-samples. tests/test_invariants.py adds properties that compare one database against another rather than against a recorded expectation: building a cohort in one pass equals building half and adding the rest; the storage layout is not observable through any query; add-then-remove returns the answers to what they were; and compact changes nothing that still has a carrier, twice over. Between them they found the chrY eligibility bug and the compact crash fixed in the two preceding commits.
The flat layout was undocumented, which is part of why a write path that only understood it survived: the data model described bucketed storage as if it were the only kind. It now describes both, says which one create-db produces, and records that a chromosome must never carry two. Troubleshooting gains the entry for databases already split by the previous add-samples: how to recognise it — AN grows, AC does not — and the remove, delete-orphans, re-add sequence that repairs it without a full rebuild. The order matters, so the recovery records the phenotype rows before removal deletes them. The ploidy page already gave chrY a female ploidy of 0; it now says what that means for eligibility, which is that n_eligible and N_HOM_REF count males only there.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #37 +/- ##
==========================================
+ Coverage 87.23% 88.67% +1.43%
==========================================
Files 21 22 +1
Lines 3017 3073 +56
Branches 479 482 +3
==========================================
+ Hits 2632 2725 +93
+ Misses 257 229 -28
+ Partials 128 119 -9 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
annotate.py computed bucket boundaries from a hardcoded 1_000_000 in three places while resolving the file path through storage. If the bucket size ever changed, storage.bucket_path would open one bucket and the WHERE pos BETWEEN range would filter another, so annotation would report AN=0 for most variants with no error. The resolver exists so readers and writers cannot drift; these were the last sites still holding their own copy of the rule.
The recovery procedure deleted every top-level Parquet file with a blanket find -delete, on the assumption that each one has a sibling bucket directory. A broken update-db that introduced a chromosome absent from the original build wrote a flat file with no sibling: check does not flag it because the layout is not mixed, queries read it correctly through the flat fallback, and it is the only copy of that chromosome. Deleting it would drop the chromosome silently. Both the listing and the delete now test for the sibling directory.
add-samples recomputed filtered_bitmap only on the chromosomes the batch carried variants on. But filtered_bitmap is (tech_bm - carriers) row by row, so enlarging a WES tech moves it at every row of every chromosome, not only where the new rows landed. A sample added with variants on chr1 was therefore left out of the coverage bitmap on chr2, and counted as homozygous reference wherever its capture BED reached — silently, and biasing the frequency down. The answer depended on nothing but which chromosome the batch happened to mention. The merge now walks every chromosome the database holds whenever a WES tech grows, and _merge_chromosome_parquet accepts a chromosome with no new rows as a recompute-only pass. The per-file dirty guard already suppresses writes where no bitmap moved, so a WGS-only batch still rewrites nothing. A test asserted the old behaviour was intended; it now covers the case where leaving a chromosome alone is genuinely correct, which is when no coverage threshold is set.
The recompute across every chromosome keyed off wes_tech_bitmaps, which is built from every existing sample that has a BED. In any database holding one capture sample that made a WGS-only batch read and deserialize the entire store, when no bitmap off its own chromosomes could move: filtered_bitmap is (tech_bm - carriers), and on an untouched chromosome carriers cannot change, so only a batch that puts a sample into a capture technology moves anything. Gate on that instead. storage.existing_bucket_ids now treats the chromosome as a directory rather than interpolating it into a glob pattern. It was previously only ever reached with names filtered against the known set; it is now also reached with names read straight off disk, where a contig spelled with an asterisk — GRCh38 does this for HLA — would pull in the buckets of its neighbours.
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
2 participants
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.
update-db --add-sampleshas been writing added samples into a Parquet file the queryengine never opens. On any database built by
create-db, the added samples countedtoward
ANbut never as carriers, so they were reported as homozygous reference at everyposition and pulled allele frequencies down across the whole database. Nothing failed and
afquery checkhad nothing to say.What happened
_merge_chromosome_parquethardcoded the flat layout:It had no branch for the bucketed layout. But
run_preprocesscallsbuild_all_parquetswithout passingpartitioned(preprocess/__init__.py:164), andthat defaults to
True(build.py:520) — so every real database is bucketed. Theflat layout is only produced by unit tests that pass
partitioned=Falseand by thehand-built fixture in
tests/conftest.py.The merge therefore found no
variants/chrN.parquet, started from an empty set, andwrote a brand-new flat file holding only the incoming samples. The bucket files were
never touched.
query.py:315and:326resolve the bucket directory first and neverfall through to the flat file, so that file is dead weight — while eligibility, which
comes from the capture index rather than the variant store, still counted the new
samples.
Observed on a 1225-sample production database after adding two samples:
variants/chr1/bucket_0.parquetat(887801,'A','G')holds 900 hom + 149 het + 55 failwith
max(sample_id) = 1222, and the orphanvariants/chr1.parquetholds exactly{1223, 1224}. Twenty-three orphan files in total, one per chromosome the batch touched.Why no test caught it
The suite was fully green — 533 passing — throughout. Every
add_samplestest copies thefixture from
conftest.py, which writesvariants/<chrom>.parquetby hand. That is theone layout
create-dbdoes not produce, so a merge that understood only the flat layoutpassed everything.
More generally the suite asked whether the code agreed with a fixture the project had
written itself, never whether the numbers agreed with the VCF, BED and manifest text the
database was built from.
Changes
storage.py— one home for the layout rule, which had been reimplemented at eightcall sites across the query engine, the dump and annotate workers, the benchmark
helper, the compactor,
remove_samples,check_databaseand the build resume check.Added first as a standalone module, then the readers moved onto it in a separate
commit. It deliberately does not validate chromosome directory names against the
canonical list: a database may hold unplaced or alt contigs, and filtering them in the
read path would make their data unreadable.
rows by
pos // BUCKET_SIZEand merges each affected bucket.filtered_bitmapis afunction of the whole cohort rather than of which bucket received rows, so when
min_coveredis active every existing bucket of an affected chromosome is revisited; aper-row comparison keeps unchanged buckets from being rewritten.
checkerrors on a chromosome held in both layouts, andadd-samplesrefuses sucha database before ingest rather than after.
remove-samplesstays unguarded onpurpose: it is the first step of the recovery.
sample_countin the manifest was receiving the next free sample id. The two agreeonly until something is removed — which is exactly what the recovery procedure asks for.
genotype;
ANhad always excluded them, so the counts disagreed with it. Fixed whereeligibility is decided, so all six read paths pick it up. Only chrY moves.
compactdied withArrowNotImplementedErrorwhen every row of a file was dropped —remove the only carrier of everything in a bucket and the command failed partway
through, leaving the database half compacted.
Tests
tests/oracle.pyderives the expected counts from the raw VCF, BED and manifest text andimports nothing from afquery. The coverage half is the point:
ACand the genotypecounts come from the VCFs, but
AN,N_HOM_REFandN_NO_COVERAGEcome from captureintervals crossed with sex-dependent ploidy, which is where every silent counting bug in
this project has lived. It runs against a bucketed build, a flat build, a phenotype
subset, and a database grown by
add-samples.tests/test_invariants.pycompares databases against each other rather than againstrecorded numbers: building a cohort in one pass equals building half and adding the rest;
the layout is not observable through any query; add-then-remove restores the answers; and
compact changes nothing that still has a carrier, twice over. The first of those would
have caught this bug on its own.
tests/test_update_partitioned.pycoversadd-samplesagainst a real bucketed database.Eight of its ten tests fail against the previous merge; the other two guard against a fix
that rewrites more than it should.
593 tests pass. The read-path refactor was verified separately by capturing
dump(plain,--all-variants, region,--by-sex --by-tech),annotate,check,info, and point,batch and region queries for a bucketed and a flat database built from the same cohort,
before and after the change: 14 of 22 outputs byte-identical, and the eight that differ do
so on one line each, all of them the intended chrY correction.
Recovering a database already split by this
docs/troubleshooting.mdgains the full procedure. In short: record the affected samples'metadata first — removal deletes their phenotype rows — then
remove-samples, delete theorphan
variants/*.parquet, confirm withcheck, and re-add with the fixed version.Databases only ever built by
create-dbare unaffected.