fix: add-samples wrote to a flat Parquet the query engine never reads

docs/advanced/ploidy-and-sex-chroms.md

@@ -19,6 +19,8 @@ AFQuery computes ploidy-aware AN for sex chromosomes (chrX, chrY) and the mitoch
 
 For each eligible sample at a given position, AFQuery adds the appropriate ploidy count to AN based on the sample's sex and the chromosome/position.
 
+A sample contributing ploidy 0 is not eligible at that position at all: a female has no chrY to genotype, so she is neither a carrier nor homozygous reference there. On chrY, `n_eligible`, `N_HOM_REF` and `AN` therefore count males only.
+
 ---
 
 ## Pseudoautosomal Regions (PAR)

docs/guides/update-database.md

@@ -38,6 +38,8 @@ afquery update-db \
 
 The new manifest follows the same format as the original (see [Manifest Format](manifest-format.md)). New samples are assigned monotonically increasing sample IDs.
 
+New variants are merged into the per-chromosome bucket files the database already uses; buckets are created on demand, and no new top-level Parquet file is produced. See [Data Model](../reference/data-model.md#storage-layouts).
+
 To add multiple manifests at once:
 
 ```bash
@@ -67,9 +69,15 @@ decisions are comparable across batches.
 When new carriers push a partially-covered tech above the `--min-covered`
 threshold at positions that were previously below it, those positions are
 re-evaluated and their non-carrier samples once again count as `N_HOM_REF`
-instead of `N_NO_COVERAGE`. The recomputation runs only for chromosomes
-touched by the new samples; existing rows on other chromosomes are not
-rewritten.
+instead of `N_NO_COVERAGE`. Because that value derives from the tech bitmaps of
+the whole cohort rather than from which file received rows, the recomputation
+spans **every bucket of every chromosome the database holds**, not only the ones
+the new samples carry variants on. Skipping the rest would leave an added WES
+sample counted as `N_HOM_REF` everywhere else its capture BED reaches. Only a
+batch that puts a sample into a capture-based technology can move those bitmaps,
+so a WGS-only batch still visits nothing beyond its own chromosomes; adding a
+WES sample to a large database rewrites broadly and takes minutes rather than
+seconds.
 
 VCFs added via `update-db` should preserve `FORMAT/DP` and `FORMAT/GQ` (the
 bundled `resources/normalize_vcf.sh` does so by default). Samples without

docs/reference/data-model.md

@@ -22,6 +22,18 @@ This page documents the on-disk layout of an AFQuery database, including file fo
     └── wes_v2.pkl
 ```
 
+### Storage layouts
+
+`create-db` always writes the bucketed layout shown above. A single-file-per-chromosome
+layout (`variants/chr1.parquet`, with no bucket directory) is also readable, and appears
+in small hand-built databases and in test fixtures. `update-db` merges into whichever
+layout a chromosome already uses, and creates buckets for a chromosome new to the
+database.
+
+A chromosome must never have both. `afquery check` reports that as an error, because
+queries resolve the bucket directory first and would silently ignore the flat file — and
+every sample stored only in it.
+
 ---
 
 ## manifest.json
@@ -156,6 +168,10 @@ Variants are partitioned into 1-Mbp buckets:
 bucket_id = pos // 1_000_000
 ```
 
+`update-db --add-samples` writes new positions into the bucket that owns
+them, creating `bucket_N.parquet` when a batch extends a chromosome past its
+previous last bucket.
+
 !!! warning "DuckDB integer arithmetic"
     When computing bucket IDs in DuckDB SQL, always use the integer-division operator:
     ```sql

docs/troubleshooting.md

@@ -81,6 +81,84 @@ AFQuery requires DuckDB to use Parquet for all temporary files. Arrow IPC is not
 
 ---
 
+## Samples Added by update-db Are Missing From Queries
+
+**Symptom:** After `afquery update-db --add-samples`, the new samples show up in
+`afquery info --db ./db/ --samples` and `AN` grows by the expected amount, but they
+never appear as carriers: `AC` does not increase at variants you know they carry, and
+`afquery variant-info` does not list them. Allele frequencies across the whole database
+drift downward after every update.
+
+`afquery check --db ./db/` reports one error per affected chromosome:
+
+```
+ERROR  chr1: both variants/chr1/ (bucketed) and variants/chr1.parquet exist.
+```
+
+**Cause:** `update-db --add-samples` before 0.4.1 only understood the
+single-file-per-chromosome variant layout. `create-db` produces the bucketed layout
+(`variants/<chrom>/bucket_N.parquet`), so the merge found nothing to merge and wrote the
+new samples to a fresh `variants/<chrom>.parquet` instead. Queries read the bucketed
+files and ignore that one, so the added samples counted toward `AN` through the capture
+index but never as carriers — they were treated as homozygous reference everywhere.
+
+**Fix:** Upgrade to 0.4.1 or later, then repair the database. `afquery info --db ./db/
+--changelog` lists the samples added by each `add_samples` event.
+
+```bash
+# 0. Stop all writers and back up the small files.
+cp ./db/manifest.json ./db/manifest.json.bak
+cp ./db/metadata.sqlite ./db/metadata.sqlite.bak
+
+# 1. Record the affected samples BEFORE removing them: removal deletes their
+#    phenotype rows, and you need them to rebuild the manifest.
+sqlite3 ./db/metadata.sqlite \
+  "SELECT s.sample_name, s.sex, t.tech_name, s.vcf_path,
+          group_concat(p.phenotype_code)
+   FROM samples s
+   JOIN technologies t ON s.tech_id = t.tech_id
+   LEFT JOIN sample_phenotype p ON p.sample_id = s.sample_id
+   WHERE s.sample_name IN ('SAMPLE_1','SAMPLE_2') GROUP BY s.sample_id;"
+
+# 2. List the flat files before deleting anything. Only those with a sibling
+#    directory of the same name are orphans; one without a sibling is the only
+#    copy of that chromosome and queries still read it.
+for f in ./db/variants/*.parquet; do
+  [ -d "${f%.parquet}" ] && echo "orphan: $f" || echo "KEEP (no sibling): $f"
+done
+
+# 3. Remove the affected samples. This clears their bits from both layouts and is
+#    safe on a split database. Every bucket is read, so budget minutes, not
+#    seconds, and do not interrupt it.
+afquery update-db --db ./db/ --remove-samples SAMPLE_1 --remove-samples SAMPLE_2
+
+# 4. Delete only the orphans listed in step 2. A flat file with no sibling
+#    directory must be kept: deleting it would drop that chromosome entirely.
+for f in ./db/variants/*.parquet; do
+  [ -d "${f%.parquet}" ] && rm "$f"
+done
+
+# 5. Confirm the errors are gone.
+afquery check --db ./db/
+
+# 6. Re-add the samples with the fixed version.
+afquery update-db --db ./db/ --add-samples repair.tsv --bed-dir ./beds/
+
+# 7. Verify a variant you know they carry.
+afquery variant-info --db ./db/ --locus chr1:887801
+```
+
+Re-added samples receive new sample IDs — IDs are never reused after a removal. That is
+expected and does not affect results.
+
+If the original VCFs are no longer available, stop after step 5. The samples are then
+absent from both the metadata and `AN`, which is a correct smaller cohort rather than a
+biased larger one. Rebuilding with `create-db` is always a valid fallback.
+
+Databases that were only ever built with `create-db`, and never updated, are unaffected.
+
+---
+
 ## Compact Takes a Long Time
 
 **Symptom:** `afquery update-db --compact` runs for many minutes or hours.
@@ -116,6 +194,7 @@ afquery info --db ./db/ --samples | grep SAMP
 | `Missing Parquet for chromosome chr3` | Re-run `create-db` or investigate incomplete build |
 | `Manifest mismatch: expected N samples, found M` | Database may be partially updated; re-run `update-db` |
 | `Capture file missing for wes_v1` | BED file was not provided at build time; rebuild with `--bed-dir` |
+| `chr1: both variants/chr1/ ... and variants/chr1.parquet exist` | Samples added by a pre-0.4.1 `update-db` are invisible to queries; see [Samples Added by update-db Are Missing From Queries](#samples-added-by-update-db-are-missing-from-queries) |
 
 ---
 

src/afquery/annotate.py

@@ -3,6 +3,7 @@
 import warnings
 import duckdb
 
+from . import storage
 from .bitmaps import deserialize
 from .constants import normalize_chrom
 from .models import AfqueryWarning, SampleFilter
@@ -48,12 +49,12 @@ def _compute_chunk_annotations(
     n_bitmap_cols = 5 if engine._has_coverage_data else 3
     variant_data: dict[tuple[int, str, str], tuple] = {}
     _db = Path(db_path)
-    bucket_start = bucket_id * 1_000_000
-    bucket_end = (bucket_id + 1) * 1_000_000 - 1
+    bucket_start = bucket_id * storage.BUCKET_SIZE
+    bucket_end = (bucket_id + 1) * storage.BUCKET_SIZE - 1
     cols = ", ".join(engine._bitmap_cols(with_pos=True))
 
-    if chrom in engine._partitioned_chroms:
-        parquet_file = _db / "variants" / chrom / f"bucket_{bucket_id}.parquet"
+    if storage.chrom_layout(_db / "variants", chrom) == storage.PARTITIONED:
+        parquet_file = storage.bucket_path(_db / "variants", chrom, bucket_id)
         if valid_positions and parquet_file.exists():
             con = duckdb.connect()
             placeholders = ", ".join("?" * len(valid_positions))
@@ -67,7 +68,7 @@ def _compute_chunk_annotations(
                 pos, ref, alt = row[0], row[1], row[2]
                 variant_data[(pos, ref, alt)] = tuple(bytes(b) for b in row[3:3 + n_bitmap_cols])
     else:
-        parquet_file = _db / "variants" / f"{chrom}.parquet"
+        parquet_file = storage.flat_path(_db / "variants", chrom)
         if valid_positions and parquet_file.exists():
             con = duckdb.connect()
             rows = con.execute(
@@ -185,7 +186,7 @@ def annotate_vcf(
     for variant in vcf:
 
         norm = normalize_chrom(variant.CHROM)
-        bucket = variant.POS // 1_000_000
+        bucket = storage.bucket_id(variant.POS)
         key = (norm, bucket)
         if key not in variant_buffers:
             work_order.append(key)

src/afquery/benchmark.py

@@ -6,6 +6,7 @@
 
 import pyarrow.parquet as pq
 
+from . import storage
 from .database import Database
 
 
@@ -116,32 +117,18 @@ def _find_test_variants(
         return []
 
     results: list[tuple[str, int, str, str]] = []
-    for entry in sorted(variants_dir.iterdir()):
+    for entry in storage.iter_variant_parquets(variants_dir):
         if len(results) >= n:
             break
-        if entry.suffix == ".parquet":
-            chrom = entry.stem
-            tbl = pq.read_table(str(entry), columns=["pos", "ref", "alt"])
-            for row in range(min(len(tbl), n - len(results))):
-                results.append((
-                    chrom,
-                    int(tbl["pos"][row].as_py()),
-                    str(tbl["ref"][row].as_py()),
-                    str(tbl["alt"][row].as_py()),
-                ))
-        elif entry.is_dir():
-            chrom = entry.name
-            for bucket in sorted(entry.glob("bucket_*.parquet")):
-                if len(results) >= n:
-                    break
-                tbl = pq.read_table(str(bucket), columns=["pos", "ref", "alt"])
-                for row in range(min(len(tbl), n - len(results))):
-                    results.append((
-                        chrom,
-                        int(tbl["pos"][row].as_py()),
-                        str(tbl["ref"][row].as_py()),
-                        str(tbl["alt"][row].as_py()),
-                    ))
+        chrom = entry.stem if entry.parent == variants_dir else entry.parent.name
+        tbl = pq.read_table(str(entry), columns=["pos", "ref", "alt"])
+        for row in range(min(len(tbl), n - len(results))):
+            results.append((
+                chrom,
+                int(tbl["pos"][row].as_py()),
+                str(tbl["ref"][row].as_py()),
+                str(tbl["alt"][row].as_py()),
+            ))
     return results
 
 

src/afquery/dump.py

@@ -7,14 +7,15 @@
 
 import duckdb
 
+from . import storage
 from .bitmaps import deserialize
 from .constants import normalize_chrom, ALL_CHROMS
 from .models import SampleFilter
 from .ploidy import split_ploidy
 
 logger = logging.getLogger(__name__)
 
-BUCKET_SIZE = 1_000_000
+BUCKET_SIZE = storage.BUCKET_SIZE
 
 
 def _build_groups(engine, base_sf, by_sex, by_tech, by_phenotype, all_groups):
@@ -136,8 +137,8 @@ def _dump_bucket_worker(
     bucket_end = (bucket_id + 1) * BUCKET_SIZE - 1
 
     # Resolve parquet path and WHERE clause
-    if chrom in engine._partitioned_chroms:
-        parquet_file = _db / "variants" / chrom / f"bucket_{bucket_id}.parquet"
+    if storage.chrom_layout(_db / "variants", chrom) == storage.PARTITIONED:
+        parquet_file = storage.bucket_path(_db / "variants", chrom, bucket_id)
         if not parquet_file.exists():
             return []
         where_parts = []
@@ -150,7 +151,7 @@ def _dump_bucket_worker(
             params.append(pos_end)
         where_clause = ("WHERE " + " AND ".join(where_parts)) if where_parts else ""
     else:
-        parquet_file = _db / "variants" / f"{chrom}.parquet"
+        parquet_file = storage.flat_path(_db / "variants", chrom)
         if not parquet_file.exists():
             return []
         range_start = max(bucket_start, pos_start) if pos_start is not None else bucket_start
@@ -326,31 +327,23 @@ def dump_database(
         # All chroms that have data
         available = set()
         for chrom in ALL_CHROMS:
-            if chrom in engine._partitioned_chroms:
-                available.add(chrom)
-            elif (variants_dir / f"{chrom}.parquet").exists():
+            if storage.variant_parquet_glob(variants_dir, chrom) is not None:
                 available.add(chrom)
         chroms = [c for c in ALL_CHROMS if c in available]
 
     # Build work units: (chrom, bucket_id) in genomic order
     work_units: list[tuple[str, int]] = []
     for chrom in chroms:
-        if chrom in engine._partitioned_chroms:
-            chrom_dir = variants_dir / chrom
-            bucket_files = sorted(
-                chrom_dir.glob("bucket_*.parquet"),
-                key=lambda p: int(p.stem.split("_")[1]),
-            )
-            for bf in bucket_files:
-                bid = int(bf.stem.split("_")[1])
+        if storage.chrom_layout(variants_dir, chrom) == storage.PARTITIONED:
+            for bid in storage.existing_bucket_ids(variants_dir, chrom):
                 # Filter by region if specified
                 if pos_start is not None and (bid + 1) * BUCKET_SIZE - 1 < pos_start:
                     continue
                 if pos_end is not None and bid * BUCKET_SIZE > pos_end:
                     continue
                 work_units.append((chrom, bid))
         else:
-            flat_path = variants_dir / f"{chrom}.parquet"
+            flat_path = storage.flat_path(variants_dir, chrom)
             if not flat_path.exists():
                 continue
             bucket_ids = _discover_flat_buckets(flat_path, pos_start, pos_end)

src/afquery/preprocess/build.py

@@ -14,12 +14,13 @@
 import pyarrow.parquet as pq
 from pyroaring import BitMap
 
+from .. import storage
 from ..bitmaps import serialize
 from ..constants import ALL_CHROMS
 
 logger = logging.getLogger(__name__)
 
-BUCKET_SIZE = 1_000_000
+BUCKET_SIZE = storage.BUCKET_SIZE
 
 PARQUET_SCHEMA = pa.schema([
     ("pos",                 pa.uint32()),
@@ -685,11 +686,9 @@ def build_all_parquets(
     for chrom in valid_chroms:
         if resume:
             if partitioned:
-                chrom_dir = os.path.join(variants_dir, chrom)
-                done = (os.path.isdir(chrom_dir) and
-                        bool(glob_module.glob(os.path.join(chrom_dir, "bucket_*.parquet"))))
+                done = bool(storage.existing_bucket_ids(variants_dir, chrom))
             else:
-                done = os.path.exists(os.path.join(variants_dir, f"{chrom}.parquet"))
+                done = storage.flat_path(variants_dir, chrom).exists()
             if done:
                 skipped_chroms.append(chrom)
                 continue

src/afquery/preprocess/compact.py

@@ -1,4 +1,3 @@
-import glob as glob_module
 import json
 import logging
 import os
@@ -11,6 +10,7 @@
 import pyarrow.parquet as pq
 from pyroaring import BitMap
 
+from .. import storage
 from ..bitmaps import deserialize, serialize
 from .build import PARQUET_SCHEMA
 
@@ -43,13 +43,7 @@ def compact_database(db_path: Path) -> dict:
     active_ids = BitMap([r[0] for r in rows])
 
     # Collect all parquet files (flat + partitioned buckets)
-    all_parquets: list[Path] = []
-    for f in sorted(variants_dir.glob("*.parquet")):
-        all_parquets.append(f)
-    for chrom_dir in sorted(variants_dir.iterdir()):
-        if chrom_dir.is_dir():
-            for f in sorted(chrom_dir.glob("bucket_*.parquet")):
-                all_parquets.append(f)
+    all_parquets: list[Path] = list(storage.iter_variant_parquets(variants_dir))
 
     logger.info("[compact] Compacting %d Parquet file(s) against %d active sample(s)...",
                 len(all_parquets), len(active_ids))
@@ -115,8 +109,11 @@ def compact_database(db_path: Path) -> dict:
             logger.debug("  [compact] %s: no changes", parquet_file.name)
             continue
 
-        # Build new table with kept rows and updated bitmaps
-        orig_keep = table.take(keep_indices)
+        # Build new table with kept rows and updated bitmaps.
+        # The index type is spelled out: a bare [] makes pyarrow infer a null
+        # array, which has no take kernel, and every row of a file can legitimately
+        # be dropped when the removed samples were the only carriers in it.
+        orig_keep = table.take(pa.array(keep_indices, type=pa.int64()))
         new_table = pa.table(
             {
                 "pos":                 orig_keep["pos"],