diff --git a/docs/pipelines/association.md b/docs/pipelines/association.md index f055e45..69ee4b6 100644 --- a/docs/pipelines/association.md +++ b/docs/pipelines/association.md @@ -342,3 +342,12 @@ The point of the publication stage is not just file creation. It is the point wh ## Current additive metadata path The association export manifest framework adds reserved `_datahub` metadata blocks to published outputs. These blocks are additive and allow the project to carry forward provenance and future coverage fields without breaking existing payload consumers. + +## Association evidence v2 sidecar + +The release-candidate v2 model separates exact provider rows, source association +observations, canonical variants, consequence annotations, clinical assertions, +population observations, and variant-phenotype summaries. It replaces the +legacy representative-row annotation contract only when the backend is explicitly +configured for v2. See [Association Evidence v2](../schemas/association_evidence_v2.md) +for counting, source completeness, build, migration, and rollback semantics. diff --git a/docs/reference/scripts.md b/docs/reference/scripts.md index 7c68c9e..d0351be 100644 --- a/docs/reference/scripts.md +++ b/docs/reference/scripts.md @@ -620,3 +620,39 @@ DuckDB before the current points-table ingest/publish stages run. ## Script philosophy The scripts directory is intentionally operational. Business/scientific logic lives in `src/datahub/` when it can. Scripts compose that logic and add environment/runtime concerns such as CLI parsing, checkpoints, and scheduler integration. + +## Association evidence v2 + +### `scripts/dataset_specific_scripts/unified/build_association_evidence_v2.py` + +Builds the resumable non-composite association sidecar from archived CVD/trait +provider rows and registers retained variant-index artifacts for gene-scoped runtime +access when provider detail is unavailable. It writes only to explicit +output, checkpoint, manifest, temporary, and log paths. Use `--dry-run` to inspect +the plan and `--checksum-inputs` only when a second full read of source files is +acceptable. `--provider-chunk-rows 0` uses direct inserts for inputs up to 2 GB +and bounded one-million-row chunks for oversized inputs; set a positive value to +force bounded chunks for every provider file. + +### `scripts/dataset_specific_scripts/unified/build_association_evidence_v2_serving.py` + +Publishes the checksum-linked, bounded Parquet serving package from a verified normalized sidecar. The production release candidate uses 4,096 variant buckets for exact variant evidence and gene-keyed projections for summary requests. + +### `scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_index.py` + +Streams every registered MVP compact artifact into a resumable gene-keyed Parquet association index. It emits deterministic association-record IDs, record_kind=source_summary_association, evidence_granularity=source_summary, exact phenotype paths, reported p-values, retained compact fields, and immutable artifact provenance. It never fabricates provider or study records. + +### build_association_evidence_v2_source_summary_rollup.py + +Builds resumable gene-keyed default-summary and exact-phenotype projections from +the completed source-summary index. It records exact phenotype variant-ID sets so +provider-level legacy and source-summary MVP evidence can be unioned without duplicate +inflation. It advances the serving manifest to schema **2.10.0-rc2**; it does not +replace the source-summary index or create inferred provider records. + +### `scripts/dataset_specific_scripts/unified/audit_association_evidence_v2_release.py` + +Reproduces portal-wide v1/v2 count and category comparisons using the normalized sidecar plus distinct compact-source membership, with a resumable runtime checkpoint. + +See the [v2 schema](../schemas/association_evidence_v2.md) for the full +command and scientific contract. diff --git a/docs/runbooks/association-evidence-v2-release.md b/docs/runbooks/association-evidence-v2-release.md new file mode 100644 index 0000000..354a499 --- /dev/null +++ b/docs/runbooks/association-evidence-v2-release.md @@ -0,0 +1,82 @@ +# Association Evidence v2 Release Candidate + +Use this runbook only after the staged sidecar, backend, and frontend have passed +scientific, reconciliation, performance, and browser gates. + +## Pre-deployment gates + +- DataHub manifest and output SHA-256 verify. +- Every DataHub QA check passes. +- The manifest records the exact `--provider-chunk-rows` policy used for oversized sources. +- Source-completeness statuses match the archived inputs. +- The complete gene-keyed MVP source-summary association index, file manifest, and checksums verify; compact rows are first-class association evidence and no provider or study rows are fabricated. +- The schema 2.10 default-summary rollups, exact phenotype variant-ID sets, file manifest, and checksums verify against the complete source-summary index. +- Five bounded multiplicity fixtures reconcile provider rows, source + observations, annotations, assertions, summaries, and exports. +- All v1/v2 category differences are classified and reviewed. +- The full comparison audit reuses the verified schema 2.10 distinct-membership rollup; it does not rescan compact source artifacts. +- Summary, detail, selected-rsID drill-down, and ZIP export use the same filter scope. +- Default search and explicit p-value searches pass performance limits. +- Browser tests pass for HMGCR and at least one high-volume gene. +- A timestamped environment/configuration backup and exact v1 rollback command exist. + +## Deployment shape + +Evidence v2 is an immutable sidecar selected by environment: + +```text +HBP_ASSOC_EVIDENCE_MODEL=v2 +HBP_ASSOC_EVIDENCE_V2_PATH=/data/DataHub/releases/association-evidence-v2/RELEASE_ID/association_evidence_v2.duckdb +HBP_ASSOC_EVIDENCE_V2_SERVING_MANIFEST=/data/DataHub/releases/association-evidence-v2/RELEASE_ID/serving/serving-manifest.json +HBP_ASSOC_ARTIFACT_ROOT=/data/DataHub/analyzed_data/association_new/final +``` + +Do not overwrite `association_serving_slim.duckdb`. Keep the prior backend +environment and v1 artifacts available until post-deployment validation completes. +The backend opens the normalized sidecar and serving package read-only. It resolves first-class MVP source-summary associations from the staged gene-keyed index and retains `HBP_ASSOC_ARTIFACT_ROOT` only as a pre-release fallback when that index is absent. Production v2 validation requires the complete index; request-time parsing of large compact JSON artifacts is not an accepted deployment state. The backend retains an explicit v1 fallback. MVP records use evidence_granularity=source_summary and provider_detail_status=not_applicable; provider and study rows are never fabricated. + +Build the release-candidate serving tree with +`--variant-bucket-characters 3` and a checksum-verified +`--coarse-serving-root`. Provider and association records then use +4,096 bounded SHA-256 prefix buckets; the table-specific width is recorded in +the serving manifest and read by the backend. This changes only physical query +partitioning, not normalized evidence, counts, filters, or provenance. + +After the complete source-summary index verifies, run **build_association_evidence_v2_source_summary_rollup.py**. Production validation requires **schema_version=2.10.0-rc2**. The backend may use these rollups only for an unfiltered default summary and exact CVD phenotype counting. Any phenotype-path or explicit p-value filter, drill-down, export, or provenance request continues to use the full indexed rows. + +## Approved migration sequence + +1. Verify the staged SHA-256 manifest and DataHub QA report without opening the file + for writing. +2. Copy the immutable sidecar, complete serving partition tree, and serving manifest + into a new release-ID directory; never reuse an existing release path. Preserve or + dereference hard links deliberately during the copy and verify every copied + checksum. +3. Save a timestamped copy of the backend service environment and current v1/v2 + selector values. Record the deployed backend/frontend commits and build checksum. +4. Set the two environment variables above and restart only the backend during an + approved maintenance window. No DuckDB import or schema migration is performed. +5. Require exact HMGCR, PCSK9, LDLR, TTN, and ANK2 summary/detail/export + reconciliation, source-completeness status, API latency, browser, CAPTCHA export, + and log-health gates before declaring the release live. +6. Preserve the staged release, copied release, v1 serving DB, previous environment, + and all validation reports for the rollback window. + +## Rollback + +Restore the timestamped prior environment or set: + +```text +HBP_ASSOC_EVIDENCE_MODEL=v1 +``` + +Then restart only the backend in an approved window and repeat the recorded v1 +health/search/download checks. The v2 sidecar can remain on disk because v1 mode does +not open it. No DataHub artifact deletion, reverse transformation, or database +conversion is required. If rollback validation fails, keep the service stopped and +restore the complete timestamped service environment rather than editing data. + +The application default remains v1. A deployment enters v2 only when the selector, +normalized sidecar, and serving manifest are all explicitly configured. Missing or +invalid serving partitions fail the v2 validation gate; they are not replaced by a +request-time global scan. diff --git a/docs/schemas/association_evidence_v2.md b/docs/schemas/association_evidence_v2.md new file mode 100644 index 0000000..e51b46b --- /dev/null +++ b/docs/schemas/association_evidence_v2.md @@ -0,0 +1,286 @@ +# Association Evidence v2 + +Association evidence v2 separates source observations from variant annotations and +visualization summaries. It is an additive DuckDB sidecar; it does not replace the +legacy association serving database during release-candidate validation. + +## Scientific entity model + +| Table | Unit | Contract | +| --- | --- | --- | +| `variants` | Canonical variant identifier | Stores unordered source allele and coordinate contexts. REF, ALT, and effect allele remain unresolved unless an authoritative association source supplies them. | +| `provider_records` | Exact archived input row | Stable ID from source, logical file, line number, and row checksum. Raw source columns are preserved without cross-row field filling. | +| `association_records` | Source association evidence | Recoverable legacy rows use `record_kind=source_association_observation` and link to exact provider rows. MVP compact rows use `record_kind=source_summary_association`, `evidence_granularity=source_summary`, source `million_veteran_program`, display name `MVP`, and `provider_detail_status=not_applicable`. Neither contract fills absent fields across records. | +| `source_summary_artifacts` | Retained MVP gene/dataset artifact | Registers the immutable compact artifact supplying MVP source-summary associations. Provider-row lineage is not applicable, but every retained variant-phenotype summary is valid association evidence. The gene-keyed serving index preserves deterministic record IDs, exact phenotype paths, p-values, variation type and other retained compact fields without creating study/provider records. | +| `consequence_annotations` | Source/version-scoped annotation | Preserves each consequence/transcript/source combination. The version field is explicit and may be unavailable for legacy snapshots. No lexical minimum or implicit severity selection is performed. | +| `clinical_assertions` | Variant + exact normalized term + assertion source | Preserves all 15 audited legacy ClinVar-derived terms without selecting a representative classification. Each assertion links to every supporting provider row; raw values remain unchanged in provider provenance. Condition, accession, review status, classification date, and ClinVar release are unavailable. | +| `population_observations` | Source/cohort/allele/build observation | Keeps recoverable association-source frequency fields separate. The API joins the full dbSNP frequency sidecar by rsID and labels the association allele unresolved. | +| `variant_phenotype_summaries` | Dataset type + gene + exact phenotype path + variant ID | Derived visualization summary with source/association counts and minimum reported p-value plus the supplying record ID. | + +## Counting and filters + +Variant Annotation Composition counts distinct `variant_id` values under the +active CVD/trait phenotype paths and an explicitly selected p-value threshold. +Variation type contributes one derived or unresolved class per variant. A variant +may contribute to more than one consequence or clinical-significance term because +all source-linked annotations are retained. Clinical display groups are nonexclusive: +a variant can contribute to more than one exact term and more than one group. The +legacy annotations may not refer to the selected HBP phenotype and are not interpreted +as condition-specific ClinVar assertions or clinical actionability. + +Clinical values are parsed from single terms, Python/JSON lists, unquoted bracketed +lists, comma- or semicolon-delimited values, and audited spelling variants. Recognized +slash classifications remain atomic. Assertion identity is deterministic over +`variant_id`, normalized exact term, and assertion source. Repeated identical terms in +one provider value are deduplicated only in the normalized relationship; the original +raw value remains in `provider_records` and exports. + +The same serialized association-filter object is used for summary, detail, and +association export. Provider association rows pass an explicit p-value filter only +when their reported p-value is present and below the threshold. Variant annotations +are version-scoped entities and remain attached to a selected variant rather than being +recast as association measurements. If their exact provider row does not pass the +association p-value filter, API and export bundles retain it only as +`annotation_provenance_only`; it is never presented as a matching association. All +provider-link tables in a bundle resolve to an included provider row. + +## Summary count fields + +- `source_observation_count` counts exact recoverable provider rows. +- `association_record_count` counts both provider-level legacy associations and first-class MVP source-summary associations. +- `provider_association_record_count` counts only associations with provider-row lineage. +- `source_summary_association_count` counts MVP summary associations; `retained_source_summary_count` remains a compatibility alias during migration. +- `minimum_reported_p_value` is an aggregate; + `minimum_reported_p_value_association_record_id` and the tied-record ID list identify + its exact supplier. + +## MVP evidence contract + +MVP is one source: machine identifier million_veteran_program and display label MVP. +Each natural key consists of source, dataset type, gene, variant ID, and exact +phenotype path, plus retained artifact context. MVP and legacy records remain separate +even when they share a variant and phenotype. Distinct-variant charts may count the +variant once within an active scope, but drill-downs and exports retain every +source-specific record. Missing effect alleles, effects, standard errors, sample +sizes, study ancestry, fine-mapping values, study IDs, and provider rows are recorded +as not provided in the current MVP summary; they do not make MVP evidence unavailable. +Source-summary consequence, clinical-significance, disease/trait, source-context, +ancestry, and metadata fields are exposed explicitly. Any unforeseen compact fields +remain in a small extras object, so mapped values are not duplicated in every row. + +## Source completeness + +The archived AWS CVD and trait files are GWAS Catalog-derived and preserve exact +provider rows, but their exact upstream snapshot release is unavailable. They do not +contain effect allele, effect size, standard error, sample size, association ancestry, +or fine-mapping fields. Those values and statuses are explicitly unavailable. + +Provider-level MVP inputs are not present in the AWS snapshot, but the retained +compact MVP variant-phenotype rows are first-class source-summary association evidence. +They use `record_kind=source_summary_association`, +`evidence_granularity=source_summary`, and +`provider_detail_status=not_applicable`. We do not fabricate provider, study, allele, +effect, sample-size, ancestry, or fine-mapping fields that are absent from the retained +summary. + +Every public MVP row identifies `dataset_id=hbp_mvp_association`, preserves the exact +phenotype path, exposes its leaf as `phenotype_slug`, and keeps the source +classification in `phenotype_kind` (`disease` or `trait`). Public API and export +records retain the logical artifact identifier, such as +`variant_index/CVD/PCSK9.json.gz`, but omit internal host and scratch paths. Raw source +path metadata remains confined to controlled build inputs and is not part of the public +serving contract. + +MVP rows participate in primary association responses, variant-phenotype summaries, +source filters, drill-downs, and exports. MVP and legacy records remain separate when +they share a variant and phenotype, while distinct-variant charts count the variant once +inside the active scope. Canonical consequence and clinical entities may be linked as +annotation provenance when independently available; they do not become fabricated MVP +provider records. + +## Build + +The production-scale release candidate is built with: + +```bash +PYTHONPATH=src .venv/bin/python \ + scripts/dataset_specific_scripts/unified/build_association_evidence_v2.py \ + --cvd-root /data/DataHub/raw_data/cvd \ + --trait-root /data/DataHub/raw_data/trait \ + --variant-index-root /data/DataHub/analyzed_data/association_new/final/variant_index \ + --output-db /path/to/staging/association_evidence_v2.duckdb \ + --checkpoint-path /path/to/staging/checkpoint.json \ + --manifest-path /path/to/staging/manifest.json \ + --release-id RELEASE_ID \ + --threads 4 \ + --memory-limit 8GB \ + --provider-chunk-rows 0 \ + --temp-directory /path/to/staging/duckdb_tmp \ + --progress-interval 1 \ + --log-path /path/to/staging/build.log +``` + +The command is resumable by source file and derived phase. With +`--provider-chunk-rows 0`, inputs up to 2 GB use DuckDB's direct single-pass insert; +larger inputs are streamed once through CSV-aware one-million-row temporary chunks. +A positive value forces that bounded row count for every provider file. Temporary +chunks are deleted immediately after their transaction commits, and provider IDs +remain based on the original logical path, record ordinal, and canonical row hash. + +Derived publication is also bounded and resumable. Provider-to-association +relationships are partitioned by stable association-record hash; variants, +consequences, clinical assertions, and recoverable population observations are +partitioned by stable variant-ID hash. Summary inputs and annotation joins use the +same variant buckets. Each bucket is one atomic transaction and is recorded in a +phase-specific checkpoint, so an interrupted run resumes only uncommitted buckets +without rescanning or deleting completed buckets. MVP compact artifacts are first registered from immutable file metadata without creating provider records. A separate resumable serving step streams each registered artifact and writes only deduplicated source-summary IDs, variant IDs, phenotype paths, representative compact p-values, completeness fields, and artifact provenance into gene-keyed Parquet. Primary-key and serving indexes +are deferred until bulk publication completes; deterministic ID uniqueness is a +release QA gate. + +The builder writes a standard sibling `.duckdb.sha256` file from the same full-file digest recorded in the JSON manifest. The serving publisher records this digest so normalized and serving artifacts share one immutable source identity. + +The manifest records the +Git base, dirty status, transformation-file checksums, phenotype-tree checksum, +input sizes/checksum status, runtime versions, counts, output checksum, and exact +command. Input checksums are optional because requesting them performs a second +full read of large source files. + +## Bounded serving package + +The backend does not use global joins across the production-scale sidecar during +requests. After the normalized database passes build QA, create a separate immutable +serving package: + +```bash +PYTHONPATH=src .venv/bin/python \\ + scripts/dataset_specific_scripts/unified/build_association_evidence_v2_serving.py \\ + --source-db /path/to/staging/association_evidence_v2.duckdb \\ + --intermediate-root /path/to/staging/duckdb_tmp \\ + --output-root /path/to/staging/association_evidence_v2_serving \\ + --coarse-serving-root /path/to/staging/association_evidence_v2_serving_256 \\ + --threads 1 \\ + --memory-limit 8GB \\ + --variant-bucket-characters 3 \\ + --verbose +``` + +Then build the first-class MVP source-summary association index from the registered immutable artifacts: + +```bash +PYTHONPATH=src .venv/bin/python \ + scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_index.py \ + --source-db /path/to/staging/association_evidence_v2.duckdb \ + --association-artifact-root /data/DataHub/analyzed_data/association_new/final \ + --serving-root /path/to/staging/association_evidence_v2_serving \ + --workers 3 \ + --progress-interval 1000 \ + --verbose +``` + +This step is resumable per registered artifact. It validates input size and modification time, preserves one stable source-summary ID per source/dataset/gene/variant/phenotype/artifact identity, and records selected and deduplicated row counts plus per-file SHA-256 checksums. It does not infer studies, alleles, sample sizes, effects, ancestries, or fine-mapping fields. + +After the complete 2.7 source-summary index verifies, build the default-search rollups: + +~~~bash +PYTHONPATH=src .venv/bin/python \ + scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_rollup.py \ + --serving-root /path/to/staging/association_evidence_v2_serving \ + --threads 2 \ + --memory-limit 4GB \ + --progress-interval 1 \ + --verbose +~~~ + +The rollup advances the serving schema to **2.10.0-rc2** and writes two immutable, +gene-keyed physical projections. **source_summary_association_base_by_gene** +contains one row per gene, dataset type, and variant for default chart summaries. +**source_summary_association_phenotype_counts_by_gene** contains one row per gene, +dataset type, and exact phenotype path, including the exact sorted variant-ID set. +Those sets permit a deduplicated union with recoverable provider evidence; counts +must never be added across the two sources. The rollups are query accelerators only. +The complete source-summary index remains authoritative for phenotype and explicit +p-value filters, row drill-downs, exports, and source provenance. + + +The serving manifest records table paths, partition keys, bucket functions, +row/file counts, runtime versions, and known limitations. Request-time projections +include `summary_base_by_gene`, `consequence_annotations_by_gene`, and +variant-keyed phenotype summaries, association records, provider records, +clinical assertions, population observations, and provenance partitions. + +The production release candidate uses a three-character SHA-256 hexadecimal +prefix for provider, association, phenotype-summary, consequence, clinical, and +population records, yielding 4,096 variant buckets. +The publisher reads each verified two-character coarse partition once, writes +its 16 fine partitions under one temporary unit, and atomically checkpoints the +completed coarse unit. It does not rescan the normalized database for every +fine bucket. An interrupted run removes only its incomplete temporary unit and +resumes after the last committed coarse bucket. +Each table declares its own prefix width in the manifest; gene-scoped and reused +annotation projections retain their declared two-character prefix. The backend +must read the table-specific value rather than assume one global bucket width. + +Summary requests read one gene bucket, apply exact phenotype paths and an explicit +p-value threshold, and count distinct variant IDs. Detail and export requests read +only selected variant buckets. The backend must be configured with both the +normalized database and serving manifest. + +Provider rows lacking `variant_id_raw` remain preserved in the normalized sidecar but are excluded from variant-keyed serving partitions. The serving manifest records the normalized total, excluded count, and exclusion reason; no unkeyed row is silently presented as variant evidence. + + +## Full runtime count audit + +The portal-wide comparison must include MVP source-summary association membership +without converting compact entries into association records. Run the resumable audit +in staging: + +~~~bash +PYTHONPATH=src .venv/bin/python \ + scripts/dataset_specific_scripts/unified/audit_association_evidence_v2_release.py \ + --v2-db /path/to/staging/association_evidence_v2.duckdb \ + --v1-overall-root /data/DataHub/analyzed_data/association_new/final/overall \ + --v1-variant-index-root /data/DataHub/analyzed_data/association_new/final/variant_index \ + --artifact-root /data/DataHub/analyzed_data/association_new/final \ + --source-summary-serving-root /path/to/staging/association_evidence_v2_serving \ + --output-dir /path/to/staging/release-audit \ + --progress-interval 1000 \ + --workers 3 +~~~ + +When a completed schema 2.8 serving root is supplied, the audit verifies and +reuses its exact gene/dataset/variant rollup in 256 resumable gene buckets. A compatible +schema 2.7 package falls back to the complete gene-keyed source-summary index. Without +either serving package, the audit streams every registered compact JSON array once. +All modes store only distinct `(gene, dataset_type, variant_id)` membership in a +resumable staging DuckDB, and +combines that membership with recoverable provider summaries. Variation type comes +only from normalized provider allele contexts; consequence and clinical categories +come only from source/versioned provider annotations. Compact representative +annotation fields are never counted as v2 annotations. + +## Migration and rollback + +1. Build and validate the sidecar outside `datamart/`. +2. Compare every v1/v2 chart category and document expected semantic differences. +3. Test the backend with `HBP_ASSOC_EVIDENCE_V2_PATH` and + `HBP_ASSOC_EVIDENCE_V2_SERVING_MANIFEST` pointing to staged artifacts. +4. Copy the verified immutable sidecar and complete serving package to one versioned + release path. +5. Set `HBP_ASSOC_EVIDENCE_MODEL=v2`, the versioned sidecar path, and versioned + serving-manifest path during an approved deployment. +6. Retain the v1 serving DB and release manifest. + +Rollback does not transform data. Set `HBP_ASSOC_EVIDENCE_MODEL=v1` (or restore +the prior service environment) and restart only during an approved deployment +window. The legacy serving DB remains intact. + +## Known limitations + +- Association REF/ALT/effect-allele identity is unresolved in current AWS inputs. +- Full dbSNP population observations are joined by rsID only; frequency-source + REF/ALT/build values are not confirmed as the tested association allele. +- Legacy SnpEff- and ClinVar-derived fields do not retain authoritative upstream + release versions in the archived rows. +- MVP provider-level study rows cannot be reconstructed from compact summaries. +- This sidecar does not claim that the entire production portal can be rebuilt + from public inputs with one command. diff --git a/mkdocs.yml b/mkdocs.yml index 5dec23f..5014d62 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -61,6 +61,7 @@ nav: - Runbooks: - Local Smoke Test: runbooks/local-smoke-test.md - Release Checklist: runbooks/release-checklist.md + - Association Evidence v2 RC: runbooks/association-evidence-v2-release.md - Reference: - Data Sources: reference/data-sources.md - Module Guide: reference/module-guide.md @@ -68,6 +69,7 @@ nav: - Config Guide: reference/config-guide.md - Schemas: - Overview: schemas/index.md + - Association Evidence v2: schemas/association_evidence_v2.md - Expression v3: schemas/expression_v3.md - Gene Profile: schemas/gene_profile.md - Protein Context: schemas/protein_context.md diff --git a/scripts/dataset_specific_scripts/unified/audit_association_evidence_v2_release.py b/scripts/dataset_specific_scripts/unified/audit_association_evidence_v2_release.py new file mode 100755 index 0000000..c36aaf9 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/audit_association_evidence_v2_release.py @@ -0,0 +1,1233 @@ +#!/usr/bin/env python3 +"""Compare association evidence v2 category counts with published v1 artifacts.""" + +from __future__ import annotations + +import argparse +import concurrent.futures +import csv +import gzip +import hashlib +import json +import logging +import os +import time +from collections import defaultdict +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import duckdb + + +EXAMPLES = ( + ("TTN", "TRAIT", "qrs_interval_europeans", "rs202094100"), + ("PCSK9", "TRAIT", "qrs_interval_europeans", "rs141502002"), + ("ANK2", "CVD", "mitral_annular_calcification", "rs10013743"), + ("HMGCR", "TRAIT", "rr_interval", "rs2303151"), + ("BMPR2", "CVD", "mitral_annular_calcification", "rs1048829"), +) +DIMENSIONS = ("vc", "msc", "cs") +REASONS = { + "vc": ( + "v2 derives a canonical class only from parseable unordered source alleles; " + "provider-unavailable variants remain unresolved" + ), + "msc": ( + "v2 preserves every source/transcript consequence and counts distinct variants " + "per term; v1 selected one compact representative term" + ), + "cs": ( + "v2 preserves every source clinical assertion and counts distinct variants per " + "term; v1 selected one compact representative term" + ), +} + + +def category_difference_reason( + *, + gene: str, + dataset_type: str, + dimension: str, + v1_count: int, + v2_count: int, + v2_totals: dict[tuple[str, str], dict[str, Any]], +) -> str: + if v1_count == v2_count: + return "" + if (gene, dataset_type) not in v2_totals: + return ( + "No v2 gene/dataset summary was recoverable. The v1 identifier may be invalid, " + "composite, unsupported by the archived provider gene field, or absent from the " + "retained source-summary artifacts; review this gene row explicitly." + ) + if v1_count == 0 and v2_count > 0: + return ( + f"{REASONS[dimension]}; this term was not exposed by the v1 representative " + "entry but is present in separated v2 evidence." + ) + if v2_count == 0 and v1_count > 0: + return ( + f"{REASONS[dimension]}; this v1 representative category is not supported by " + "the v2 derivation or retained source-linked annotation set." + ) + return REASONS[dimension] + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--v2-db", required=True, type=Path) + parser.add_argument("--v1-overall-root", required=True, type=Path) + parser.add_argument("--v1-variant-index-root", required=True, type=Path) + parser.add_argument("--artifact-root", type=Path, help="Root used to resolve logical source-summary artifacts; defaults to the variant-index parent.") + parser.add_argument( + "--source-summary-serving-root", + type=Path, + help=( + "Completed association-evidence-v2 serving root containing the " + "gene-keyed unavailable-provider source-summary index. When supplied, " + "the audit reuses its verified Parquet partitions instead of rereading " + "the registered compact JSON artifacts." + ), + ) + parser.add_argument("--output-dir", required=True, type=Path) + parser.add_argument("--v2-sha256", help="Verified manifest checksum; avoids a second full database read.") + parser.add_argument("--checkpoint-path", type=Path) + parser.add_argument("--log-path", type=Path) + parser.add_argument("--progress-interval", type=int, default=1000) + parser.add_argument( + "--workers", + type=int, + default=1, + help="Parallel immutable artifact readers; DuckDB writes remain single-process.", + ) + parser.add_argument("--max-v1-files", type=int, help="Smoke-test limit per dataset type.") + parser.add_argument("--reset", action="store_true", help="Reset only this audit's checkpoint and intermediate files.") + return parser.parse_args() + + +def configure_logging(path: Path) -> logging.Logger: + path.parent.mkdir(parents=True, exist_ok=True) + logger = logging.getLogger("datahub.association_evidence_v2.audit") + logger.handlers.clear() + logger.setLevel(logging.INFO) + formatter = logging.Formatter("%(asctime)s | %(levelname)s | %(name)s | %(message)s") + stream = logging.StreamHandler() + stream.setFormatter(formatter) + logger.addHandler(stream) + file_handler = logging.FileHandler(path) + file_handler.setFormatter(formatter) + logger.addHandler(file_handler) + return logger + + +def configuration_hash(args: argparse.Namespace) -> str: + value = { + "v2_db": str(args.v2_db.resolve()), + "v1_overall_root": str(args.v1_overall_root.resolve()), + "v1_variant_index_root": str(args.v1_variant_index_root.resolve()), + "artifact_root": str((args.artifact_root or args.v1_variant_index_root.parent).resolve()), + "source_summary_serving_root": ( + str(args.source_summary_serving_root.resolve()) + if args.source_summary_serving_root else None + ), + "runtime_compact_membership": str( + (args.output_dir / "runtime-compact-membership.duckdb").resolve() + ), + "max_v1_files": args.max_v1_files, + } + return hashlib.sha256(json.dumps(value, sort_keys=True).encode()).hexdigest() + + +def save_checkpoint(path: Path, configuration: str, phase: str, details: dict[str, Any]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text(json.dumps({ + "version": 1, + "configuration_hash": configuration, + "phase": phase, + "updated_at": datetime.now(timezone.utc).isoformat(), + "details": details, + }, indent=2, sort_keys=True) + "\n") + temporary.replace(path) + + +def load_checkpoint(path: Path, configuration: str) -> dict[str, Any]: + if not path.exists(): + return {} + value = json.loads(path.read_text()) + if value.get("configuration_hash") != configuration: + raise RuntimeError("Audit checkpoint configuration does not match this run; use --reset.") + return value + + +def encode_count_bundle(category: dict, totals: dict) -> dict[str, Any]: + return { + "category": [ + {"gene": key[0], "dataset_type": key[1], "dimension": key[2], "label": key[3], "count": count} + for key, count in sorted(category.items()) + ], + "totals": [ + {"gene": key[0], "dataset_type": key[1], **value} + if isinstance(value, dict) + else {"gene": key[0], "dataset_type": key[1], "count": value} + for key, value in sorted(totals.items()) + ], + } + + +def decode_count_bundle(value: dict[str, Any], *, rich_totals: bool) -> tuple[dict, dict]: + category = { + (row["gene"], row["dataset_type"], row["dimension"], row["label"]): int(row["count"]) + for row in value.get("category") or [] + } + if rich_totals: + totals = { + (row["gene"], row["dataset_type"]): { + key: int(item) for key, item in row.items() + if key not in {"gene", "dataset_type"} + } + for row in value.get("totals") or [] + } + else: + totals = { + (row["gene"], row["dataset_type"]): int(row["count"]) + for row in value.get("totals") or [] + } + return category, totals + + +def normalized_label(value: Any) -> str: + text = str(value or "").strip().replace("_", " ").casefold() + text = " ".join(text.split()) + if text in {"snp", "single nucleotide variant"}: + return "snv" + return text or "unannotated" + + +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def write_csv(path: Path, rows: list[dict[str, Any]]) -> None: + fieldnames = sorted({key for row in rows for key in row}) + with path.open("w", newline="", encoding="utf-8") as stream: + writer = csv.DictWriter(stream, fieldnames=fieldnames) + writer.writeheader() + writer.writerows(rows) + + +def iter_json_array(path: Path, *, chunk_size: int = 1024 * 1024): + """Stream objects from a JSON array without loading a gene artifact in memory.""" + + opener = gzip.open if path.name.endswith(".gz") else open + decoder = json.JSONDecoder() + with opener(path, "rt", encoding="utf-8") as stream: + buffer = "" + position = 0 + started = False + eof = False + while True: + if not eof: + chunk = stream.read(chunk_size) + if chunk: + buffer += chunk + else: + eof = True + while True: + while position < len(buffer) and buffer[position].isspace(): + position += 1 + if not started: + if position >= len(buffer): + break + if buffer[position] != "[": + raise ValueError(f"Expected JSON array in {path}") + started = True + position += 1 + continue + while position < len(buffer) and ( + buffer[position].isspace() or buffer[position] == "," + ): + position += 1 + if position < len(buffer) and buffer[position] == "]": + return + if position >= len(buffer): + break + try: + value, end = decoder.raw_decode(buffer, position) + except json.JSONDecodeError: + if eof: + raise + break + if not isinstance(value, dict): + raise ValueError(f"Expected object entries in {path}") + yield value + position = end + if position: + buffer = buffer[position:] + position = 0 + if eof: + if buffer.strip(): + raise ValueError(f"Truncated JSON array in {path}") + return + + +def _flush_membership_rows( + connection: duckdb.DuckDBPyConnection, + rows: list[tuple[str, str, str]], + batch_csv: Path, +) -> None: + if not rows: + return + with batch_csv.open("w", newline="", encoding="utf-8") as stream: + writer = csv.writer(stream) + writer.writerow(("gene_id", "dataset_type", "variant_id")) + writer.writerows(rows) + connection.execute( + """CREATE OR REPLACE TEMP TABLE __membership_batch AS + SELECT * FROM read_csv(?, header=true, columns={ + 'gene_id':'VARCHAR','dataset_type':'VARCHAR','variant_id':'VARCHAR' + })""", + [str(batch_csv)], + ) + connection.execute( + """INSERT OR IGNORE INTO compact_variant_membership + SELECT gene_id, dataset_type, variant_id FROM __membership_batch + WHERE variant_id IS NOT NULL AND trim(variant_id)<>''""" + ) + rows.clear() + + +def _read_membership_artifact( + task: tuple[str, str, str, str, int, int, str], +) -> dict[str, Any]: + """Read one immutable artifact without writing shared audit state.""" + + source, gene, dataset_type, logical, size_bytes, modified_at_ns, path_text = task + path = Path(path_text) + stat = path.stat() + if int(stat.st_size) != int(size_bytes) or int(stat.st_mtime_ns) != int(modified_at_ns): + raise RuntimeError(f"Registered source-summary artifact changed: {path}") + selected_rows = 0 + variant_ids: set[str] = set() + for entry in iter_json_array(path): + entry_sources = entry.get("sources") + if not isinstance(entry_sources, list): + entry_sources = [entry.get("source")] + normalized_sources = { + str(value).strip().casefold() for value in entry_sources if value + } + if str(source).casefold() not in normalized_sources: + continue + variant_id = str(entry.get("variant_id") or "").strip() + if not variant_id: + continue + selected_rows += 1 + variant_ids.add(variant_id) + return { + "source": str(source), + "gene": str(gene), + "dataset_type": str(dataset_type), + "logical": str(logical), + "size_bytes": int(size_bytes), + "modified_at_ns": int(modified_at_ns), + "selected_rows": selected_rows, + "membership_rows": [ + (str(gene), str(dataset_type), variant_id) + for variant_id in sorted(variant_ids) + ], + } + + +def materialize_compact_membership( + registry_connection: duckdb.DuckDBPyConnection, + *, + artifact_root: Path, + runtime_db_path: Path, + logger: logging.Logger, + progress_interval: int, + batch_rows: int = 100_000, + workers: int = 1, +) -> dict[str, int]: + """Build a resumable narrow index of unavailable-source variant membership.""" + + runtime_db_path.parent.mkdir(parents=True, exist_ok=True) + batch_csv = runtime_db_path.with_suffix(".membership-batch.csv") + runtime = duckdb.connect(str(runtime_db_path)) + runtime.execute( + """CREATE TABLE IF NOT EXISTS compact_variant_membership ( + gene_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + variant_id VARCHAR NOT NULL, + PRIMARY KEY (gene_id, dataset_type, variant_id) + ); + CREATE TABLE IF NOT EXISTS processed_source_summary_artifacts ( + logical_artifact VARCHAR PRIMARY KEY, + source VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + size_bytes UBIGINT NOT NULL, + modified_at_ns UBIGINT NOT NULL, + selected_source_rows UBIGINT NOT NULL, + processed_at TIMESTAMP NOT NULL + );""" + ) + processed = { + row[0] + for row in runtime.execute( + "SELECT logical_artifact FROM processed_source_summary_artifacts" + ).fetchall() + } + artifacts = registry_connection.execute( + """SELECT source, gene_id, dataset_type, logical_artifact, size_bytes, + modified_at_ns + FROM source_summary_artifacts + WHERE provider_detail_status='unavailable' + AND publication_mode='gene_scoped_on_demand' + ORDER BY dataset_type, gene_id, logical_artifact""" + ).fetchall() + started = time.perf_counter() + completed = len(processed) + selected_rows_total = int( + runtime.execute( + "SELECT coalesce(sum(selected_source_rows), 0) FROM processed_source_summary_artifacts" + ).fetchone()[0] + ) + tasks = [] + for source, gene, dataset_type, logical, size_bytes, modified_at_ns in artifacts: + if logical in processed: + continue + path = (artifact_root / str(logical)).resolve() + try: + path.relative_to(artifact_root.resolve()) + except ValueError as exc: + raise RuntimeError(f"Artifact escapes configured root: {logical}") from exc + tasks.append(( + str(source), str(gene), str(dataset_type), str(logical), + int(size_bytes), int(modified_at_ns), str(path), + )) + + def commit_result(result: dict[str, Any]) -> None: + nonlocal completed, selected_rows_total + runtime.execute("BEGIN TRANSACTION") + try: + rows = result["membership_rows"] + for offset in range(0, len(rows), max(1, batch_rows)): + batch = rows[offset : offset + max(1, batch_rows)] + _flush_membership_rows(runtime, batch, batch_csv) + runtime.execute( + """INSERT INTO processed_source_summary_artifacts VALUES + (?, ?, ?, ?, ?, ?, ?, current_timestamp)""", + [ + result["logical"], result["source"], result["gene"], + result["dataset_type"], result["size_bytes"], + result["modified_at_ns"], result["selected_rows"], + ], + ) + runtime.execute("COMMIT") + except Exception: + runtime.execute("ROLLBACK") + raise + finally: + batch_csv.unlink(missing_ok=True) + completed += 1 + selected_rows_total += int(result["selected_rows"]) + processed.add(str(result["logical"])) + if completed % max(1, progress_interval) == 0: + logger.info( + "runtime membership progress artifacts=%d/%d selected_rows=%d elapsed=%.1fs workers=%d", + completed, len(artifacts), selected_rows_total, + time.perf_counter() - started, max(1, workers), + ) + + worker_count = max(1, int(workers)) + if worker_count == 1: + for task in tasks: + commit_result(_read_membership_artifact(task)) + else: + executor = concurrent.futures.ProcessPoolExecutor(max_workers=worker_count) + task_iterator = iter(tasks) + futures: dict[concurrent.futures.Future, tuple] = {} + + def submit_next() -> bool: + try: + task = next(task_iterator) + except StopIteration: + return False + futures[executor.submit(_read_membership_artifact, task)] = task + return True + + try: + for _ in range(worker_count * 2): + if not submit_next(): + break + while futures: + done, _ = concurrent.futures.wait( + futures, return_when=concurrent.futures.FIRST_COMPLETED + ) + for future in done: + futures.pop(future) + commit_result(future.result()) + submit_next() + except BaseException: + for future in futures: + future.cancel() + executor.shutdown(wait=True, cancel_futures=True) + raise + else: + executor.shutdown(wait=True) + distinct_membership = int( + runtime.execute("SELECT count(*) FROM compact_variant_membership").fetchone()[0] + ) + runtime.close() + logger.info( + "runtime membership complete artifacts=%d distinct_gene_dataset_variants=%d selected_rows=%d elapsed=%.1fs", + completed, distinct_membership, selected_rows_total, + time.perf_counter() - started, + ) + return { + "artifacts": completed, + "distinct_gene_dataset_variants": distinct_membership, + "selected_source_rows": selected_rows_total, + } + + + +def materialize_indexed_compact_membership( + *, + serving_root: Path, + runtime_db_path: Path, + logger: logging.Logger, +) -> dict[str, int]: + """Build resumable membership from a completed gene-keyed serving index.""" + + serving_root = serving_root.resolve() + manifest_path = serving_root / "serving-manifest.json" + manifest_checksum_path = serving_root / "serving-manifest.json.sha256" + file_manifest_path = serving_root / "source-summary-index-files.jsonl" + manifest = json.loads(manifest_path.read_text()) + table = (manifest.get("tables") or {}).get( + "source_summary_associations_by_gene" + ) or {} + compatible_schemas = {"2.9.0-rc2", "2.10.0-rc2"} + if manifest.get("schema_version") not in compatible_schemas: + raise RuntimeError( + "Source-summary serving index is not a completed compatible release " + f"({', '.join(sorted(compatible_schemas))})." + ) + if table.get("contract") != "source_summary_association_index_v2": + raise RuntimeError("Unexpected source-summary serving-index contract.") + if int(table.get("files") or 0) != int(table.get("source_artifacts") or -1): + raise RuntimeError("Source-summary index file/artifact counts do not reconcile.") + expected_manifest_checksum = manifest_checksum_path.read_text().split()[0] + if sha256_file(manifest_path) != expected_manifest_checksum: + raise RuntimeError("Source-summary serving manifest checksum failed.") + if not file_manifest_path.exists(): + raise RuntimeError("Source-summary index file manifest is missing.") + + rollup_table = (manifest.get("tables") or {}).get( + "source_summary_association_base_by_gene" + ) or {} + use_rollup = ( + manifest.get("schema_version") == "2.10.0-rc2" + and rollup_table.get("contract") + == "source_summary_association_rollup_v2" + ) + if use_rollup and not ( + serving_root / "source-summary-rollup-files.jsonl" + ).is_file(): + raise RuntimeError("Source-summary rollup file manifest is missing.") + membership_table = rollup_table if use_rollup else table + logger.info( + "indexed membership source=%s", + "exact_default_rollup" if use_rollup else "full_source_summary_index", + ) + + runtime_db_path.parent.mkdir(parents=True, exist_ok=True) + runtime = duckdb.connect(str(runtime_db_path)) + runtime.execute("SET preserve_insertion_order=false") + runtime.execute("SET memory_limit='4GB'") + runtime.execute( + f"SET temp_directory='{runtime_db_path.parent.as_posix()}/duckdb-audit-tmp'" + ) + runtime.execute( + """CREATE TABLE IF NOT EXISTS compact_variant_membership ( + gene_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + variant_id VARCHAR NOT NULL, + PRIMARY KEY (gene_id, dataset_type, variant_id) + ); + CREATE TABLE IF NOT EXISTS processed_source_summary_artifacts ( + logical_artifact VARCHAR PRIMARY KEY, + source VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + size_bytes UBIGINT NOT NULL, + modified_at_ns UBIGINT NOT NULL, + selected_source_rows UBIGINT NOT NULL, + processed_at TIMESTAMP NOT NULL + ); + CREATE TABLE IF NOT EXISTS indexed_source_summary_buckets ( + gene_bucket VARCHAR PRIMARY KEY, + membership_rows_after UBIGINT NOT NULL, + processed_at TIMESTAMP NOT NULL + );""" + ) + runtime.execute( + """INSERT OR IGNORE INTO processed_source_summary_artifacts + SELECT logical_artifact, source, gene_id, dataset_type, + 0, 0, selected_source_rows, current_timestamp + FROM read_json_auto(?, format='newline_delimited')""", + [str(file_manifest_path)], + ) + processed_buckets = { + row[0] + for row in runtime.execute( + "SELECT gene_bucket FROM indexed_source_summary_buckets" + ).fetchall() + } + table_root = serving_root / str(membership_table["path"]) + started = time.perf_counter() + for bucket_dir in sorted(table_root.glob("gene_bucket=*")): + bucket = bucket_dir.name.split("=", 1)[1] + if bucket in processed_buckets: + continue + parquet_glob = str(bucket_dir / "gene_key=*" / "*.parquet") + runtime.execute("BEGIN TRANSACTION") + try: + runtime.execute( + """INSERT OR IGNORE INTO compact_variant_membership + SELECT gene_id, dataset_type, variant_id + FROM read_parquet(?, union_by_name=true) + WHERE variant_id IS NOT NULL AND trim(variant_id)<>''""", + [parquet_glob], + ) + membership_rows = int( + runtime.execute( + "SELECT count(*) FROM compact_variant_membership" + ).fetchone()[0] + ) + runtime.execute( + "INSERT INTO indexed_source_summary_buckets VALUES (?, ?, current_timestamp)", + [bucket, membership_rows], + ) + runtime.execute("COMMIT") + except Exception: + runtime.execute("ROLLBACK") + raise + processed_buckets.add(bucket) + logger.info( + "indexed membership progress buckets=%d/256 membership_rows=%d elapsed=%.1fs", + len(processed_buckets), membership_rows, time.perf_counter() - started, + ) + artifacts = int( + runtime.execute( + "SELECT count(*) FROM processed_source_summary_artifacts" + ).fetchone()[0] + ) + selected_rows = int( + runtime.execute( + "SELECT coalesce(sum(selected_source_rows),0) FROM processed_source_summary_artifacts" + ).fetchone()[0] + ) + distinct_membership = int( + runtime.execute("SELECT count(*) FROM compact_variant_membership").fetchone()[0] + ) + bucket_count = int( + runtime.execute("SELECT count(*) FROM indexed_source_summary_buckets").fetchone()[0] + ) + runtime.close() + if artifacts != int(table["source_artifacts"]): + raise RuntimeError("Indexed artifact count does not match the serving manifest.") + if selected_rows != int(table["selected_source_rows"]): + raise RuntimeError("Indexed selected-row count does not match the serving manifest.") + if bucket_count != 256: + raise RuntimeError(f"Expected 256 completed gene buckets, observed {bucket_count}.") + return { + "artifacts": artifacts, + "distinct_gene_dataset_variants": distinct_membership, + "selected_source_rows": selected_rows, + "indexed_gene_buckets": bucket_count, + } + +def v2_counts( + connection: duckdb.DuckDBPyConnection, + *, + runtime_db_path: Path | None = None, + v2_db_path: Path | None = None, +) -> tuple[ + dict[tuple[str, str, str, str], int], + dict[tuple[str, str], dict[str, int]], +]: + query_connection = connection + prefix = "" + membership_cte = """ + SELECT DISTINCT gene_id, dataset_type, variant_id + FROM variant_phenotype_summaries + """ + runtime = None + if runtime_db_path is not None: + if v2_db_path is None: + raise ValueError("v2_db_path is required with runtime_db_path") + runtime = duckdb.connect(str(runtime_db_path)) + escaped_v2_path = str(v2_db_path).replace("'", "''") + runtime.execute(f"ATTACH '{escaped_v2_path}' AS v2 (READ_ONLY)") + query_connection = runtime + prefix = "v2." + membership_cte = """ + SELECT gene_id, dataset_type, variant_id + FROM v2.variant_phenotype_summaries + UNION + SELECT gene_id, dataset_type, variant_id + FROM compact_variant_membership + """ + category: dict[tuple[str, str, str, str], int] = {} + queries = { + "vc": f""" + WITH membership AS ({membership_cte}) + SELECT m.gene_id, m.dataset_type, + coalesce(v.variation_type, 'Unresolved'), + count(DISTINCT m.variant_id) + FROM membership m + LEFT JOIN {prefix}variants v USING (variant_id) + GROUP BY 1, 2, 3 + """, + "msc": f""" + WITH membership AS ({membership_cte}) + SELECT m.gene_id, m.dataset_type, a.consequence, + count(DISTINCT m.variant_id) + FROM membership m + JOIN {prefix}consequence_annotations a + ON a.variant_id=m.variant_id AND a.gene_id=m.gene_id + GROUP BY 1, 2, 3 + """, + "cs": f""" + WITH membership AS ({membership_cte}) + SELECT m.gene_id, m.dataset_type, a.clinical_significance, + count(DISTINCT m.variant_id) + FROM membership m + JOIN {prefix}clinical_assertions a USING (variant_id) + GROUP BY 1, 2, 3 + """, + } + try: + for dimension, query in queries.items(): + for gene, dataset_type, label, count in query_connection.execute(query).fetchall(): + category[(gene, dataset_type, dimension, normalized_label(label))] = int(count) + for dimension, annotation_table, join_clause in ( + ("msc", "consequence_annotations", "a.variant_id=m.variant_id AND a.gene_id=m.gene_id"), + ("cs", "clinical_assertions", "a.variant_id=m.variant_id"), + ): + for gene, dataset_type, count in query_connection.execute( + f""" + WITH membership AS ({membership_cte}) + SELECT m.gene_id, m.dataset_type, count(DISTINCT m.variant_id) + FROM membership m + WHERE NOT EXISTS ( + SELECT 1 FROM {prefix}{annotation_table} a WHERE {join_clause} + ) + GROUP BY 1, 2 + """ + ).fetchall(): + category[(gene, dataset_type, dimension, "unannotated")] = int(count) + + provider_totals = { + (row[0], row[1]): { + "source_observation_count": int(row[2] or 0), + "association_record_count": int(row[3] or 0), + "provider_association_record_count": int(row[4] or 0), + } + for row in query_connection.execute( + f""" + SELECT gene_id, dataset_type, sum(source_observation_count), + sum(association_record_count), + sum(provider_association_record_count) + FROM {prefix}variant_phenotype_summaries + GROUP BY 1, 2 + """ + ).fetchall() + } + compact_totals = {} + if runtime is not None: + compact_totals = { + (row[0], row[1]): { + "retained_source_summary_count": int(row[2] or 0), + "retained_source_summary_artifacts": int(row[3] or 0), + } + for row in query_connection.execute( + """SELECT gene_id, dataset_type, sum(selected_source_rows), count(*) + FROM processed_source_summary_artifacts GROUP BY 1, 2""" + ).fetchall() + } + totals = {} + for gene, dataset_type, count in query_connection.execute( + f"""WITH membership AS ({membership_cte}) + SELECT gene_id, dataset_type, count(DISTINCT variant_id) + FROM membership GROUP BY 1, 2""" + ).fetchall(): + key = (gene, dataset_type) + totals[key] = { + "v2_distinct_variants": int(count), + **provider_totals.get(key, { + "source_observation_count": 0, + "association_record_count": 0, + "provider_association_record_count": 0, + }), + **compact_totals.get(key, { + "retained_source_summary_count": 0, + "retained_source_summary_artifacts": 0, + }), + } + return category, totals + finally: + if runtime is not None: + runtime.close() + +def v1_counts( + root: Path, + *, + logger: logging.Logger, + progress_interval: int, + max_files: int | None, +) -> tuple[ + dict[tuple[str, str, str, str], int], + dict[tuple[str, str], int], +]: + category: dict[tuple[str, str, str, str], int] = {} + totals: dict[tuple[str, str], int] = {} + files_seen = 0 + started = time.perf_counter() + for dataset_type in ("CVD", "TRAIT"): + directory = root / dataset_type + paths = sorted(directory.glob("*.json*")) + if max_files is not None: + paths = paths[:max_files] + for path in paths: + files_seen += 1 + opener = gzip.open if path.name.endswith(".gz") else open + with opener(path, "rt", encoding="utf-8") as stream: + payload = json.load(stream) + gene = path.name.removesuffix(".gz").removesuffix(".json") + data = payload.get("data") if isinstance(payload, dict) else {} + data = data if isinstance(data, dict) else {} + for dimension in DIMENSIONS: + values = data.get(dimension) if isinstance(data.get(dimension), dict) else {} + for label, count in values.items(): + key = (gene, dataset_type, dimension, normalized_label(label)) + category[key] = category.get(key, 0) + int(count or 0) + totals[(gene, dataset_type)] = sum( + int(value or 0) + for value in (data.get("vc") or {}).values() + ) + if files_seen % max(1, progress_interval) == 0: + logger.info( + "v1 count progress files=%d current=%s elapsed=%.1fs", + files_seen, + path, + time.perf_counter() - started, + ) + logger.info("v1 count complete files=%d elapsed=%.1fs", files_seen, time.perf_counter() - started) + return category, totals + + +def compact_variant_index_entries( + root: Path, + gene: str, + dataset_type: str, + phenotype_slug: str, + variant_id: str, +) -> list[dict[str, Any]]: + directory = root / dataset_type + candidates = (directory / f"{gene}.json.gz", directory / f"{gene}.json") + path = next((candidate for candidate in candidates if candidate.exists()), None) + if path is None: + return [] + opener = gzip.open if path.name.endswith(".gz") else open + with opener(path, "rt", encoding="utf-8") as stream: + payload = json.load(stream) + matches = [] + for row in payload if isinstance(payload, list) else []: + row_path = row.get("phenotype_path") or row.get("disease") or row.get("trait") or [] + row_slug = str(row.get("phenotype") or (row_path[-1] if row_path else "")).strip().lower() + if str(row.get("variant_id") or "") == variant_id and row_slug == phenotype_slug: + matches.append(row) + return matches + + +def example_audit( + connection: duckdb.DuckDBPyConnection, + variant_index_root: Path, +) -> list[dict[str, Any]]: + results = [] + for gene, dataset_type, phenotype_slug, variant_id in EXAMPLES: + summary = connection.execute( + """ + SELECT * FROM variant_phenotype_summaries + WHERE gene_id=? AND dataset_type=? AND phenotype_slug=? AND variant_id=? + """, + [gene, dataset_type, phenotype_slug, variant_id], + ).fetchone() + if summary is None: + results.append( + { + "gene": gene, + "dataset_type": dataset_type, + "phenotype_slug": phenotype_slug, + "variant_id": variant_id, + "status": "missing", + } + ) + continue + columns = [ + row[0] + for row in connection.execute( + "DESCRIBE variant_phenotype_summaries" + ).fetchall() + ] + value = dict(zip(columns, summary)) + summary_id = value["variant_phenotype_summary_id"] + source_rows = connection.execute( + """ + SELECT ar.source, ar.record_kind, count(*), sum(ar.provider_record_count), + min(ar.reported_p_value) + FROM summary_association_records l + JOIN association_records ar USING (association_record_id) + WHERE l.variant_phenotype_summary_id=? + GROUP BY 1, 2 ORDER BY 1, 2 + """, + [summary_id], + ).fetchall() + provider_ids = connection.execute( + """ + SELECT p.provider_record_id + FROM summary_association_records s + JOIN association_record_provider_records l USING (association_record_id) + JOIN provider_records p USING (provider_record_id) + WHERE s.variant_phenotype_summary_id=? + ORDER BY p.provider_record_id + """, + [summary_id], + ).fetchall() + compact_entries = compact_variant_index_entries( + variant_index_root, gene, dataset_type, phenotype_slug, variant_id + ) + results.append( + { + "gene": gene, + "dataset_type": dataset_type, + "phenotype_slug": phenotype_slug, + "variant_id": variant_id, + "status": "verified", + "summary_id": summary_id, + "source_observation_count": int(value["source_observation_count"]), + "association_record_count": int(value["association_record_count"]), + "provider_record_count": len(provider_ids), + "minimum_reported_p_value": value["minimum_reported_p_value"], + "minimum_p_supplier": value[ + "minimum_reported_p_value_association_record_id" + ], + "consequence_annotation_count": int( + value["consequence_annotation_count"] + ), + "clinical_assertion_count": int(value["clinical_assertion_count"]), + "source_breakdown": [ + { + "source": row[0], + "record_kind": row[1], + "association_records": int(row[2]), + "provider_rows": int(row[3]), + "minimum_reported_p_value": row[4], + } + for row in source_rows + ], + "provider_record_ids": [row[0] for row in provider_ids], + "published_v1_compact_entries": compact_entries, + "published_v1_representative_sources": sorted( + {str(row.get("source") or "") for row in compact_entries if row.get("source")} + ), + "published_v1_all_sources": sorted( + { + str(source) + for row in compact_entries + for source in (row.get("sources") or [row.get("source")]) + if source + } + ), + "source_priority_interpretation": ( + "The v1 compact entry exposes one representative source/annotation. " + "The v2 result preserves every recoverable provider row and versioned " + "annotation; minimum_reported_p_value is aggregate-only and records its supplier." + ), + } + ) + return results + + +SUPPLEMENTAL_QA_QUERIES = { + "association_gene_links_resolve": """ + SELECT count(*) FROM association_record_genes g + LEFT JOIN association_records a USING (association_record_id) + WHERE a.association_record_id IS NULL + """, + "consequence_link_associations_resolve": """ + SELECT count(*) FROM consequence_annotation_provider_records l + LEFT JOIN association_records a USING (association_record_id) + WHERE a.association_record_id IS NULL + """, + "clinical_link_associations_resolve": """ + SELECT count(*) FROM clinical_assertion_provider_records l + LEFT JOIN association_records a USING (association_record_id) + WHERE a.association_record_id IS NULL + """, + "population_link_associations_resolve": """ + SELECT count(*) FROM population_observation_provider_records l + LEFT JOIN association_records a USING (association_record_id) + WHERE a.association_record_id IS NULL + """, + "summary_consequence_links_resolve": """ + SELECT count(*) FROM summary_consequence_annotations l + LEFT JOIN variant_phenotype_summaries s USING (variant_phenotype_summary_id) + LEFT JOIN consequence_annotations a USING (consequence_annotation_id) + WHERE s.variant_phenotype_summary_id IS NULL OR a.consequence_annotation_id IS NULL + """, + "summary_clinical_links_resolve": """ + SELECT count(*) FROM summary_clinical_assertions l + LEFT JOIN variant_phenotype_summaries s USING (variant_phenotype_summary_id) + LEFT JOIN clinical_assertions a USING (clinical_assertion_id) + WHERE s.variant_phenotype_summary_id IS NULL OR a.clinical_assertion_id IS NULL + """, + "summary_variants_resolve": """ + SELECT count(*) FROM variant_phenotype_summaries s + LEFT JOIN variants v USING (variant_id) + WHERE v.variant_id IS NULL + """, + "no_compact_rows_promoted_to_associations": """ + SELECT count(*) FROM association_records + WHERE record_kind='source_summary' OR provider_detail_status='unavailable' + """, + "no_compact_rows_promoted_to_provider_records": """ + SELECT count(*) FROM provider_records + WHERE source='million_veteran_program' + """, +} + + +def supplemental_qa(connection: duckdb.DuckDBPyConnection) -> list[dict[str, Any]]: + rows = [] + for name, query in SUPPLEMENTAL_QA_QUERIES.items(): + violations = int(connection.execute(query).fetchone()[0]) + rows.append({ + "check_name": name, + "status": "passed" if violations == 0 else "failed", + "violations": violations, + "expected": 0, + }) + return rows + + +def main() -> int: + args = parse_args() + args.output_dir.mkdir(parents=True, exist_ok=True) + checkpoint_path = args.checkpoint_path or args.output_dir / "audit-checkpoint.json" + log_path = args.log_path or args.output_dir / "audit.log" + v2_intermediate = args.output_dir / "intermediate-v2-counts.json" + v1_intermediate = args.output_dir / "intermediate-v1-counts.json" + runtime_db_path = args.output_dir / "runtime-compact-membership.duckdb" + artifact_root = (args.artifact_root or args.v1_variant_index_root.parent).resolve() + if args.reset: + for path in (checkpoint_path, v2_intermediate, v1_intermediate, runtime_db_path): + path.unlink(missing_ok=True) + logger = configure_logging(log_path) + configuration = configuration_hash(args) + checkpoint = load_checkpoint(checkpoint_path, configuration) + logger.info("audit start v2=%s v1=%s output=%s", args.v2_db, args.v1_overall_root, args.output_dir) + + connection = duckdb.connect(str(args.v2_db), read_only=True) + if checkpoint.get("phase") in {"v2_counts_complete", "v1_counts_complete", "complete"} and v2_intermediate.exists(): + logger.info("resume: loading v2 count intermediate %s", v2_intermediate) + v2_category, v2_totals = decode_count_bundle(json.loads(v2_intermediate.read_text()), rich_totals=True) + else: + started = time.perf_counter() + logger.info("v2 runtime membership phase start artifact_root=%s", artifact_root) + if args.source_summary_serving_root: + logger.info( + "reusing completed source-summary serving index root=%s", + args.source_summary_serving_root, + ) + membership_stats = materialize_indexed_compact_membership( + serving_root=args.source_summary_serving_root, + runtime_db_path=runtime_db_path, + logger=logger, + ) + else: + membership_stats = materialize_compact_membership( + connection, + artifact_root=artifact_root, + runtime_db_path=runtime_db_path, + logger=logger, + progress_interval=args.progress_interval, + workers=args.workers, + ) + logger.info("v2 count phase start") + v2_category, v2_totals = v2_counts( + connection, + runtime_db_path=runtime_db_path, + v2_db_path=args.v2_db, + ) + v2_intermediate.write_text(json.dumps(encode_count_bundle(v2_category, v2_totals), sort_keys=True) + "\n") + save_checkpoint(checkpoint_path, configuration, "v2_counts_complete", { + "category_rows": len(v2_category), "gene_dataset_rows": len(v2_totals), + "runtime_membership": membership_stats, + "elapsed_seconds": round(time.perf_counter() - started, 3), + }) + checkpoint = load_checkpoint(checkpoint_path, configuration) + logger.info("v2 count phase complete categories=%d elapsed=%.1fs", len(v2_category), time.perf_counter() - started) + examples = example_audit(connection, args.v1_variant_index_root) + qa = [ + dict(zip(("check_name", "status", "observed", "expected"), row)) + for row in connection.execute( + "SELECT check_name, status, observed_value, expected_value " + "FROM qa_results ORDER BY check_name" + ).fetchall() + ] + supplemental_qa_rows: list[dict[str, Any]] = [] + supplemental_qa_execution = { + "status": "not_reexecuted", + "reason": ( + "The release audit reuses the persisted full-build QA results. " + "A second set of referential joins over the normalized link tables was " + "intentionally not rerun because it duplicates build-time checks and " + "requires an unbounded rescan of up to one billion link rows." + ), + "persisted_build_qa_checks": len(qa), + } + metadata = connection.execute("SELECT * FROM build_metadata").fetchone() + metadata_columns = [ + row[0] for row in connection.execute("DESCRIBE build_metadata").fetchall() + ] + connection.close() + + if checkpoint.get("phase") in {"v1_counts_complete", "complete"} and v1_intermediate.exists(): + logger.info("resume: loading v1 count intermediate %s", v1_intermediate) + v1_category, v1_totals = decode_count_bundle(json.loads(v1_intermediate.read_text()), rich_totals=False) + else: + v1_category, v1_totals = v1_counts( + args.v1_overall_root, + logger=logger, + progress_interval=args.progress_interval, + max_files=args.max_v1_files, + ) + v1_intermediate.write_text(json.dumps(encode_count_bundle(v1_category, v1_totals), sort_keys=True) + "\n") + save_checkpoint(checkpoint_path, configuration, "v1_counts_complete", { + "category_rows": len(v1_category), "gene_dataset_rows": len(v1_totals), + }) + comparison = [] + exact = defaultdict(int) + changed = defaultdict(int) + aggregate = defaultdict(lambda: {"v1": 0, "v2": 0}) + for key in sorted(set(v1_category) | set(v2_category)): + gene, dataset_type, dimension, label = key + v1 = v1_category.get(key, 0) + v2 = v2_category.get(key, 0) + matches = v1 == v2 + exact[dimension] += int(matches) + changed[dimension] += int(not matches) + aggregate[dimension]["v1"] += v1 + aggregate[dimension]["v2"] += v2 + comparison.append( + { + "gene": gene, + "dataset_type": dataset_type, + "dimension": dimension, + "normalized_label": label, + "v1_count": v1, + "v2_count": v2, + "delta": v2 - v1, + "status": "unchanged" if matches else "changed_expected_review", + "difference_contract": "" if matches else REASONS[dimension], + "difference_reason": category_difference_reason( + gene=gene, + dataset_type=dataset_type, + dimension=dimension, + v1_count=v1, + v2_count=v2, + v2_totals=v2_totals, + ), + } + ) + + gene_rows = [] + for key in sorted(set(v1_totals) | set(v2_totals)): + details = v2_totals.get(key, {}) + v1 = v1_totals.get(key, 0) + v2 = details.get("v2_distinct_variants", 0) + gene_rows.append( + { + "gene": key[0], + "dataset_type": key[1], + "v1_variation_axis_total": v1, + "v2_distinct_variants": v2, + "delta": v2 - v1, + **{k: value for k, value in details.items() if k != "v2_distinct_variants"}, + } + ) + + write_csv(args.output_dir / "category-count-comparison.csv", comparison) + write_csv(args.output_dir / "gene-count-comparison.csv", gene_rows) + report = { + "generated_at": datetime.now(timezone.utc).isoformat(), + "v2_database": str(args.v2_db.resolve()), + "v2_database_sha256": args.v2_sha256 or sha256_file(args.v2_db), + "v1_overall_root": str(args.v1_overall_root.resolve()), + "v1_variant_index_root": str(args.v1_variant_index_root.resolve()), + "artifact_root": str(artifact_root), + "source_summary_serving_root": ( + str(args.source_summary_serving_root.resolve()) + if args.source_summary_serving_root else None + ), + "runtime_compact_membership": str(runtime_db_path.resolve()), + "build_metadata": dict(zip(metadata_columns, metadata)), + "category_rows": len(comparison), + "gene_dataset_rows": len(gene_rows), + "exact_category_rows": dict(exact), + "changed_category_rows": dict(changed), + "aggregate_counts": dict(aggregate), + "five_example_validation": examples, + "qa_results": qa, + "supplemental_qa_results": supplemental_qa_rows, + "supplemental_qa_execution": supplemental_qa_execution, + "difference_contracts": REASONS, + } + report_path = args.output_dir / "release-audit.json" + report_path.write_text(json.dumps(report, indent=2, sort_keys=True, default=str) + "\n") + checksum_lines = [] + for path in sorted(args.output_dir.iterdir()): + if path.is_file() and path.name != "checksums.sha256": + checksum_lines.append(f"{sha256_file(path)} {path.name}") + (args.output_dir / "checksums.sha256").write_text( + "\n".join(checksum_lines) + "\n" + ) + save_checkpoint(checkpoint_path, configuration, "complete", { + "report": str(report_path), + "category_rows": len(comparison), + "gene_dataset_rows": len(gene_rows), + }) + logger.info("audit complete report=%s categories=%d genes=%d", report_path, len(comparison), len(gene_rows)) + print(json.dumps({ + "output_dir": str(args.output_dir), + "category_rows": len(comparison), + "gene_dataset_rows": len(gene_rows), + "five_examples_verified": sum(row["status"] == "verified" for row in examples), + "qa_passed": ( + all(row["status"] == "passed" for row in qa) + ), + }, indent=2)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2.py new file mode 100644 index 0000000..0442826 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python3 +"""Build the non-composite association evidence v2 sidecar. + +The builder is resumable at source-file and derived-phase boundaries. It writes +only to the explicitly supplied sidecar database, checkpoint, manifest, and log +paths. MVP compact evidence is represented as first-class source-summary associations +without fabricating provider or study records. +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import logging +import shlex +import subprocess +import sys +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +REPO_ROOT = Path(__file__).resolve().parents[3] +sys.path.insert(0, str(REPO_ROOT / "src")) + +from datahub.association_evidence_v2.builder import EvidenceV2Builder # noqa: E402 + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--cvd-root", required=True, type=Path) + parser.add_argument("--trait-root", required=True, type=Path) + parser.add_argument("--variant-index-root", required=True, type=Path) + parser.add_argument("--output-db", required=True, type=Path) + parser.add_argument("--checkpoint-path", required=True, type=Path) + parser.add_argument("--manifest-path", required=True, type=Path) + parser.add_argument( + "--phenotype-tree", + type=Path, + default=REPO_ROOT / "config" / "phenotype_tree.json", + ) + parser.add_argument("--release-id", required=True) + parser.add_argument("--threads", type=int, default=4) + parser.add_argument("--memory-limit", default="12GB") + parser.add_argument("--temp-directory", type=Path) + parser.add_argument("--progress-interval", type=int, default=1) + parser.add_argument( + "--provider-chunk-rows", + type=int, + default=250_000, + help=( + "Provider rows per bounded transaction; 0 uses direct inserts up to 2 GB " + "and 1,000,000-row one-pass chunks for larger sources." + ), + ) + parser.add_argument( + "--smoke-files", + type=int, + help="Process at most this many source files per dataset type.", + ) + parser.add_argument( + "--smoke-rows-per-file", + type=int, + help="Process at most this many rows from each selected source file.", + ) + parser.add_argument("--no-resume", action="store_true") + parser.add_argument("--log-path", type=Path) + parser.add_argument( + "--checksum-inputs", + action="store_true", + help="SHA-256 source files after the build. This performs a second full read.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Validate paths and print the planned configuration without writing.", + ) + return parser.parse_args() + + +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def write_sha256_sidecar(path: Path, digest: str) -> Path: + """Write the checksum contract consumed by the serving publisher.""" + checksum_path = path.with_suffix(path.suffix + ".sha256") + checksum_path.write_text(f"{digest} {path.name}\n") + return checksum_path + + +def git_commit() -> str | None: + result = subprocess.run( + ["git", "-C", str(REPO_ROOT), "rev-parse", "HEAD"], + capture_output=True, + text=True, + check=False, + ) + return result.stdout.strip() or None + + +def git_dirty() -> bool: + result = subprocess.run( + ["git", "-C", str(REPO_ROOT), "status", "--porcelain"], + capture_output=True, + text=True, + check=False, + ) + return bool(result.stdout.strip()) + + +def transformation_files() -> list[Path]: + package_root = REPO_ROOT / "src" / "datahub" / "association_evidence_v2" + return [ + Path(__file__).resolve(), + *sorted(package_root.glob("*.py")), + ] + + +def source_files(root: Path, limit: int | None) -> list[Path]: + files = sorted(path for path in root.iterdir() if path.is_file()) + return files[:limit] if limit is not None else files + + +def configure_logging(path: Path | None) -> None: + handlers: list[logging.Handler] = [logging.StreamHandler()] + if path: + path.parent.mkdir(parents=True, exist_ok=True) + handlers.append(logging.FileHandler(path)) + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + handlers=handlers, + ) + + +def main() -> int: + args = parse_args() + configure_logging(args.log_path) + logger = logging.getLogger("datahub.association_evidence_v2.cli") + for path in (args.cvd_root, args.trait_root, args.variant_index_root, args.phenotype_tree): + if not path.exists(): + raise FileNotFoundError(path) + command = " ".join(shlex.quote(item) for item in sys.argv) + datahub_commit = git_commit() + plan = { + "release_id": args.release_id, + "cvd_root": str(args.cvd_root.resolve()), + "trait_root": str(args.trait_root.resolve()), + "variant_index_root": str(args.variant_index_root.resolve()), + "output_db": str(args.output_db.resolve()), + "checkpoint_path": str(args.checkpoint_path.resolve()), + "manifest_path": str(args.manifest_path.resolve()), + "phenotype_tree": str(args.phenotype_tree.resolve()), + "source_files": { + "CVD": len(source_files(args.cvd_root, args.smoke_files)), + "TRAIT": len(source_files(args.trait_root, args.smoke_files)), + }, + "smoke_files": args.smoke_files, + "smoke_rows_per_file": args.smoke_rows_per_file, + "threads": args.threads, + "memory_limit": args.memory_limit, + "provider_chunk_rows": args.provider_chunk_rows, + "temp_directory": str(args.temp_directory.resolve()) if args.temp_directory else None, + "command": command, + "datahub_git_commit": datahub_commit, + "datahub_git_dirty": git_dirty(), + } + logger.info("Evidence v2 build plan: %s", json.dumps(plan, sort_keys=True)) + if args.dry_run: + print(json.dumps({"mode": "dry_run", **plan}, indent=2, sort_keys=True)) + return 0 + + builder = EvidenceV2Builder( + cvd_root=args.cvd_root, + trait_root=args.trait_root, + variant_index_root=args.variant_index_root, + output_db=args.output_db, + checkpoint_path=args.checkpoint_path, + phenotype_tree=args.phenotype_tree, + release_id=args.release_id, + command=command, + builder_git_commit=datahub_commit, + threads=args.threads, + memory_limit=args.memory_limit, + temp_directory=args.temp_directory, + smoke_files=args.smoke_files, + smoke_rows_per_file=args.smoke_rows_per_file, + provider_chunk_rows=args.provider_chunk_rows, + progress_interval=args.progress_interval, + ) + summary = builder.run(resume=not args.no_resume) + inputs: list[dict[str, Any]] = [] + for dataset_type, root in (("CVD", args.cvd_root), ("TRAIT", args.trait_root)): + for path in source_files(root, args.smoke_files): + item = { + "dataset_type": dataset_type, + "path": str(path.resolve()), + "size_bytes": path.stat().st_size, + "sha256": None, + "checksum_status": "not_requested", + } + if args.checksum_inputs: + item["sha256"] = sha256_file(path) + item["checksum_status"] = "verified" + inputs.append(item) + variant_index_files = sorted( + path + for path in args.variant_index_root.rglob("*") + if path.is_file() and (path.name.endswith(".json") or path.name.endswith(".json.gz")) + ) + transformation_inventory = [ + { + "path": str(path.relative_to(REPO_ROOT)), + "size_bytes": path.stat().st_size, + "sha256": sha256_file(path), + } + for path in transformation_files() + ] + output_sha256 = sha256_file(args.output_db) + output_checksum_path = write_sha256_sidecar(args.output_db, output_sha256) + manifest = { + **summary, + "generated_at": datetime.now(timezone.utc).isoformat(), + "datahub_git_commit": datahub_commit, + "datahub_git_dirty": git_dirty(), + "builder_command": command, + "transformation_files": transformation_inventory, + "phenotype_tree": { + "path": str(args.phenotype_tree.resolve()), + "size_bytes": args.phenotype_tree.stat().st_size, + "sha256": sha256_file(args.phenotype_tree), + }, + "variant_index_input": { + "path": str(args.variant_index_root.resolve()), + "file_count": len(variant_index_files), + "size_bytes": sum(path.stat().st_size for path in variant_index_files), + "checksum_status": "not_requested_to_avoid_second_full_read", + }, + "inputs": inputs, + "output": { + "path": str(args.output_db.resolve()), + "size_bytes": args.output_db.stat().st_size, + "sha256": output_sha256, + "checksum_path": str(output_checksum_path.resolve()), + "checksum_file_sha256": sha256_file(output_checksum_path), + }, + "checkpoint_sha256": sha256_file(args.checkpoint_path), + } + args.manifest_path.parent.mkdir(parents=True, exist_ok=True) + args.manifest_path.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n") + print(json.dumps(manifest, indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_serving.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_serving.py new file mode 100644 index 0000000..865d3d9 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_serving.py @@ -0,0 +1,69 @@ +#!/usr/bin/env python3 +"""Publish bounded serving partitions for an association evidence v2 sidecar.""" + +from __future__ import annotations + +import argparse +import json +import logging +import sys +from pathlib import Path + +ROOT = Path(__file__).resolve().parents[3] +SRC = ROOT / "src" +if str(SRC) not in sys.path: + sys.path.insert(0, str(SRC)) + +from datahub.association_evidence_v2.serving import AssociationEvidenceV2ServingBuilder + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--source-db", type=Path, required=True) + parser.add_argument("--intermediate-root", type=Path, required=True) + parser.add_argument("--output-root", type=Path, required=True) + parser.add_argument("--memory-limit", default="8GB") + parser.add_argument("--threads", type=int, default=1) + parser.add_argument( + "--variant-bucket-characters", + type=int, + choices=(2, 3, 4), + default=3, + help="Hexadecimal SHA-256 prefix width for variant-scoped runtime partitions.", + ) + parser.add_argument("--progress-interval", type=int, default=60) + parser.add_argument( + "--coarse-serving-root", + type=Path, + help=( + "Existing verified lower-width serving package used as the bounded " + "input when --variant-bucket-characters is greater than 2." + ), + ) + parser.add_argument("--verbose", action="store_true") + args = parser.parse_args() + args.output_root.mkdir(parents=True, exist_ok=True) + logging.basicConfig( + level=logging.DEBUG if args.verbose else logging.INFO, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + handlers=[ + logging.StreamHandler(), + logging.FileHandler(args.output_root / "serving-build.log", mode="a"), + ], + ) + builder = AssociationEvidenceV2ServingBuilder( + source_db=args.source_db, + intermediate_root=args.intermediate_root, + output_root=args.output_root, + memory_limit=args.memory_limit, + threads=args.threads, + variant_bucket_characters=args.variant_bucket_characters, + progress_interval=args.progress_interval, + coarse_serving_root=args.coarse_serving_root, + ) + print(json.dumps(builder.build(), indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_index.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_index.py new file mode 100644 index 0000000..838e821 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_index.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 +"""Build the resumable first-class MVP source-summary association index.""" + +from __future__ import annotations + +import argparse +import json +import logging +from pathlib import Path + +from datahub.association_evidence_v2.source_summary_index import SourceSummaryIndexBuilder + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--source-db", type=Path, required=True) + parser.add_argument("--association-artifact-root", type=Path, required=True) + parser.add_argument("--serving-root", type=Path, required=True) + parser.add_argument("--workers", type=int, default=3) + parser.add_argument("--progress-interval", type=int, default=1000) + parser.add_argument("--verbose", action="store_true") + args = parser.parse_args() + logging.basicConfig( + level=logging.DEBUG if args.verbose else logging.INFO, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + ) + result = SourceSummaryIndexBuilder( + source_db=args.source_db, + artifact_root=args.association_artifact_root, + serving_root=args.serving_root, + workers=args.workers, + progress_interval=args.progress_interval, + ).build() + print(json.dumps(result, indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_rollup.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_rollup.py new file mode 100755 index 0000000..7883933 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_source_summary_rollup.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python3 +"""Build default-search rollups from a completed evidence-v2 source-summary index.""" + +from __future__ import annotations + +import argparse +import json +import logging +from pathlib import Path + +from datahub.association_evidence_v2.source_summary_rollup import SourceSummaryRollupBuilder + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--serving-root", required=True, type=Path) + parser.add_argument("--threads", type=int, default=2) + parser.add_argument("--memory-limit", default="4GB") + parser.add_argument("--progress-interval", type=int, default=1) + parser.add_argument("--verbose", action="store_true") + args = parser.parse_args() + logging.basicConfig( + level=logging.INFO if args.verbose else logging.WARNING, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + ) + result = SourceSummaryRollupBuilder( + serving_root=args.serving_root, + threads=args.threads, + memory_limit=args.memory_limit, + progress_interval=args.progress_interval, + logger=logging.getLogger("datahub.association_evidence_v2.source_summary_rollup"), + ).build() + print(json.dumps(result, indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup.py new file mode 100755 index 0000000..8b50621 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup.py @@ -0,0 +1,239 @@ +#!/usr/bin/env python3 +"""Build a bounded canonical-variant lookup for evidence-model-v2 serving.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import logging +import shutil +import threading +import time +from datetime import datetime, timezone +from pathlib import Path + + +LOGGER = logging.getLogger("datahub.association_evidence_v2.variant_lookup") + + +def utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def write_json(path: Path, payload: dict) -> None: + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text(json.dumps(payload, indent=2, sort_keys=True) + "\n") + temporary.replace(path) + + +def monitor_progress(root: Path, interval: int, stop: threading.Event) -> None: + started = time.perf_counter() + while not stop.wait(interval): + files = 0 + size = 0 + if root.exists(): + for path in root.rglob("*.parquet"): + try: + stat = path.stat() + except FileNotFoundError: + continue + files += 1 + size += stat.st_size + LOGGER.info( + "variant lookup progress elapsed_seconds=%.1f parquet_files=%d bytes=%d", + time.perf_counter() - started, + files, + size, + ) + + +def build(args: argparse.Namespace) -> dict: + import duckdb + + source_db = args.source_db.resolve() + output_root = args.output_root.resolve() + checkpoint = output_root / "variant-lookup-checkpoint.json" + manifest_path = output_root / "serving-manifest.json" + checksum_path = output_root / "checksums.sha256" + table_root = output_root / "tables" / "variants" + incomplete_root = output_root / ".incomplete-variants" + + if not source_db.is_file(): + raise FileNotFoundError(f"Source sidecar not found: {source_db}") + if args.reset and output_root.exists(): + shutil.rmtree(output_root) + output_root.mkdir(parents=True, exist_ok=True) + + configuration = { + "source_db": str(source_db), + "source_db_sha256": args.source_sha256, + "bucket_characters": args.bucket_characters, + "smoke_limit": args.smoke_limit, + } + if checkpoint.is_file(): + state = json.loads(checkpoint.read_text()) + if state.get("configuration") != configuration: + raise RuntimeError("Existing checkpoint configuration does not match this run") + if state.get("complete") and table_root.is_dir() and manifest_path.is_file(): + LOGGER.info("variant lookup already complete output_root=%s", output_root) + return json.loads(manifest_path.read_text()) + else: + state = { + "configuration": configuration, + "complete": False, + "started_at": utc_now(), + } + write_json(checkpoint, state) + + if incomplete_root.exists(): + shutil.rmtree(incomplete_root) + if table_root.exists(): + raise RuntimeError(f"Uncheckpointed final table already exists: {table_root}") + incomplete_table = incomplete_root / "tables" / "variants" + incomplete_table.parent.mkdir(parents=True, exist_ok=True) + + connection = duckdb.connect( + str(source_db), + read_only=True, + config={"threads": str(args.threads), "memory_limit": args.memory_limit}, + ) + connection.execute("SET preserve_insertion_order=false") + available_rows = int( + connection.execute( + "SELECT estimated_size FROM duckdb_tables() WHERE table_name='variants'" + ).fetchone()[0] + ) + limit_sql = f" LIMIT {int(args.smoke_limit)}" if args.smoke_limit else "" + LOGGER.info( + "variant lookup build start source=%s output=%s available_rows=%d smoke_limit=%s", + source_db, + output_root, + available_rows, + args.smoke_limit, + ) + stop = threading.Event() + monitor = threading.Thread( + target=monitor_progress, + args=(incomplete_table, args.progress_interval, stop), + daemon=True, + ) + monitor.start() + started = time.perf_counter() + try: + copied = connection.execute( + f""" +COPY ( + SELECT *, substr(sha256(variant_id), 1, {args.bucket_characters}) AS variant_bucket + FROM variants + {limit_sql} +) TO '{str(incomplete_table).replace("'", "''")}' +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ).fetchone() + finally: + stop.set() + monitor.join(timeout=max(1, args.progress_interval)) + connection.close() + rows = int(copied[0]) if copied else 0 + files = sorted(incomplete_table.rglob("*.parquet")) + if rows <= 0 or not files: + raise RuntimeError("Variant lookup build produced no rows") + + table_root.parent.mkdir(parents=True, exist_ok=True) + incomplete_table.rename(table_root) + shutil.rmtree(incomplete_root) + elapsed = time.perf_counter() - started + manifest = { + "contract": "association_evidence_v2_partitioned_serving", + "schema_version": "2.8.1-optimized-rc1", + "built_at": utc_now(), + "source_db": str(source_db), + "source_db_sha256": args.source_sha256, + "runtime": { + "python": __import__("platform").python_version(), + "duckdb": duckdb.__version__, + "threads": args.threads, + "memory_limit": args.memory_limit, + "progress_interval": args.progress_interval, + }, + "tables": { + "variants": { + "path": "tables/variants", + "partition_key": "variant_id", + "bucket_function": f"sha256(variant_id)[0:{args.bucket_characters}]", + "bucket_characters": args.bucket_characters, + "directory_layout": "flat_variant_bucket", + "rows": rows, + "files": len(files), + "complete": True, + } + }, + "build": { + "elapsed_seconds": round(elapsed, 3), + "smoke_limit": args.smoke_limit, + "checkpoint_granularity": "atomic complete variants table", + }, + } + write_json(manifest_path, manifest) + + checksum_rows = [] + for path in sorted([manifest_path, *table_root.rglob("*.parquet")]): + checksum_rows.append(f"{sha256_file(path)} {path.relative_to(output_root)}") + checksum_path.write_text("\n".join(checksum_rows) + "\n") + state.update( + { + "complete": True, + "completed_at": utc_now(), + "rows": rows, + "files": len(files), + "elapsed_seconds": round(elapsed, 3), + "checksum_entries": len(checksum_rows), + } + ) + write_json(checkpoint, state) + LOGGER.info( + "variant lookup build complete rows=%d files=%d elapsed_seconds=%.3f", + rows, + len(files), + elapsed, + ) + return manifest + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--source-db", type=Path, required=True) + parser.add_argument("--source-sha256", required=True) + parser.add_argument("--output-root", type=Path, required=True) + parser.add_argument("--bucket-characters", type=int, choices=(2, 3, 4), default=3) + parser.add_argument("--threads", type=int, default=4) + parser.add_argument("--memory-limit", default="4GB") + parser.add_argument("--progress-interval", type=int, default=30) + parser.add_argument("--smoke-limit", type=int) + parser.add_argument("--reset", action="store_true") + parser.add_argument("--verbose", action="store_true") + return parser.parse_args() + + +def main() -> int: + args = parse_args() + logging.basicConfig( + level=logging.DEBUG if args.verbose else logging.INFO, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + ) + result = build(args) + print(json.dumps(result, indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup_duckdb.py b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup_duckdb.py new file mode 100755 index 0000000..e92b1b0 --- /dev/null +++ b/scripts/dataset_specific_scripts/unified/build_association_evidence_v2_variant_lookup_duckdb.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +"""Build an indexed DuckDB lookup from verified canonical-variant Parquet.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import logging +import platform +import threading +import time +from datetime import datetime, timezone +from pathlib import Path + + +LOGGER = logging.getLogger("datahub.association_evidence_v2.variant_lookup_duckdb") + + +def utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def write_json(path: Path, payload: dict) -> None: + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text(json.dumps(payload, indent=2, sort_keys=True) + "\n") + temporary.replace(path) + + +def monitor(path: Path, interval: int, stop: threading.Event) -> None: + started = time.perf_counter() + while not stop.wait(interval): + size = path.stat().st_size if path.exists() else 0 + LOGGER.info( + "variant lookup DuckDB progress elapsed_seconds=%.1f bytes=%d", + time.perf_counter() - started, + size, + ) + + +def build(args: argparse.Namespace) -> dict: + import duckdb + + parquet_root = args.parquet_root.resolve() + output_db = args.output_db.resolve() + checkpoint = output_db.with_suffix(output_db.suffix + ".checkpoint.json") + manifest_path = output_db.with_suffix(output_db.suffix + ".manifest.json") + checksum_path = output_db.with_suffix(output_db.suffix + ".sha256") + temporary_db = output_db.with_suffix(output_db.suffix + ".incomplete") + configuration = { + "parquet_root": str(parquet_root), + "parquet_manifest_sha256": args.parquet_manifest_sha256, + "expected_rows": args.expected_rows, + "smoke_limit": args.smoke_limit, + } + if not parquet_root.is_dir(): + raise FileNotFoundError(f"Canonical variant Parquet root not found: {parquet_root}") + output_db.parent.mkdir(parents=True, exist_ok=True) + if args.reset: + for path in (output_db, checkpoint, manifest_path, checksum_path, temporary_db): + path.unlink(missing_ok=True) + if checkpoint.is_file(): + state = json.loads(checkpoint.read_text()) + if state.get("configuration") != configuration: + raise RuntimeError("Existing checkpoint configuration does not match this run") + if state.get("complete") and output_db.is_file() and manifest_path.is_file(): + LOGGER.info("indexed variant lookup already complete output=%s", output_db) + return json.loads(manifest_path.read_text()) + else: + state = {"configuration": configuration, "started_at": utc_now(), "complete": False} + write_json(checkpoint, state) + temporary_db.unlink(missing_ok=True) + if output_db.exists(): + raise RuntimeError(f"Uncheckpointed lookup database exists: {output_db}") + + glob = str(parquet_root / "**" / "*.parquet").replace("'", "''") + limit_sql = f" LIMIT {int(args.smoke_limit)}" if args.smoke_limit else "" + connection = duckdb.connect( + str(temporary_db), + config={"threads": str(args.threads), "memory_limit": args.memory_limit}, + ) + connection.execute("SET preserve_insertion_order=false") + stop = threading.Event() + progress = threading.Thread( + target=monitor, + args=(temporary_db, args.progress_interval, stop), + daemon=True, + ) + progress.start() + started = time.perf_counter() + try: + connection.execute( + f"CREATE TABLE variants AS SELECT * EXCLUDE (variant_bucket) " + f"FROM read_parquet('{glob}', hive_partitioning=true){limit_sql}" + ) + connection.execute("CREATE UNIQUE INDEX idx_variants_variant_id ON variants(variant_id)") + connection.execute("ANALYZE variants") + rows = int(connection.execute("SELECT count(*) FROM variants").fetchone()[0]) + duplicate_count = int( + connection.execute( + "SELECT count(*)-count(DISTINCT variant_id) FROM variants" + ).fetchone()[0] + ) + connection.execute("CHECKPOINT") + finally: + connection.close() + stop.set() + progress.join(timeout=max(1, args.progress_interval)) + if duplicate_count: + raise RuntimeError(f"Variant lookup contains {duplicate_count} duplicate IDs") + if args.expected_rows is not None and rows != args.expected_rows: + raise RuntimeError( + f"Variant lookup row count mismatch: expected {args.expected_rows}, found {rows}" + ) + temporary_db.rename(output_db) + elapsed = time.perf_counter() - started + digest = sha256_file(output_db) + checksum_path.write_text(f"{digest} {output_db.name}\n") + manifest = { + "contract": "association_evidence_v2_canonical_variant_lookup", + "schema_version": "2.8.1-optimized-rc1", + "built_at": utc_now(), + "database": str(output_db), + "database_sha256": digest, + "source_parquet_root": str(parquet_root), + "source_parquet_manifest_sha256": args.parquet_manifest_sha256, + "table": "variants", + "primary_lookup_key": "variant_id", + "unique_index": "idx_variants_variant_id", + "rows": rows, + "duplicate_variant_ids": duplicate_count, + "runtime": { + "python": platform.python_version(), + "duckdb": duckdb.__version__, + "threads": args.threads, + "memory_limit": args.memory_limit, + }, + "build": { + "elapsed_seconds": round(elapsed, 3), + "smoke_limit": args.smoke_limit, + "checkpoint_granularity": "atomic indexed lookup database", + }, + } + write_json(manifest_path, manifest) + state.update( + { + "complete": True, + "completed_at": utc_now(), + "rows": rows, + "database_sha256": digest, + "elapsed_seconds": round(elapsed, 3), + } + ) + write_json(checkpoint, state) + LOGGER.info( + "indexed variant lookup complete rows=%d bytes=%d elapsed_seconds=%.3f", + rows, + output_db.stat().st_size, + elapsed, + ) + return manifest + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--parquet-root", type=Path, required=True) + parser.add_argument("--parquet-manifest-sha256", required=True) + parser.add_argument("--output-db", type=Path, required=True) + parser.add_argument("--expected-rows", type=int) + parser.add_argument("--threads", type=int, default=2) + parser.add_argument("--memory-limit", default="4GB") + parser.add_argument("--progress-interval", type=int, default=30) + parser.add_argument("--smoke-limit", type=int) + parser.add_argument("--reset", action="store_true") + parser.add_argument("--verbose", action="store_true") + return parser.parse_args() + + +def main() -> int: + args = parse_args() + logging.basicConfig( + level=logging.DEBUG if args.verbose else logging.INFO, + format="%(asctime)s | %(levelname)s | %(name)s | %(message)s", + ) + print(json.dumps(build(args), indent=2, sort_keys=True)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/datahub/association_evidence_v2/__init__.py b/src/datahub/association_evidence_v2/__init__.py new file mode 100644 index 0000000..058cfa7 --- /dev/null +++ b/src/datahub/association_evidence_v2/__init__.py @@ -0,0 +1,5 @@ +"""Scientifically separated association evidence model v2.""" + +from .contracts import SCHEMA_VERSION, SOURCE_DETAIL_UNAVAILABLE_FIELDS + +__all__ = ["SCHEMA_VERSION", "SOURCE_DETAIL_UNAVAILABLE_FIELDS"] diff --git a/src/datahub/association_evidence_v2/builder.py b/src/datahub/association_evidence_v2/builder.py new file mode 100644 index 0000000..89222af --- /dev/null +++ b/src/datahub/association_evidence_v2/builder.py @@ -0,0 +1,2182 @@ +"""Resumable production-scale builder for association evidence model v2.""" + +from __future__ import annotations + +import csv +import gzip +import hashlib +import json +import logging +import os +import platform +import sys +import time +from dataclasses import dataclass +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Iterable + +from datahub.phenotype_paths import PhenotypePathResolver + +from .clinical_classification import ( + CONDITION_STATUS, + EVIDENCE_GRANULARITY, + PROVENANCE_LIMITATION, + display_group_for_term, + normalized_clinical_terms, +) +from .contracts import ( + ASSOCIATION_RECORD_KIND_PROVIDER, + ASSOCIATION_RECORD_KIND_SOURCE_SUMMARY, + POPULATION_FIELDS, + PROVIDER_DETAIL_AVAILABLE, + PROVIDER_DETAIL_NOT_APPLICABLE, + PROVIDER_DETAIL_UNAVAILABLE, + SCHEMA_VERSION, + SOURCE_DETAIL_UNAVAILABLE_FIELDS, +) +from .normalization import normalized_alleles, variation_type_from_alleles +from .schema import create_indexes, create_schema + + +LOGGER = logging.getLogger("datahub.association_evidence_v2") + +RAW_COLUMNS: tuple[tuple[str, str], ...] = ( + ("MarkerID", "marker_id_raw"), + ("pval", "p_value_raw"), + ("GWASsummaryPath", "gwas_summary_path_raw"), + ("studyID", "study_id_raw"), + ("Phenotype", "phenotype_raw"), + ("Study", "study_raw"), + ("PMID", "pmid_raw"), + ("StudyGenomeBuild", "study_genome_build_raw"), + ("dbsnp.rsid", "variant_id_raw"), + ("dbsnp.dbsnp_build", "dbsnp_build_raw"), + ("dbsnp.alleles.allele", "allele_string_raw"), + ("dbsnp.chrom", "chromosome_raw"), + ("dbsnp.hg19.start", "hg19_start_raw"), + ("dbsnp.hg19.end", "hg19_end_raw"), + ("dbsnp.vartype", "source_variation_type_raw"), + ("gnomad_genome.af.af", "frequency_total_raw"), + ("gnomad_genome.af.af_afr", "frequency_afr_raw"), + ("gnomad_genome.af.af_amr", "frequency_amr_raw"), + ("gnomad_genome.af.af_asj", "frequency_asj_raw"), + ("gnomad_genome.af.af_eas", "frequency_eas_raw"), + ("gnomad_genome.af.af_fin", "frequency_fin_raw"), + ("gnomad_genome.af.af_nfe", "frequency_nfe_raw"), + ("gnomad_genome.af.af_oth", "frequency_oth_raw"), + ("snpeff.ann.gene_id", "gene_id_raw"), + ("snpeff.ann.effect", "consequence_raw"), + ("snpeff.ann.putative_impact", "putative_impact_raw"), + ("snpeff.ann.feature_id", "feature_id_raw"), + ("snpeff.ann.hgvs_p", "hgvs_p_raw"), + ("snpeff.ann.protein.length", "protein_length_raw"), + ("dbnsfp.chrom", "dbnsfp_chrom_raw"), + ("dbnsfp.hg18.start", "dbnsfp_hg18_start_raw"), + ("dbnsfp.hg18.end", "dbnsfp_hg18_end_raw"), + ("dbnsfp.hg19.start", "dbnsfp_hg19_start_raw"), + ("dbnsfp.hg19.end", "dbnsfp_hg19_end_raw"), + ("dbnsfp.hg38.start", "dbnsfp_hg38_start_raw"), + ("dbnsfp.hg38.end", "dbnsfp_hg38_end_raw"), + ("dbnsfp.ensembl.proteinid", "ensembl_protein_id_raw"), + ("dbnsfp.ensembl.transcriptid", "ensembl_transcript_id_raw"), + ("clinvar.rcv.clinical_significance", "clinical_significance_raw"), +) + +POPULATION_ALIASES = { + "gnomad_genome.af.af": "frequency_total_raw", + "gnomad_genome.af.af_afr": "frequency_afr_raw", + "gnomad_genome.af.af_amr": "frequency_amr_raw", + "gnomad_genome.af.af_asj": "frequency_asj_raw", + "gnomad_genome.af.af_eas": "frequency_eas_raw", + "gnomad_genome.af.af_fin": "frequency_fin_raw", + "gnomad_genome.af.af_nfe": "frequency_nfe_raw", + "gnomad_genome.af.af_oth": "frequency_oth_raw", +} + + +@dataclass(frozen=True) +class SourceJob: + path: Path + logical_path: str + dataset_type: str + dataset_id: str + source: str + delimiter: str + + @property + def key(self) -> str: + return f"{self.dataset_type}:{self.logical_path}" + + +class BuildCheckpoint: + VERSION = 1 + + def __init__(self, path: Path, configuration_hash: str) -> None: + self.path = path + self.configuration_hash = configuration_hash + self.payload: dict[str, Any] = { + "version": self.VERSION, + "configuration_hash": configuration_hash, + "completed_source_jobs": {}, + "completed_phases": {}, + "updated_at": None, + } + + def load(self) -> None: + if not self.path.exists(): + return + loaded = json.loads(self.path.read_text()) + if loaded.get("version") != self.VERSION: + raise RuntimeError("Checkpoint version does not match this builder") + if loaded.get("configuration_hash") != self.configuration_hash: + raise RuntimeError("Checkpoint configuration does not match this run") + self.payload = loaded + + def save(self) -> None: + self.payload["updated_at"] = _utc_now() + self.path.parent.mkdir(parents=True, exist_ok=True) + temp = self.path.with_suffix(self.path.suffix + ".tmp") + temp.write_text(json.dumps(self.payload, indent=2, sort_keys=True) + "\n") + temp.replace(self.path) + + def source_complete(self, key: str) -> bool: + return key in self.payload["completed_source_jobs"] + + def mark_source_complete(self, job: SourceJob, rows: int, elapsed: float) -> None: + self.payload["completed_source_jobs"][job.key] = { + "rows": rows, + "elapsed_seconds": round(elapsed, 3), + "completed_at": _utc_now(), + } + self.save() + + def phase_complete(self, phase: str) -> bool: + return phase in self.payload["completed_phases"] + + def mark_phase_complete(self, phase: str, details: dict[str, Any]) -> None: + self.payload["completed_phases"][phase] = { + **details, + "completed_at": _utc_now(), + } + self.save() + + +def _utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def _canonical_json(value: Any) -> str: + return json.dumps(value, sort_keys=True, separators=(",", ":"), ensure_ascii=True) + + +def _configuration_hash(value: dict[str, Any]) -> str: + return hashlib.sha256(_canonical_json(value).encode("utf-8")).hexdigest() + + +def _sql_literal(value: str) -> str: + return "'" + value.replace("'", "''") + "'" + + +def _read_header(path: Path, delimiter: str) -> list[str]: + opener = gzip.open if path.name.endswith(".gz") else open + with opener(path, "rt", encoding="utf-8", errors="strict", newline="") as stream: + return next(csv.reader(stream, delimiter=delimiter)) + + +def _input_select(header: set[str]) -> str: + seen_aliases: set[str] = set() + expressions: list[str] = [] + for raw_name, alias in RAW_COLUMNS: + if alias in seen_aliases: + continue + seen_aliases.add(alias) + if raw_name in header: + safe_name = raw_name.replace('"', '""') + expressions.append(f'CAST("{safe_name}" AS VARCHAR) AS {alias}') + else: + expressions.append(f"NULL::VARCHAR AS {alias}") + return ",\n ".join(expressions) + + +def _concat_columns(prefix: str = "") -> str: + seen: set[str] = set() + expressions = [] + for _raw_name, alias in RAW_COLUMNS: + if alias in seen: + continue + seen.add(alias) + expressions.append(f"coalesce({prefix}{alias}, '')") + return ", ".join(expressions) + + +class EvidenceV2Builder: + def __init__( + self, + *, + cvd_root: Path, + trait_root: Path, + variant_index_root: Path, + output_db: Path, + checkpoint_path: Path, + phenotype_tree: Path, + release_id: str, + command: str, + builder_git_commit: str | None = None, + threads: int = 4, + memory_limit: str = "12GB", + temp_directory: Path | None = None, + smoke_files: int | None = None, + smoke_rows_per_file: int | None = None, + provider_chunk_rows: int = 250_000, + progress_interval: int = 1, + ) -> None: + self.cvd_root = cvd_root.resolve() + self.trait_root = trait_root.resolve() + self.variant_index_root = variant_index_root.resolve() + self.output_db = output_db.resolve() + self.phenotype_tree = phenotype_tree.resolve() + self.release_id = release_id + self.command = command + self.builder_git_commit = builder_git_commit + self.threads = max(1, int(threads)) + self.memory_limit = memory_limit + self.temp_directory = temp_directory.resolve() if temp_directory else None + self.smoke_files = smoke_files + self.smoke_rows_per_file = smoke_rows_per_file + self.provider_chunk_rows = int(provider_chunk_rows) + self.progress_interval = max(1, int(progress_interval)) + configuration = { + "schema_version": SCHEMA_VERSION, + "cvd_root": str(self.cvd_root), + "trait_root": str(self.trait_root), + "variant_index_root": str(self.variant_index_root), + "output_db": str(self.output_db), + "phenotype_tree": str(self.phenotype_tree), + "smoke_files": smoke_files, + "smoke_rows_per_file": smoke_rows_per_file, + } + self.checkpoint = BuildCheckpoint( + checkpoint_path.resolve(), _configuration_hash(configuration) + ) + self.connection: Any | None = None + + def run(self, *, resume: bool = True) -> dict[str, Any]: + import duckdb + + started = time.perf_counter() + if resume: + self.checkpoint.load() + self.output_db.parent.mkdir(parents=True, exist_ok=True) + if self.temp_directory: + self.temp_directory.mkdir(parents=True, exist_ok=True) + self.connection = duckdb.connect(str(self.output_db)) + self.connection.execute(f"SET threads={self.threads}") + self.connection.execute(f"SET memory_limit={_sql_literal(self.memory_limit)}") + # Provider source lines come from a single-threaded sequence; retaining global + # insertion order makes wide source-file inserts consume unbounded memory. + self.connection.execute("SET preserve_insertion_order=false") + if self.temp_directory: + self.connection.execute( + f"SET temp_directory={_sql_literal(str(self.temp_directory))}" + ) + create_schema(self.connection) + self._record_build_start() + self._load_phenotype_registry() + self._record_source_completeness() + + jobs = self._source_jobs() + for index, job in enumerate(jobs, start=1): + if resume and self.checkpoint.source_complete(job.key): + continue + file_started = time.perf_counter() + rows = self._ingest_source_file(job) + elapsed = time.perf_counter() - file_started + self.checkpoint.mark_source_complete(job, rows, elapsed) + if index % self.progress_interval == 0 or index == len(jobs): + LOGGER.info( + "Provider ingest progress files=%d/%d rows=%d file=%s elapsed=%.1fs", + index, + len(jobs), + rows, + job.logical_path, + elapsed, + ) + + self._run_phase("provider_relations", self._build_provider_relations) + self._run_phase("source_summaries", self._load_source_summary_records) + self._run_phase("variants", self._build_variants) + self._run_phase("consequences", self._build_consequences) + self._run_phase("clinical_assertions", self._build_clinical_assertions) + self._run_phase("population_observations", self._build_population_observations) + self._run_phase("summaries", self._build_summaries) + self._run_phase("qa", self._run_qa) + self._run_phase("indexes", lambda: create_indexes(self.connection)) + counts = self._counts() + completed_at = _utc_now() + self.connection.execute( + "UPDATE build_metadata SET status='complete', completed_at=?, counts_json=?", + [completed_at, _canonical_json(counts)], + ) + self.connection.execute("CHECKPOINT") + elapsed = time.perf_counter() - started + summary = { + "release_id": self.release_id, + "schema_version": SCHEMA_VERSION, + "output_db": str(self.output_db), + "checkpoint": str(self.checkpoint.path), + "completed_at": completed_at, + "elapsed_seconds": round(elapsed, 3), + "counts": counts, + } + LOGGER.info("Evidence v2 build complete: %s", _canonical_json(summary)) + self.connection.close() + self.connection = None + return summary + + def _record_build_start(self) -> None: + assert self.connection is not None + attempt = { + "started_at": datetime.now(timezone.utc).isoformat(), + "python": sys.version, + "platform": platform.platform(), + "duckdb": self.connection.execute("SELECT version()").fetchone()[0], + "threads": self.threads, + "memory_limit": self.memory_limit, + "temp_directory": str(self.temp_directory) if self.temp_directory else None, + "command": self.command, + } + existing = self.connection.execute( + "SELECT runtime_json FROM build_metadata ORDER BY built_at LIMIT 1" + ).fetchone() + if existing: + runtime = json.loads(str(existing[0])) if existing[0] else {} + history = list(runtime.get("run_attempts") or []) + if not history: + history.append({key: runtime.get(key) for key in ( + "python", "platform", "duckdb", "threads", "memory_limit", + "temp_directory" + )}) + history.append(attempt) + runtime["run_attempts"] = history + runtime["active_settings"] = attempt + self.connection.execute( + "UPDATE build_metadata SET schema_version=?, runtime_json=?, builder_git_commit=?, builder_command=?", + [SCHEMA_VERSION, _canonical_json(runtime), self.builder_git_commit, self.command], + ) + return + runtime = { + **attempt, + "run_attempts": [attempt], + "active_settings": attempt, + } + self.connection.execute( + "INSERT INTO build_metadata VALUES (?, ?, now(), NULL, 'building', ?, ?, ?, ?, ?, NULL, ?)", + [ + SCHEMA_VERSION, + self.release_id, + str(self.cvd_root.parent), + str(self.variant_index_root), + self.builder_git_commit, + self.command, + _canonical_json(runtime), + _canonical_json( + [ + "MVP provider-row lineage is not applicable to compact source-summary associations.", + "Legacy gnomAD frequency fields do not identify REF, ALT, or association effect allele.", + "Legacy SnpEff and ClinVar enrichment versions are not retained in source rows.", + ] + ), + ], + ) + + def _source_jobs(self) -> list[SourceJob]: + jobs = [] + for dataset_type, root, delimiter, dataset_id, source in ( + ("CVD", self.cvd_root, ",", "hbp_legacy_cvd_association", "legacy_cvd_raw"), + ("TRAIT", self.trait_root, "\t", "hbp_legacy_trait_association", "legacy_trait_raw"), + ): + paths = sorted(path for path in root.iterdir() if path.is_file()) + if self.smoke_files is not None: + paths = paths[: self.smoke_files] + for path in paths: + jobs.append( + SourceJob( + path=path, + logical_path=f"raw_data/{root.name}/{path.name}", + dataset_type=dataset_type, + dataset_id=dataset_id, + source=source, + delimiter=delimiter, + ) + ) + return jobs + + def _load_phenotype_registry(self) -> None: + assert self.connection is not None + if self.connection.execute("SELECT count(*) FROM phenotype_registry").fetchone()[0]: + return + resolver = PhenotypePathResolver.from_tree_json(self.phenotype_tree) + rows = [] + for (dataset_type, phenotype), path in sorted(resolver.leaf_paths.items()): + normalized_path = [str(item) for item in path] + rows.append( + ( + dataset_type, + phenotype, + _canonical_json(normalized_path), + " > ".join(normalized_path), + ) + ) + if rows: + self.connection.executemany("INSERT INTO phenotype_registry VALUES (?, ?, ?, ?)", rows) + + def _record_source_completeness(self) -> None: + assert self.connection is not None + self.connection.execute("DELETE FROM source_completeness") + rows = [ + ( + "legacy_cvd_raw", + '["CVD"]', + PROVIDER_DETAIL_AVAILABLE, + None, + _canonical_json(SOURCE_DETAIL_UNAVAILABLE_FIELDS), + None, + _canonical_json([str(self.cvd_root)]), + "GWAS Catalog-derived AWS snapshot; exact upstream release unavailable", + None, + "GWAS Catalog-derived provider rows are recoverable from the archived CVD files.", + ), + ( + "legacy_trait_raw", + '["TRAIT"]', + PROVIDER_DETAIL_AVAILABLE, + None, + _canonical_json(SOURCE_DETAIL_UNAVAILABLE_FIELDS), + None, + _canonical_json([str(self.trait_root)]), + "GWAS Catalog-derived AWS snapshot; exact upstream release unavailable", + None, + "GWAS Catalog-derived provider rows are recoverable from the archived trait files.", + ), + ( + "million_veteran_program", + '["CVD", "TRAIT"]', + PROVIDER_DETAIL_NOT_APPLICABLE, + None, + _canonical_json(SOURCE_DETAIL_UNAVAILABLE_FIELDS), + "variant_index", + _canonical_json([str(self.variant_index_root)]), + None, + None, + "MVP compact variant-phenotype summaries are first-class source-summary associations; provider-row lineage does not apply.", + ), + ] + self.connection.executemany("INSERT INTO source_completeness VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)", rows) + + def _ingest_source_file(self, job: SourceJob) -> int: + assert self.connection is not None + header = set(_read_header(job.path, job.delimiter)) + select_columns = _input_select(header) + raw_concat = _concat_columns("r.") + delimiter = "\\t" if job.delimiter == "\t" else job.delimiter + select_sql_template = f""" +WITH source_rows AS ( + SELECT + nextval('evidence_v2_source_line_seq') AS source_line, + {select_columns} + FROM read_csv( + __HBP_SOURCE_PATH__, + header=true, + all_varchar=true, + delim={_sql_literal(delimiter)}, + quote='"', + strict_mode=false, + ignore_errors=false, + null_padding=true, + parallel=false, + maximum_line_size=8388608 + ) +), +normalized AS ( + SELECT + *, + lower(regexp_replace(replace(trim(coalesce(phenotype_raw, '')), '/', '_'), '\\s+', '_', 'g')) AS phenotype_slug, + sha256(concat_ws(chr(31), {raw_concat})) AS source_row_sha256 + FROM source_rows r +) +SELECT + 'provider:' || sha256(concat_ws(chr(31), {_sql_literal(job.source)}, {_sql_literal(job.logical_path)}, cast(source_line AS VARCHAR), source_row_sha256)) AS provider_record_id, + 'association:' || sha256(concat_ws(chr(31), + {_sql_literal(job.dataset_type)}, {_sql_literal(job.source)}, + coalesce(study_id_raw, ''), coalesce(study_raw, ''), coalesce(pmid_raw, ''), + coalesce(gwas_summary_path_raw, ''), coalesce(study_genome_build_raw, ''), + coalesce(nullif(trim(variant_id_raw), ''), nullif(trim(marker_id_raw), ''), ''), + phenotype_slug, coalesce(p_value_raw, '') + )) AS association_record_id, + {_sql_literal(job.dataset_id)} AS dataset_id, + {_sql_literal(job.dataset_type)} AS dataset_type, + {_sql_literal(job.source)} AS source, + {_sql_literal(job.logical_path)} AS source_file_logical, + {_sql_literal(str(job.path))} AS source_file_absolute, + source_line, + source_row_sha256, + marker_id_raw, p_value_raw, gwas_summary_path_raw, study_id_raw, phenotype_raw, + phenotype_slug, study_raw, pmid_raw, study_genome_build_raw, + coalesce(nullif(trim(variant_id_raw), ''), nullif(trim(marker_id_raw), '')) AS variant_id_raw, + dbsnp_build_raw, allele_string_raw, chromosome_raw, hg19_start_raw, hg19_end_raw, + source_variation_type_raw, frequency_total_raw, frequency_afr_raw, frequency_amr_raw, + frequency_asj_raw, frequency_eas_raw, frequency_fin_raw, frequency_nfe_raw, + frequency_oth_raw, gene_id_raw, consequence_raw, putative_impact_raw, feature_id_raw, + hgvs_p_raw, protein_length_raw, dbnsfp_chrom_raw, dbnsfp_hg18_start_raw, + dbnsfp_hg18_end_raw, dbnsfp_hg19_start_raw, dbnsfp_hg19_end_raw, + dbnsfp_hg38_start_raw, dbnsfp_hg38_end_raw, ensembl_protein_id_raw, + ensembl_transcript_id_raw, clinical_significance_raw, now() AS ingested_at +FROM normalized +""" + # A failed file is restarted from its first source row; completed files remain + # protected by the durable file-level checkpoint. + self.connection.execute("BEGIN") + try: + self.connection.execute( + "DELETE FROM provider_records WHERE dataset_type=? AND source_file_logical=?", + [job.dataset_type, job.logical_path], + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + + total_rows = 0 + for chunk_path, offset, expected_rows in self._iter_source_chunks(job): + insert_sql = select_sql_template.replace( + "__HBP_SOURCE_PATH__", _sql_literal(str(chunk_path)) + ) + self.connection.execute("BEGIN") + try: + # Single-threaded parsing plus a record offset preserves stable source lines + # without retaining a multi-million-row transaction in memory. + self.connection.execute("DROP SEQUENCE IF EXISTS evidence_v2_source_line_seq") + self.connection.execute( + f"CREATE SEQUENCE evidence_v2_source_line_seq START {offset + 2}" + ) + inserted = self.connection.execute( + "INSERT INTO provider_records " + insert_sql + ).fetchone() + rows = int(inserted[0]) if inserted else 0 + if expected_rows is not None and rows != expected_rows: + raise RuntimeError( + f"Provider chunk row mismatch for {job.logical_path}: " + f"expected {expected_rows}, inserted {rows}" + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + total_rows += rows + if total_rows % 1_000_000 == 0: + LOGGER.info( + "Provider file progress file=%s rows=%d", job.logical_path, total_rows + ) + return total_rows + + def _iter_source_chunks(self, job: SourceJob) -> Iterable[tuple[Path, int, int | None]]: + """Yield the original source or bounded CSV-aware chunks in one pass.""" + + oversized_threshold = 2_000_000_000 + if ( + self.provider_chunk_rows <= 0 + and self.smoke_rows_per_file is None + and job.path.stat().st_size <= oversized_threshold + ): + yield job.path, 0, None + return + configured_chunk_rows = self.provider_chunk_rows if self.provider_chunk_rows > 0 else 1_000_000 + chunk_rows = int(self.smoke_rows_per_file or configured_chunk_rows) + chunk_rows = min(chunk_rows, configured_chunk_rows) + chunk_root = (self.temp_directory or self.output_db.parent / "duckdb_tmp") / "source_chunks" + chunk_root.mkdir(parents=True, exist_ok=True) + chunk_key = hashlib.sha256(job.logical_path.encode()).hexdigest()[:16] + csv.field_size_limit(8 * 1024 * 1024) + with job.path.open("r", encoding="utf-8", newline="") as source_handle: + reader = csv.reader( + source_handle, delimiter=job.delimiter, quotechar='"', strict=True + ) + try: + header = next(reader) + except StopIteration: + return + offset = 0 + chunk_index = 0 + while self.smoke_rows_per_file is None or offset < self.smoke_rows_per_file: + path = chunk_root / f"{chunk_key}.{chunk_index:06d}.csv" + row_limit = min( + chunk_rows, + (self.smoke_rows_per_file - offset) + if self.smoke_rows_per_file is not None + else chunk_rows, + ) + written = 0 + with path.open("w", encoding="utf-8", newline="") as chunk_handle: + writer = csv.writer( + chunk_handle, + delimiter=job.delimiter, + quotechar='"', + lineterminator="\n", + ) + writer.writerow(header) + for row in reader: + writer.writerow(row) + written += 1 + if written >= row_limit: + break + if written == 0: + path.unlink(missing_ok=True) + break + try: + yield path, offset, written + finally: + path.unlink(missing_ok=True) + offset += written + chunk_index += 1 + + def _run_phase(self, name: str, callback: Any) -> None: + if self.checkpoint.phase_complete(name): + return + started = time.perf_counter() + LOGGER.info("Evidence v2 phase start: %s", name) + callback() + elapsed = time.perf_counter() - started + self.checkpoint.mark_phase_complete(name, {"elapsed_seconds": round(elapsed, 3)}) + LOGGER.info("Evidence v2 phase complete: %s elapsed=%.1fs", name, elapsed) + + def _build_provider_relations(self) -> None: + assert self.connection is not None + partition_root = ( + self.temp_directory or self.output_db.parent / "duckdb_tmp" + ) / "provider_relation_partitions" + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".provider-relations.json" + ) + provider_count = int( + self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + state: dict[str, Any] = {} + if state_path.exists(): + state = json.loads(state_path.read_text()) + if int(state.get("provider_record_count", -1)) != provider_count: + raise RuntimeError( + "Provider-relation checkpoint does not match provider record count" + ) + + def save_state() -> None: + state["provider_record_count"] = provider_count + state["updated_at"] = _utc_now() + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + + if not state.get("partitions_complete"): + if partition_root.exists(): + stale = partition_root.with_name( + partition_root.name + ".stale-" + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + partition_root.rename(stale) + partition_root.parent.mkdir(parents=True, exist_ok=True) + LOGGER.info( + "Provider relation partitioning start rows=%d root=%s", + provider_count, + partition_root, + ) + self.connection.execute( + f""" +COPY ( + SELECT provider_record_id, association_record_id, dataset_id, dataset_type, + source, source_file_logical, variant_id_raw, phenotype_raw, phenotype_slug, + study_id_raw, study_raw, pmid_raw, gwas_summary_path_raw, + study_genome_build_raw, p_value_raw, gene_id_raw, + source_variation_type_raw, + substr(association_record_id, 13, 2) AS relation_bucket + FROM provider_records + WHERE coalesce(trim(variant_id_raw), '') <> '' + AND coalesce(trim(phenotype_slug), '') <> '' +) TO {_sql_literal(str(partition_root))} +(FORMAT PARQUET, PARTITION_BY (relation_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + parquet_files = sorted(partition_root.rglob("*.parquet")) + if not parquet_files: + raise RuntimeError("Provider relation partitioning produced no Parquet files") + state.update( + { + "partitions_complete": True, + "partition_files": len(parquet_files), + "completed_buckets": [], + } + ) + save_state() + LOGGER.info( + "Provider relation partitioning complete files=%d", len(parquet_files) + ) + + completed = set(state.get("completed_buckets") or []) + bucket_directories = sorted( + path for path in partition_root.iterdir() if path.is_dir() + ) + if not completed: + self.connection.execute( + "CREATE OR REPLACE TABLE association_record_genes AS " + "SELECT * FROM association_record_genes LIMIT 0" + ) + self.connection.execute( + "CREATE OR REPLACE TABLE association_record_provider_records AS " + "SELECT * FROM association_record_provider_records LIMIT 0" + ) + for table in ( + "consequence_annotation_provider_records", + "clinical_assertion_provider_records", + "population_observation_provider_records", + ): + self.connection.execute( + f"CREATE OR REPLACE TABLE {table} AS SELECT * FROM {table} LIMIT 0" + ) + self.connection.execute( + "DELETE FROM association_records WHERE record_kind=?", + [ASSOCIATION_RECORD_KIND_PROVIDER], + ) + + for index, bucket_directory in enumerate(bucket_directories, start=1): + bucket = bucket_directory.name.split("=", 1)[-1] + if bucket in completed: + continue + source_sql = ( + "read_parquet(" + + _sql_literal(str(bucket_directory / "*.parquet")) + + ")" + ) + self.connection.execute("BEGIN") + try: + # Each bucket is committed atomically. A bucket absent from the + # durable checkpoint has no retained rows, so rescanning growing + # output tables for defensive deletes would be quadratic. + self.connection.execute( + f""" +INSERT INTO association_records +SELECT + p.association_record_id, + {_sql_literal(ASSOCIATION_RECORD_KIND_PROVIDER)}, + 'provider_record', + p.source, + arg_min(p.dataset_id, p.provider_record_id), + p.dataset_type, + p.source, + NULL, + p.variant_id_raw, + arg_min(nullif(trim(p.phenotype_raw), ''), p.provider_record_id), + p.phenotype_slug, + coalesce(arg_min(pr.phenotype_path_json, p.provider_record_id), to_json([p.phenotype_slug])), + coalesce(arg_min(pr.phenotype_path_key, p.provider_record_id), p.phenotype_slug), + arg_min(nullif(trim(p.study_id_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.study_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.pmid_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.gwas_summary_path_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.study_genome_build_raw), ''), p.provider_record_id), + min(try_cast(p.p_value_raw AS DOUBLE)), + arg_min(nullif(trim(p.source_variation_type_raw), ''), p.provider_record_id), + NULL, 'unavailable', + NULL, NULL, 'unavailable', + NULL, 'unavailable', + NULL, 'unavailable', + NULL, 'unavailable', + NULL, 'unavailable', + 'available', + NULL, + '["effect_allele","effect_size","standard_error","sample_size","ancestry","fine_mapping"]', + NULL, + NULL, + count(*), + count(DISTINCT p.source_file_logical) +FROM {source_sql} p +LEFT JOIN phenotype_registry pr + ON pr.dataset_type=p.dataset_type AND pr.phenotype_slug=p.phenotype_slug +GROUP BY p.association_record_id, p.dataset_type, p.source, p.variant_id_raw, p.phenotype_slug +""" + ) + self.connection.execute( + f""" +INSERT INTO association_record_provider_records +SELECT association_record_id, provider_record_id FROM {source_sql} +""" + ) + self.connection.execute( + f""" +INSERT INTO association_record_genes +SELECT association_record_id, trim(gene_id_raw), + 'source_reported_annotation_gene', count(*) +FROM {source_sql} +WHERE coalesce(trim(gene_id_raw), '') <> '' + AND regexp_matches(trim(gene_id_raw), '^[A-Za-z][A-Za-z0-9_.-]*$') +GROUP BY association_record_id, trim(gene_id_raw) +""" + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + save_state() + if index % 16 == 0 or index == len(bucket_directories): + LOGGER.info( + "Provider relation progress buckets=%d/%d", + len(completed), + len(bucket_directories), + ) + + def _load_source_summary_records(self) -> None: + """Register MVP source-summary associations without inventing study rows.""" + + assert self.connection is not None + files: list[tuple[str, Path]] = [] + if self.smoke_files is not None: + requested = {"TTN", "PCSK9", "ANK2", "HMGCR", "BMPR2"} + for dataset_type in ("CVD", "TRAIT"): + root = self.variant_index_root / dataset_type + for gene in sorted(requested): + for suffix in (".json.gz", ".json"): + path = root / f"{gene}{suffix}" + if path.exists(): + files.append((dataset_type, path)) + break + else: + for dataset_type in ("CVD", "TRAIT"): + root = self.variant_index_root / dataset_type + if root.exists(): + files.extend( + (dataset_type, Path(entry.path)) + for entry in sorted(os.scandir(root), key=lambda item: item.name) + if entry.is_file() + and (entry.name.endswith(".json") or entry.name.endswith(".json.gz")) + ) + + logical_files = [ + f"variant_index/{dataset_type}/{path.name}" + for dataset_type, path in files + ] + fingerprint = hashlib.sha256(_canonical_json(logical_files).encode()).hexdigest() + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".source-summaries.json" + ) + state: dict[str, Any] = {} + if state_path.exists(): + state = json.loads(state_path.read_text()) + prior_fingerprint = state.get("input_fingerprint") + if prior_fingerprint and prior_fingerprint != fingerprint: + raise RuntimeError("Source-summary checkpoint input mismatch") + + def save_state() -> None: + state["input_fingerprint"] = fingerprint + state.pop("registered_files", None) + state.pop("completed_files", None) + state["completed_file_count"] = len(logical_files) + state["total_file_count"] = len(files) + state["payload_contract_version"] = "gene-scoped-artifact-registry-v1" + state["source_summary_engine_version"] = "artifact-registry-v1" + state["updated_at"] = _utc_now() + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + + missing_fields = _canonical_json(SOURCE_DETAIL_UNAVAILABLE_FIELDS) + reason = None + rows = [] + for dataset_type, path in files: + logical = f"variant_index/{dataset_type}/{path.name}" + gene = path.name.removesuffix(".gz").removesuffix(".json") + stat = path.stat() + rows.append( + ( + "source-summary-artifact:" + + hashlib.sha256( + _canonical_json( + ["million_veteran_program", dataset_type, gene, logical] + ).encode() + ).hexdigest(), + "million_veteran_program", + dataset_type, + gene, + logical, + "json.gz" if path.name.endswith(".json.gz") else "json", + stat.st_size, + stat.st_mtime_ns, + PROVIDER_DETAIL_NOT_APPLICABLE, + reason, + missing_fields, + "retained_compact_variant_index_v1", + "gene_scoped_on_demand", + ) + ) + + self.connection.execute("BEGIN") + try: + # Remove rows left by the superseded per-entry materialization attempt. + self.connection.execute( + "CREATE OR REPLACE TABLE association_records AS " + "SELECT * FROM association_records " + "WHERE record_kind<>'source_summary_association'" + ) + self.connection.execute( + "CREATE OR REPLACE TABLE association_record_genes AS " + "SELECT g.* FROM association_record_genes g " + "JOIN association_records a USING (association_record_id)" + ) + self.connection.execute("DELETE FROM source_summary_artifacts") + if rows: + self.connection.executemany( + "INSERT INTO source_summary_artifacts VALUES " + "(?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)", + rows, + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + + state["superseded_per_entry_rows_removed"] = True + state["registered_artifact_count"] = len(rows) + save_state() + LOGGER.info("MVP source-summary artifact registry complete files=%d", len(rows)) + + def _provider_variant_partitions(self) -> list[Path]: + """Materialize one reusable provider projection partitioned by variant ID.""" + + assert self.connection is not None + root = ( + self.temp_directory or self.output_db.parent / "duckdb_tmp" + ) / "provider_variant_partitions" + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".provider-variants.json" + ) + provider_count = int( + self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + if state_path.exists(): + state = json.loads(state_path.read_text()) + if ( + int(state.get("provider_record_count", -1)) == provider_count + and int(state.get("projection_schema_version", -1)) == 2 + and state.get("partitions_complete") + and root.exists() + ): + return sorted(path for path in root.iterdir() if path.is_dir()) + raise RuntimeError("Provider-variant partition checkpoint is inconsistent") + if root.exists(): + stale = root.with_name( + root.name + ".stale-" + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + root.rename(stale) + root.parent.mkdir(parents=True, exist_ok=True) + LOGGER.info( + "Provider variant partitioning start rows=%d root=%s", provider_count, root + ) + self.connection.execute( + f""" +COPY ( + SELECT provider_record_id, association_record_id, source, variant_id_raw, allele_string_raw, + dbsnp_build_raw, chromosome_raw, hg19_start_raw, hg19_end_raw, + source_variation_type_raw, gene_id_raw, consequence_raw, putative_impact_raw, + feature_id_raw, ensembl_transcript_id_raw, ensembl_protein_id_raw, hgvs_p_raw, + protein_length_raw, study_genome_build_raw, clinical_significance_raw, + frequency_total_raw, frequency_afr_raw, frequency_amr_raw, frequency_asj_raw, + frequency_eas_raw, frequency_fin_raw, frequency_nfe_raw, frequency_oth_raw, + substr(sha256(variant_id_raw), 1, 2) AS variant_bucket + FROM provider_records + WHERE coalesce(trim(variant_id_raw), '') <> '' +) TO {_sql_literal(str(root))} +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + files = sorted(root.rglob("*.parquet")) + if not files: + raise RuntimeError("Provider variant partitioning produced no Parquet files") + state = { + "provider_record_count": provider_count, + "projection_schema_version": 2, + "partitions_complete": True, + "partition_files": len(files), + "completed_at": _utc_now(), + } + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + LOGGER.info("Provider variant partitioning complete files=%d", len(files)) + return sorted(path for path in root.iterdir() if path.is_dir()) + + def _partition_phase_state( + self, phase: str, provider_count: int, *, phase_schema_version: int = 1 + ) -> tuple[Path, dict[str, Any]]: + path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + f".{phase}.json" + ) + state: dict[str, Any] = {} + if path.exists(): + state = json.loads(path.read_text()) + if int(state.get("provider_record_count", -1)) != provider_count: + raise RuntimeError(f"{phase} checkpoint provider count mismatch") + if int(state.get("phase_schema_version", 1)) != phase_schema_version: + state = {} + state["phase_schema_version"] = phase_schema_version + return path, state + + @staticmethod + def _save_partition_phase_state( + path: Path, state: dict[str, Any], provider_count: int + ) -> None: + state["provider_record_count"] = provider_count + state["updated_at"] = _utc_now() + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(path) + + def _build_variants(self) -> None: + assert self.connection is not None + provider_count = int( + self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + partitions = self._provider_variant_partitions() + state_path, state = self._partition_phase_state("variants", provider_count) + completed = set(state.get("completed_buckets") or []) + if not completed: + self.connection.execute("DELETE FROM variants") + variant_sql = """ +INSERT INTO variants +WITH parsed AS ( + SELECT + variant_id_raw, + allele_string_raw, + dbsnp_build_raw, + chromosome_raw, + hg19_start_raw, + hg19_end_raw, + source_variation_type_raw, + regexp_replace( + upper(trim(allele_string_raw)), + '[\\[\\]''"[:space:]]', + '', + 'g' + ) AS clean_alleles + FROM __HBP_VARIANT_SOURCE__ + WHERE coalesce(trim(variant_id_raw), '') <> '' +), +typed AS ( + SELECT + *, + CASE + WHEN regexp_full_match(clean_alleles, '^[ACGTN]+([,/|][ACGTN]+)+$') + THEN regexp_split_to_array(clean_alleles, '[,/|]') + ELSE NULL + END AS normalized_alleles + FROM parsed +), +classified AS ( + SELECT + *, + CASE + WHEN normalized_alleles IS NULL THEN NULL + WHEN list_min(list_transform(normalized_alleles, allele -> length(allele)))=1 + AND list_max(list_transform(normalized_alleles, allele -> length(allele)))=1 + THEN 'SNV' + WHEN list_min(list_transform(normalized_alleles, allele -> length(allele))) + <> list_max(list_transform(normalized_alleles, allele -> length(allele))) + THEN 'INDEL' + ELSE 'MNV' + END AS derived_variation_type, + coalesce(trim(dbsnp_build_raw), '') <> '' + OR coalesce(trim(chromosome_raw), '') <> '' + OR coalesce(trim(hg19_start_raw), '') <> '' + OR coalesce(trim(hg19_end_raw), '') <> '' AS has_coordinate + FROM typed +), +grouped AS ( + SELECT + variant_id_raw AS variant_id, + coalesce( + to_json( + list_sort( + list( + DISTINCT struct_pack( + raw := allele_string_raw, + normalized_unordered := list_sort(normalized_alleles) + ) + ) FILTER (WHERE normalized_alleles IS NOT NULL) + ) + ), + '[]' + ) AS allele_contexts, + coalesce( + to_json( + list_sort( + list( + DISTINCT struct_pack( + dbsnp_build := dbsnp_build_raw, + chromosome := chromosome_raw, + hg19_start := hg19_start_raw, + hg19_end := hg19_end_raw + ) + ) FILTER (WHERE has_coordinate) + ) + ), + '[]' + ) AS coordinate_contexts, + coalesce( + to_json( + list_sort( + list(DISTINCT source_variation_type_raw) + FILTER ( + WHERE coalesce(trim(source_variation_type_raw), '') NOT IN ('', 'NA') + ) + ) + ), + '[]' + ) AS source_variation_types, + count(DISTINCT derived_variation_type) + FILTER (WHERE derived_variation_type IS NOT NULL) AS derived_type_count, + min(derived_variation_type) + FILTER (WHERE derived_variation_type IS NOT NULL) AS only_derived_type, + coalesce( + bool_or(len(normalized_alleles) > 2) + FILTER (WHERE normalized_alleles IS NOT NULL), + false + ) AS multiallelic_observed, + count( + DISTINCT concat_ws( + chr(31), + coalesce(dbsnp_build_raw, ''), + coalesce(chromosome_raw, ''), + coalesce(hg19_start_raw, ''), + coalesce(hg19_end_raw, '') + ) + ) FILTER (WHERE has_coordinate) > 1 AS build_or_position_conflict + FROM classified + GROUP BY variant_id_raw +) +SELECT + variant_id, + allele_contexts, + coordinate_contexts, + source_variation_types, + CASE WHEN derived_type_count=1 THEN only_derived_type ELSE NULL END, + CASE WHEN derived_type_count=1 THEN 'derived' ELSE 'unresolved' END, + 'normalized_unordered_source_allele_lengths', + CASE + WHEN derived_type_count=1 + THEN 'All parseable source allele contexts agree. REF/ALT roles remain unresolved.' + WHEN derived_type_count>1 + THEN 'Source allele contexts imply conflicting variation classes.' + ELSE 'No source allele context supports an unambiguous sequence-length class.' + END, + 'unresolved', + 'unresolved', + 'unresolved', + multiallelic_observed, + build_or_position_conflict +FROM grouped +""" + for index, directory in enumerate(partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + source_sql = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + self.connection.execute("BEGIN") + try: + # Uncheckpointed buckets contain no committed rows. + self.connection.execute( + variant_sql.replace("__HBP_VARIANT_SOURCE__", source_sql) + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + self._save_partition_phase_state(state_path, state, provider_count) + if index % 16 == 0 or index == len(partitions): + LOGGER.info( + "Variant derivation progress buckets=%d/%d variants=%d", + len(completed), len(partitions), + int(self.connection.execute("SELECT count(*) FROM variants").fetchone()[0]), + ) + + LOGGER.info( + "Variant derivation complete total_variants=%d", + int(self.connection.execute("SELECT count(*) FROM variants").fetchone()[0]), + ) + + def _build_consequences(self) -> None: + assert self.connection is not None + provider_count = int(self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0]) + partitions = self._provider_variant_partitions() + state_path, state = self._partition_phase_state("consequences", provider_count) + completed = set(state.get("completed_buckets") or []) + if not completed: + self.connection.execute("DELETE FROM consequence_annotation_provider_records") + self.connection.execute("DELETE FROM consequence_annotations") + annotation_sql = """ +INSERT INTO consequence_annotations +SELECT + 'consequence:' || sha256(concat_ws(chr(31), p.variant_id_raw, coalesce(trim(p.gene_id_raw), ''), + trim(p.consequence_raw), coalesce(trim(p.putative_impact_raw), ''), + coalesce(trim(p.feature_id_raw), ''), coalesce(trim(p.ensembl_transcript_id_raw), ''), + coalesce(trim(p.ensembl_protein_id_raw), ''), coalesce(trim(p.hgvs_p_raw), ''), + coalesce(trim(p.protein_length_raw), ''), coalesce(trim(p.study_genome_build_raw), ''))) AS consequence_annotation_id, + p.variant_id_raw, + nullif(coalesce(trim(p.gene_id_raw), ''), ''), + trim(p.consequence_raw), + arg_min(nullif(trim(p.putative_impact_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.feature_id_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.ensembl_transcript_id_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.ensembl_protein_id_raw), ''), p.provider_record_id), + arg_min(nullif(trim(p.hgvs_p_raw), ''), p.provider_record_id), + arg_min(try_cast(p.protein_length_raw AS DOUBLE), p.provider_record_id), + 'legacy SnpEff-derived source fields', + NULL, + arg_min(nullif(trim(p.study_genome_build_raw), ''), p.provider_record_id), + 'available', + 'not_performed', + count(*) +FROM __HBP_VARIANT_SOURCE__ p +WHERE coalesce(trim(p.variant_id_raw), '') <> '' + AND coalesce(trim(p.consequence_raw), '') NOT IN ('', 'NA') +GROUP BY p.variant_id_raw, coalesce(trim(p.gene_id_raw), ''), trim(p.consequence_raw), + coalesce(trim(p.putative_impact_raw), ''), coalesce(trim(p.feature_id_raw), ''), + coalesce(trim(p.ensembl_transcript_id_raw), ''), coalesce(trim(p.ensembl_protein_id_raw), ''), + coalesce(trim(p.hgvs_p_raw), ''), coalesce(trim(p.protein_length_raw), ''), + coalesce(trim(p.study_genome_build_raw), '') +""" + link_sql = """ +INSERT INTO consequence_annotation_provider_records +SELECT + 'consequence:' || sha256(concat_ws(chr(31), p.variant_id_raw, coalesce(trim(p.gene_id_raw), ''), + trim(p.consequence_raw), coalesce(trim(p.putative_impact_raw), ''), + coalesce(trim(p.feature_id_raw), ''), coalesce(trim(p.ensembl_transcript_id_raw), ''), + coalesce(trim(p.ensembl_protein_id_raw), ''), coalesce(trim(p.hgvs_p_raw), ''), + coalesce(trim(p.protein_length_raw), ''), coalesce(trim(p.study_genome_build_raw), ''))), + p.provider_record_id, + p.association_record_id +FROM __HBP_VARIANT_SOURCE__ p +WHERE coalesce(trim(p.variant_id_raw), '') <> '' + AND coalesce(trim(p.consequence_raw), '') NOT IN ('', 'NA') +""" + for index, directory in enumerate(partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + source_sql = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + self.connection.execute("BEGIN") + try: + # Uncheckpointed buckets contain no committed rows. + self.connection.execute(annotation_sql.replace("__HBP_VARIANT_SOURCE__", source_sql)) + self.connection.execute(link_sql.replace("__HBP_VARIANT_SOURCE__", source_sql)) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + self._save_partition_phase_state(state_path, state, provider_count) + if index % 16 == 0 or index == len(partitions): + LOGGER.info( + "Consequence publication progress buckets=%d/%d annotations=%d", + len(completed), len(partitions), + int(self.connection.execute("SELECT count(*) FROM consequence_annotations").fetchone()[0]), + ) + + def _build_clinical_assertions(self) -> None: + assert self.connection is not None + provider_count = int(self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0]) + partitions = self._provider_variant_partitions() + state_path, state = self._partition_phase_state( + "clinical-assertions", provider_count, phase_schema_version=2 + ) + completed = set(state.get("completed_buckets") or []) + if not completed: + self.connection.execute("DELETE FROM clinical_assertion_provider_records") + self.connection.execute("DELETE FROM clinical_assertions") + dimension_sql = """ +INSERT INTO clinical_assertions +SELECT + 'clinical:' || sha256(concat_ws(chr(31), p.variant_id_raw, + t.normalized_term, p.source)), + p.variant_id_raw, + t.normalized_term, + t.normalized_term, + t.display_group, + p.source, + NULL, + NULL, + ?, + 'available', + ?, + ?, + to_json(list(DISTINCT t.raw_source_value ORDER BY t.raw_source_value)), + count(DISTINCT p.provider_record_id) +FROM __HBP_VARIANT_SOURCE__ p +JOIN __clinical_terms t + ON t.raw_source_value=trim(p.clinical_significance_raw) +WHERE coalesce(trim(p.variant_id_raw), '') <> '' +GROUP BY p.variant_id_raw, t.normalized_term, t.display_group, p.source +""" + link_sql = """ +INSERT INTO clinical_assertion_provider_records +SELECT DISTINCT + 'clinical:' || sha256(concat_ws(chr(31), p.variant_id_raw, + t.normalized_term, p.source)), + p.provider_record_id, + p.association_record_id +FROM __HBP_VARIANT_SOURCE__ p +JOIN __clinical_terms t + ON t.raw_source_value=trim(p.clinical_significance_raw) +WHERE coalesce(trim(p.variant_id_raw), '') <> '' +""" + for index, directory in enumerate(partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + source_sql = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + raw_values = self.connection.execute( + f"""SELECT DISTINCT trim(clinical_significance_raw) FROM {source_sql} + WHERE coalesce(trim(clinical_significance_raw), '') NOT IN ('', 'NA', '[]')""" + ).fetchall() + term_rows = [ + (raw_value, term, display_group_for_term(term)) + for (raw_value,) in raw_values + for term in normalized_clinical_terms(raw_value) + ] + self.connection.execute( + "CREATE OR REPLACE TEMP TABLE __clinical_terms(" + "raw_source_value VARCHAR, normalized_term VARCHAR, display_group VARCHAR)" + ) + if term_rows: + self.connection.executemany( + "INSERT INTO __clinical_terms VALUES (?, ?, ?)", term_rows + ) + self.connection.execute("BEGIN") + try: + # Uncheckpointed buckets contain no committed rows. + if term_rows: + self.connection.execute( + dimension_sql.replace("__HBP_VARIANT_SOURCE__", source_sql), + [CONDITION_STATUS, EVIDENCE_GRANULARITY, PROVENANCE_LIMITATION], + ) + self.connection.execute( + link_sql.replace("__HBP_VARIANT_SOURCE__", source_sql) + ) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + self._save_partition_phase_state(state_path, state, provider_count) + if index % 16 == 0 or index == len(partitions): + LOGGER.info( + "Clinical assertion progress buckets=%d/%d assertions=%d", + len(completed), len(partitions), + int(self.connection.execute( + "SELECT count(*) FROM clinical_assertions" + ).fetchone()[0]), + ) + + def _build_population_observations(self) -> None: + assert self.connection is not None + provider_count = int( + self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + partitions = self._provider_variant_partitions() + state_path, state = self._partition_phase_state( + "population-observations", provider_count + ) + completed = set(state.get("completed_buckets") or []) + if not completed: + self.connection.execute("DELETE FROM population_observation_provider_records") + self.connection.execute("DELETE FROM population_observations") + population_aliases = [ + POPULATION_ALIASES[raw_name] for raw_name, _ in POPULATION_FIELDS + ] + population_labels = " ".join( + f"WHEN {_sql_literal(POPULATION_ALIASES[raw_name])} " + f"THEN {_sql_literal(population)}" + for raw_name, population in POPULATION_FIELDS + ) + identity = "sha256(concat_ws(chr(31), variant_id_raw, population, cast(reported_frequency AS VARCHAR), coalesce(allele_string_raw, ''), coalesce(study_genome_build_raw, '')))" + observation_sql = f""" +INSERT INTO population_observations +SELECT + 'population:' || {identity}, + variant_id_raw, + 'gnomAD genome legacy enrichment', + NULL, + population, + reported_frequency, + to_json(struct_pack(raw_allele_string := allele_string_raw, ref_alt_status := 'unresolved')), + nullif(trim(study_genome_build_raw), ''), + 'source_row_variant_id', + 'unresolved', + 'unresolved', + 'available', + count(*) +FROM __population_long +GROUP BY variant_id_raw, population, reported_frequency, allele_string_raw, study_genome_build_raw +""" + link_sql = f""" +INSERT INTO population_observation_provider_records +SELECT 'population:' || {identity}, provider_record_id, association_record_id +FROM __population_long +""" + for index, directory in enumerate(partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + source_sql = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + # UNPIVOT scans each provider partition once. The previous UNION ALL + # formulation reread it once per population column and was + # prohibitively I/O-bound at production scale. + long_sql = f""" +SELECT + provider_record_id, + association_record_id, + variant_id_raw, + allele_string_raw, + study_genome_build_raw, + CASE population_field {population_labels} END AS population, + try_cast(reported_frequency_raw AS DOUBLE) AS reported_frequency +FROM ( + UNPIVOT ( + SELECT provider_record_id, association_record_id, variant_id_raw, + allele_string_raw, study_genome_build_raw, + {', '.join(population_aliases)} + FROM {source_sql} + ) + ON {', '.join(population_aliases)} + INTO NAME population_field VALUE reported_frequency_raw +) +WHERE try_cast(reported_frequency_raw AS DOUBLE) BETWEEN 0 AND 1 +""" + self.connection.execute( + "CREATE OR REPLACE TEMP TABLE __population_long AS " + long_sql + ) + self.connection.execute("BEGIN") + try: + # Uncheckpointed buckets contain no committed rows. + self.connection.execute(observation_sql) + self.connection.execute(link_sql) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + self._save_partition_phase_state(state_path, state, provider_count) + if index % 16 == 0 or index == len(partitions): + LOGGER.info( + "Population observation progress buckets=%d/%d observations=%d", + len(completed), + len(partitions), + int( + self.connection.execute( + "SELECT count(*) FROM population_observations" + ).fetchone()[0] + ), + ) + + def _published_variant_partitions(self, table: str) -> dict[str, Path]: + """Partition one derived table once for bounded summary joins.""" + + assert self.connection is not None + allowed = { + "consequence_annotations", + "clinical_assertions", + "population_observations", + } + if table not in allowed: + raise ValueError(f"Unsupported published variant table: {table}") + row_count = int( + self.connection.execute(f"SELECT count(*) FROM {table}").fetchone()[0] + ) + root = ( + self.temp_directory or self.output_db.parent / "duckdb_tmp" + ) / "published_variant_partitions" / table + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + f".published-{table}.json" + ) + if state_path.exists(): + state = json.loads(state_path.read_text()) + directories = { + path.name.split("=", 1)[-1]: path + for path in root.iterdir() + if path.is_dir() + } if root.exists() else {} + if ( + int(state.get("row_count", -1)) == row_count + and state.get("partitions_complete") + and (row_count == 0 or directories) + ): + return directories + + if root.exists(): + stale = root.with_name( + root.name + + ".stale-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + root.rename(stale) + root.parent.mkdir(parents=True, exist_ok=True) + if row_count: + LOGGER.info( + "Published variant partitioning start table=%s rows=%d", + table, + row_count, + ) + self.connection.execute( + f""" +COPY ( + SELECT *, substr(sha256(variant_id), 1, 2) AS variant_bucket + FROM {table} + WHERE coalesce(trim(variant_id), '') <> '' +) TO {_sql_literal(str(root))} +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + directories = { + path.name.split("=", 1)[-1]: path + for path in root.iterdir() + if path.is_dir() + } if root.exists() else {} + state = { + "table": table, + "row_count": row_count, + "partitions_complete": True, + "partition_files": len(list(root.rglob("*.parquet"))) if root.exists() else 0, + "completed_at": _utc_now(), + } + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + LOGGER.info( + "Published variant partitioning complete table=%s buckets=%d", + table, + len(directories), + ) + return directories + + def _summary_base_partitions(self) -> list[Path]: + """Materialize a resumable variant-bucketed summary input projection.""" + + assert self.connection is not None + counts = { + "association_records": int( + self.connection.execute("SELECT count(*) FROM association_records").fetchone()[0] + ), + "association_record_genes": int( + self.connection.execute("SELECT count(*) FROM association_record_genes").fetchone()[0] + ), + "variants": int( + self.connection.execute("SELECT count(*) FROM variants").fetchone()[0] + ), + } + root = ( + self.temp_directory or self.output_db.parent / "duckdb_tmp" + ) / "summary_base_partitions" + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".summary-base.json" + ) + if state_path.exists(): + state = json.loads(state_path.read_text()) + directories = sorted(path for path in root.iterdir() if path.is_dir()) if root.exists() else [] + if ( + state.get("input_counts") == counts + and state.get("partitions_complete") + and directories + ): + return directories + + if root.exists(): + stale = root.with_name( + root.name + + ".stale-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + root.rename(stale) + root.parent.mkdir(parents=True, exist_ok=True) + LOGGER.info("Summary base partitioning start inputs=%s", counts) + self.connection.execute( + f""" +COPY ( + SELECT + 'summary:' || sha256(concat_ws(chr(31), ar.dataset_type, arg.gene_id, ar.variant_id, ar.phenotype_path_key)) AS summary_id, + ar.dataset_type, + arg.gene_id, + ar.variant_id, + ar.phenotype_slug, + ar.phenotype_path_json, + ar.phenotype_path_key, + ar.association_record_id, + ar.source, + ar.provider_detail_status, + ar.provider_record_count, + ar.study_id, + ar.study_title, + ar.reported_p_value, + v.variation_type, + coalesce(v.variation_type_status, 'unresolved') AS variation_type_status, + substr(sha256(ar.variant_id), 1, 2) AS variant_bucket + FROM association_records ar + JOIN association_record_genes arg USING (association_record_id) + LEFT JOIN variants v USING (variant_id) +) TO {_sql_literal(str(root))} +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + directories = sorted(path for path in root.iterdir() if path.is_dir()) + if not directories: + raise RuntimeError("Summary base partitioning produced no Parquet files") + state = { + "input_counts": counts, + "partitions_complete": True, + "partition_files": len(list(root.rglob("*.parquet"))), + "completed_at": _utc_now(), + } + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + LOGGER.info("Summary base partitioning complete buckets=%d", len(directories)) + return directories + + def _build_summaries(self) -> None: + assert self.connection is not None + provider_count = int( + self.connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + summary_partitions = self._summary_base_partitions() + consequence_partitions = self._published_variant_partitions( + "consequence_annotations" + ) + clinical_partitions = self._published_variant_partitions( + "clinical_assertions" + ) + population_partitions = self._published_variant_partitions( + "population_observations" + ) + state_path, state = self._partition_phase_state("summaries", provider_count) + completed = set(state.get("completed_buckets") or []) + if not completed: + for table in ( + "summary_association_records", + "summary_consequence_annotations", + "summary_clinical_assertions", + "variant_phenotype_summaries", + ): + self.connection.execute(f"DELETE FROM {table}") + + summary_sql = """ +INSERT INTO variant_phenotype_summaries +WITH grouped AS ( + SELECT + summary_id, + first(dataset_type) AS dataset_type, + first(gene_id) AS gene_id, + first(variant_id) AS variant_id, + first(phenotype_slug) AS phenotype_slug, + first(phenotype_path_json) AS phenotype_path_json, + first(phenotype_path_key) AS phenotype_path_key, + sum(provider_record_count) AS source_observation_count, + count(DISTINCT association_record_id) AS association_record_count, + count(DISTINCT CASE + WHEN provider_detail_status='available' THEN association_record_id + END) AS provider_association_record_count, + count(DISTINCT CASE + WHEN provider_detail_status='unavailable' THEN association_record_id + END) AS retained_source_summary_count, + count(DISTINCT CASE + WHEN provider_detail_status='available' + AND coalesce(study_id, study_title) IS NOT NULL + THEN coalesce(study_id, study_title) + END) AS study_count, + min(reported_p_value) AS minimum_p, + first(variation_type) AS variation_type, + first(variation_type_status) AS variation_type_status, + to_json( + list_sort( + list( + DISTINCT struct_pack( + source := source, + provider_detail_status := provider_detail_status + ) + ) + ) + ) AS source_completeness + FROM __summary_base + GROUP BY summary_id +), +annotation_counts AS ( + SELECT + sb.summary_id, + count(DISTINCT ca.consequence_annotation_id) AS consequence_count, + count(DISTINCT ca.consequence) > 1 AS consequence_conflict + FROM __summary_base sb + JOIN __bucket_consequences ca + ON ca.variant_id=sb.variant_id AND ca.gene_id=sb.gene_id + GROUP BY sb.summary_id +), +clinical_counts AS ( + SELECT + sb.summary_id, + count(DISTINCT ca.clinical_assertion_id) AS clinical_count, + count(DISTINCT ca.clinical_significance) > 1 AS clinical_conflict + FROM __summary_base sb + JOIN __bucket_clinical ca ON ca.variant_id=sb.variant_id + GROUP BY sb.summary_id +), +population_counts AS ( + SELECT variant_id, count(*) AS population_count + FROM __bucket_population + GROUP BY variant_id +) +SELECT + g.summary_id, + g.dataset_type, + g.gene_id, + g.variant_id, + g.phenotype_slug, + g.phenotype_path_json, + g.phenotype_path_key, + g.source_observation_count, + g.association_record_count, + g.provider_association_record_count, + g.retained_source_summary_count, + g.study_count, + coalesce(pc.population_count, 0), + coalesce(ac.consequence_count, 0), + coalesce(cc.clinical_count, 0), + g.minimum_p, + ( + SELECT min(supplier.association_record_id) + FROM __summary_base supplier + WHERE supplier.summary_id=g.summary_id + AND supplier.reported_p_value=g.minimum_p + ), + coalesce( + ( + SELECT to_json(list(association_record_id ORDER BY association_record_id)) + FROM __summary_base tied + WHERE tied.summary_id=g.summary_id + AND tied.reported_p_value=g.minimum_p + ), + '[]' + ), + NULL, + 'unavailable', + NULL, + 'unavailable', + NULL, + 'unavailable', + coalesce(ac.consequence_conflict, false), + coalesce(cc.clinical_conflict, false), + g.variation_type, + g.variation_type_status, + g.source_completeness +FROM grouped g +LEFT JOIN annotation_counts ac ON ac.summary_id=g.summary_id +LEFT JOIN clinical_counts cc ON cc.summary_id=g.summary_id +LEFT JOIN population_counts pc ON pc.variant_id=g.variant_id +""" + for index, directory in enumerate(summary_partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + summary_source = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + consequence_directory = consequence_partitions.get(bucket) + clinical_directory = clinical_partitions.get(bucket) + population_directory = population_partitions.get(bucket) + consequence_source = ( + "read_parquet(" + _sql_literal(str(consequence_directory / "*.parquet")) + ")" + if consequence_directory + else "(SELECT * FROM consequence_annotations WHERE false)" + ) + clinical_source = ( + "read_parquet(" + _sql_literal(str(clinical_directory / "*.parquet")) + ")" + if clinical_directory + else "(SELECT * FROM clinical_assertions WHERE false)" + ) + population_source = ( + "read_parquet(" + _sql_literal(str(population_directory / "*.parquet")) + ")" + if population_directory + else "(SELECT * FROM population_observations WHERE false)" + ) + self.connection.execute( + f"CREATE OR REPLACE TEMP VIEW __summary_base AS SELECT * FROM {summary_source}" + ) + self.connection.execute( + f"CREATE OR REPLACE TEMP VIEW __bucket_consequences AS SELECT * FROM {consequence_source}" + ) + self.connection.execute( + f"CREATE OR REPLACE TEMP VIEW __bucket_clinical AS SELECT * FROM {clinical_source}" + ) + self.connection.execute( + f"CREATE OR REPLACE TEMP VIEW __bucket_population AS SELECT * FROM {population_source}" + ) + self.connection.execute("BEGIN") + try: + # Uncheckpointed buckets contain no committed rows. + self.connection.execute( + """ +INSERT INTO summary_association_records +SELECT DISTINCT summary_id, association_record_id +FROM __summary_base +""" + ) + self.connection.execute( + """ +INSERT INTO summary_consequence_annotations +SELECT DISTINCT sb.summary_id, ca.consequence_annotation_id +FROM __summary_base sb +JOIN __bucket_consequences ca + ON ca.variant_id=sb.variant_id AND ca.gene_id=sb.gene_id +""" + ) + self.connection.execute( + """ +INSERT INTO summary_clinical_assertions +SELECT DISTINCT sb.summary_id, ca.clinical_assertion_id +FROM __summary_base sb +JOIN __bucket_clinical ca ON ca.variant_id=sb.variant_id +""" + ) + self.connection.execute(summary_sql) + self.connection.execute("COMMIT") + except Exception: + try: + self.connection.execute("ROLLBACK") + except Exception: + pass + raise + completed.add(bucket) + state["completed_buckets"] = sorted(completed) + self._save_partition_phase_state(state_path, state, provider_count) + if index % 16 == 0 or index == len(summary_partitions): + LOGGER.info( + "Summary publication progress buckets=%d/%d summaries=%d", + len(completed), + len(summary_partitions), + int( + self.connection.execute( + "SELECT count(*) FROM variant_phenotype_summaries" + ).fetchone()[0] + ), + ) + + def _minimum_p_summary_partitions(self) -> dict[str, Path]: + """Partition persisted minimum-p suppliers once for bounded QA.""" + + assert self.connection is not None + row_count = int( + self.connection.execute( + """ +SELECT count(*) +FROM variant_phenotype_summaries +WHERE minimum_reported_p_value IS NOT NULL +""" + ).fetchone()[0] + ) + root = ( + self.temp_directory or self.output_db.parent / "duckdb_tmp" + ) / "minimum_p_summary_partitions" + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".minimum-p-summary.json" + ) + if state_path.exists(): + state = json.loads(state_path.read_text()) + directories = { + path.name.split("=", 1)[-1]: path + for path in root.iterdir() + if path.is_dir() + } if root.exists() else {} + if ( + int(state.get("row_count", -1)) == row_count + and state.get("partitions_complete") + and (row_count == 0 or directories) + ): + return directories + + if root.exists(): + stale = root.with_name( + root.name + + ".stale-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + root.rename(stale) + root.parent.mkdir(parents=True, exist_ok=True) + if row_count: + LOGGER.info( + "Minimum-p summary QA partitioning start rows=%d", row_count + ) + self.connection.execute( + f""" +COPY ( + SELECT + variant_phenotype_summary_id, + minimum_reported_p_value, + minimum_reported_p_value_association_record_id, + substr(sha256(variant_id), 1, 2) AS variant_bucket + FROM variant_phenotype_summaries + WHERE minimum_reported_p_value IS NOT NULL +) TO {_sql_literal(str(root))} +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + directories = { + path.name.split("=", 1)[-1]: path + for path in root.iterdir() + if path.is_dir() + } if root.exists() else {} + state = { + "row_count": row_count, + "partitions_complete": True, + "partition_files": len(list(root.rglob("*.parquet"))) if root.exists() else 0, + "completed_at": _utc_now(), + } + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + LOGGER.info( + "Minimum-p summary QA partitioning complete buckets=%d", + len(directories), + ) + return directories + + def _qa_minimum_p_supplier_is_deterministic_tied_minimum(self) -> bool: + """Reconcile tied minimum-p suppliers one variant bucket at a time.""" + + assert self.connection is not None + base_partitions = self._summary_base_partitions() + persisted_partitions = self._minimum_p_summary_partitions() + state_path = self.checkpoint.path.with_name( + self.checkpoint.path.stem + ".qa-minimum-p-supplier.json" + ) + state = json.loads(state_path.read_text()) if state_path.exists() else {} + completed = set(state.get("completed_buckets") or []) + for index, directory in enumerate(base_partitions, start=1): + bucket = directory.name.split("=", 1)[-1] + if bucket in completed: + continue + persisted_directory = persisted_partitions.get(bucket) + persisted_source = ( + "read_parquet(" + _sql_literal(str(persisted_directory / "*.parquet")) + ")" + if persisted_directory + else "(SELECT * FROM variant_phenotype_summaries WHERE false)" + ) + base_source = "read_parquet(" + _sql_literal(str(directory / "*.parquet")) + ")" + passed = bool( + self.connection.execute( + f""" +WITH minima AS ( + SELECT summary_id, min(reported_p_value) AS minimum_p + FROM {base_source} + WHERE reported_p_value IS NOT NULL + GROUP BY summary_id +), expected AS ( + SELECT b.summary_id, m.minimum_p, + min(b.association_record_id) AS supplier_id + FROM {base_source} b + JOIN minima m + ON m.summary_id=b.summary_id + AND b.reported_p_value=m.minimum_p + GROUP BY b.summary_id, m.minimum_p +), persisted AS ( + SELECT * FROM {persisted_source} +) +SELECT count(*)=0 +FROM expected e +FULL OUTER JOIN persisted p + ON p.variant_phenotype_summary_id=e.summary_id +WHERE e.summary_id IS NULL + OR p.variant_phenotype_summary_id IS NULL + OR p.minimum_reported_p_value IS DISTINCT FROM e.minimum_p + OR p.minimum_reported_p_value_association_record_id IS DISTINCT FROM e.supplier_id +""" + ).fetchone()[0] + ) + if not passed: + return False + completed.add(bucket) + state = { + "completed_buckets": sorted(completed), + "total_buckets": len(base_partitions), + "updated_at": _utc_now(), + } + temporary = state_path.with_suffix(state_path.suffix + ".tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(state_path) + if index % 16 == 0 or index == len(base_partitions): + LOGGER.info( + "Minimum-p supplier QA progress buckets=%d/%d", + len(completed), + len(base_partitions), + ) + return True + + def _run_qa(self) -> None: + assert self.connection is not None + checks = { + "provider_record_ids_are_unique": "SELECT count(*)=count(DISTINCT provider_record_id) FROM provider_records", + "provider_source_lines_are_unique": "SELECT count(*)=0 FROM (SELECT source_file_logical, source_line FROM provider_records GROUP BY ALL HAVING count(*)<>1) duplicates", + "association_provider_links_resolve": "SELECT count(*)=0 FROM association_record_provider_records l LEFT JOIN provider_records p USING(provider_record_id) WHERE p.provider_record_id IS NULL", + "association_links_resolve": "SELECT count(*)=0 FROM association_record_provider_records l LEFT JOIN association_records a USING(association_record_id) WHERE a.association_record_id IS NULL", + "consequence_provider_links_resolve": "SELECT count(*)=0 FROM consequence_annotation_provider_records l LEFT JOIN provider_records p USING(provider_record_id) WHERE p.provider_record_id IS NULL", + "clinical_provider_links_resolve": "SELECT count(*)=0 FROM clinical_assertion_provider_records l LEFT JOIN provider_records p USING(provider_record_id) WHERE p.provider_record_id IS NULL", + "population_provider_links_resolve": "SELECT count(*)=0 FROM population_observation_provider_records l LEFT JOIN provider_records p USING(provider_record_id) WHERE p.provider_record_id IS NULL", + "population_observation_integrity": """ + WITH observations AS ( + SELECT population_observation_id, count(*) AS row_count, + min(provider_record_count) AS expected_links, + max(provider_record_count) AS maximum_expected_links + FROM population_observations + GROUP BY population_observation_id + ), + links AS ( + SELECT population_observation_id, count(*) AS link_count + FROM population_observation_provider_records + GROUP BY population_observation_id + ) + SELECT count(*)=0 + FROM observations o + FULL OUTER JOIN links l USING(population_observation_id) + WHERE coalesce(o.row_count, 0)<>1 + OR o.expected_links<>o.maximum_expected_links + OR coalesce(l.link_count, 0)<>coalesce(o.expected_links, 0) + """, + "summary_minimum_p_reconciles": "SELECT count(*)=0 FROM variant_phenotype_summaries s LEFT JOIN association_records a ON a.association_record_id=s.minimum_reported_p_value_association_record_id WHERE s.minimum_reported_p_value IS NOT NULL AND (a.reported_p_value IS NULL OR a.reported_p_value<>s.minimum_reported_p_value)", + "summary_associations_resolve": "SELECT count(*)=0 FROM summary_association_records l LEFT JOIN association_records a USING(association_record_id) WHERE a.association_record_id IS NULL", + "summary_ids_are_unique": "SELECT count(*)=count(DISTINCT variant_phenotype_summary_id) FROM variant_phenotype_summaries", + "summary_association_links_are_unique": "SELECT count(*)=count(DISTINCT variant_phenotype_summary_id || chr(31) || association_record_id) FROM summary_association_records", + "summary_consequence_links_are_unique": "SELECT count(*)=count(DISTINCT variant_phenotype_summary_id || chr(31) || consequence_annotation_id) FROM summary_consequence_annotations", + "summary_clinical_links_are_unique": "SELECT count(*)=count(DISTINCT variant_phenotype_summary_id || chr(31) || clinical_assertion_id) FROM summary_clinical_assertions", + "not_applicable_sources_have_no_provider_links": "SELECT count(*)=0 FROM association_records a JOIN association_record_provider_records l USING(association_record_id) WHERE a.provider_detail_status='not_applicable'", + "all_association_variants_resolve": "SELECT count(*)=0 FROM association_records a LEFT JOIN variants v ON v.variant_id=a.variant_id WHERE v.variant_id IS NULL", + "source_summary_paths_are_logical": "SELECT count(*)=0 FROM source_summary_artifacts WHERE logical_artifact LIKE '/%'", + "source_summary_registry_is_unique": "SELECT count(*)=count(DISTINCT source || chr(31) || dataset_type || chr(31) || gene_id || chr(31) || logical_artifact) FROM source_summary_artifacts", + "source_summary_registry_is_first_class": "SELECT count(*)=0 FROM source_summary_artifacts WHERE provider_detail_status<>'not_applicable' OR publication_mode<>'gene_scoped_on_demand'", + "no_fake_mvp_provider_rows": "SELECT count(*)=0 FROM association_records WHERE source='million_veteran_program' AND record_kind<>'source_summary_association'", + "clinical_assertion_links_resolve": "SELECT count(*)=0 FROM clinical_assertion_provider_records l LEFT JOIN clinical_assertions c USING(clinical_assertion_id) WHERE c.clinical_assertion_id IS NULL", + "clinical_assertion_keys_are_unique": "SELECT count(*)=count(DISTINCT variant_id || chr(31) || normalized_term || chr(31) || assertion_source) FROM clinical_assertions", + "clinical_assertions_use_approved_terms": "SELECT count(*)=0 FROM clinical_assertions WHERE normalized_term NOT IN ('Benign', 'Likely benign', 'Benign/Likely benign', 'Pathogenic', 'Likely pathogenic', 'Pathogenic/Likely pathogenic', 'Uncertain significance', 'Conflicting interpretations of pathogenicity', 'association', 'risk factor', 'protective', 'drug response', 'Affects', 'other', 'not provided')", + "clinical_assertion_labels_are_not_arrays": "SELECT count(*)=0 FROM clinical_assertions WHERE regexp_matches(normalized_term, '^\\s*\\[')", + "clinical_assertions_have_provider_support": "SELECT count(*)=0 FROM clinical_assertions WHERE provider_record_count=0", + "available_associations_have_provider_links": "SELECT count(*)=0 FROM association_records a WHERE a.provider_detail_status='available' AND NOT EXISTS (SELECT 1 FROM association_record_provider_records l WHERE l.association_record_id=a.association_record_id)", + "consequence_link_counts_reconcile": """ + WITH links AS ( + SELECT consequence_annotation_id, count(*) AS link_count + FROM consequence_annotation_provider_records + GROUP BY consequence_annotation_id + ) + SELECT count(*)=0 + FROM consequence_annotations c + LEFT JOIN links l USING(consequence_annotation_id) + WHERE c.provider_record_count<>coalesce(l.link_count, 0) + """, + "clinical_link_counts_reconcile": """ + WITH links AS ( + SELECT clinical_assertion_id, count(*) AS link_count + FROM clinical_assertion_provider_records + GROUP BY clinical_assertion_id + ) + SELECT count(*)=0 + FROM clinical_assertions c + LEFT JOIN links l USING(clinical_assertion_id) + WHERE c.provider_record_count<>coalesce(l.link_count, 0) + """, + "minimum_p_supplier_is_deterministic_tied_minimum": None, + } + completed = { + str(row[0]) + for row in self.connection.execute( + "SELECT check_name FROM qa_results WHERE status='passed'" + ).fetchall() + } + for name, query in checks.items(): + if name in completed: + LOGGER.info("Evidence v2 QA check already passed: %s", name) + continue + started = time.perf_counter() + LOGGER.info("Evidence v2 QA check start: %s", name) + if name == "minimum_p_supplier_is_deterministic_tied_minimum": + passed = ( + self._qa_minimum_p_supplier_is_deterministic_tied_minimum() + ) + else: + assert query is not None + passed = bool(self.connection.execute(query).fetchone()[0]) + row = ( + name, + "passed" if passed else "failed", + str(passed), + "true", + "{}", + datetime.now(timezone.utc), + ) + self.connection.execute( + "DELETE FROM qa_results WHERE check_name=?", [name] + ) + self.connection.execute( + "INSERT INTO qa_results VALUES (?, ?, ?, ?, ?, ?)", row + ) + LOGGER.info( + "Evidence v2 QA check complete: %s passed=%s elapsed=%.1fs", + name, + passed, + time.perf_counter() - started, + ) + if not passed: + raise RuntimeError(f"Evidence v2 QA failed: {name}") + + def _counts(self) -> dict[str, int]: + assert self.connection is not None + tables = ( + "provider_records", + "source_summary_artifacts", + "variants", + "association_records", + "association_record_genes", + "consequence_annotations", + "clinical_assertions", + "population_observations", + "variant_phenotype_summaries", + ) + return { + table: int(self.connection.execute(f"SELECT count(*) FROM {table}").fetchone()[0]) + for table in tables + } diff --git a/src/datahub/association_evidence_v2/clinical_classification.py b/src/datahub/association_evidence_v2/clinical_classification.py new file mode 100644 index 0000000..8ae7705 --- /dev/null +++ b/src/datahub/association_evidence_v2/clinical_classification.py @@ -0,0 +1,172 @@ +"""Audited normalization for legacy ClinVar-derived classification values.""" + +from __future__ import annotations + +import ast +import json +import re +import unicodedata +from collections import Counter +from dataclasses import dataclass +from typing import Any + + +PROVENANCE_LIMITATION = ( + "Legacy ClinVar-derived variant-level classification. The retained source " + "does not include ClinVar condition, accession, review status, classification " + "date, or release version and may not refer to the selected HBP phenotype." +) +EVIDENCE_GRANULARITY = "legacy_variant_level_annotation" +CONDITION_STATUS = "unavailable" + +TERM_TO_DISPLAY_GROUP = { + "Benign": "Benign spectrum", + "Likely benign": "Benign spectrum", + "Benign/Likely benign": "Benign spectrum", + "Pathogenic": "Pathogenic spectrum", + "Likely pathogenic": "Pathogenic spectrum", + "Pathogenic/Likely pathogenic": "Pathogenic spectrum", + "Uncertain significance": "Uncertain or conflicting", + "Conflicting interpretations of pathogenicity": "Uncertain or conflicting", + "association": "Association and risk context", + "risk factor": "Association and risk context", + "protective": "Association and risk context", + "drug response": "Drug response", + "Affects": "Other phenotype effect", + "other": "Details only", + "not provided": "Details only", +} +APPROVED_TERMS = tuple(TERM_TO_DISPLAY_GROUP) +MISSING_MARKERS = {"", "na", "n/a", "null", "none", "[]"} + + +def _normalized_text(value: Any) -> str: + text = unicodedata.normalize("NFKC", str(value or "")) + text = text.replace("\u2013", "-").replace("\u2014", "-") + text = text.strip().strip("\"'").strip() + text = re.sub(r"\s*/\s*", "/", text) + text = re.sub(r"\s+", " ", text) + return text + + +def _term_key(value: Any) -> str: + text = _normalized_text(value).casefold() + text = text.replace("_", " ") + text = re.sub(r"(?<=\w)-(?=\w)", " ", text) + text = re.sub(r"[.]$", "", text) + return re.sub(r"\s+", " ", text).strip() + + +TERM_ALIASES = {_term_key(term): term for term in APPROVED_TERMS} + + +@dataclass(frozen=True) +class ParsedClinicalClassification: + raw_value: str + parse_form: str + terms: tuple[str, ...] + display_terms: tuple[str, ...] + duplicate_terms: dict[str, int] + unparsed_terms: tuple[str, ...] + + @property + def meaningful(self) -> bool: + return bool(self.terms) + + @property + def confident(self) -> bool: + return self.meaningful and not self.unparsed_terms + + +def _canonical_term(value: Any) -> str | None: + return TERM_ALIASES.get(_term_key(value)) + + +def _split_values(text: str) -> tuple[list[Any], str]: + if text.startswith("[") and text.endswith("]"): + try: + value = json.loads(text) + if isinstance(value, list): + return value, "json_array" + except (TypeError, ValueError, json.JSONDecodeError): + pass + try: + value = ast.literal_eval(text) + if isinstance(value, (list, tuple)): + return list(value), "python_list" + except (TypeError, ValueError, SyntaxError): + pass + return re.split(r"\s*[,;]\s*", text[1:-1]), "unquoted_bracketed_list" + + if _canonical_term(text): + return [text], "single" + if ";" in text: + return re.split(r"\s*;\s*", text), "semicolon_delimited" + if "," in text: + return re.split(r"\s*,\s*", text), "comma_delimited" + return [text], "single" + + +def parse_clinical_classification(value: Any) -> ParsedClinicalClassification: + """Parse one raw value without modifying it or choosing a representative term.""" + + raw_value = "" if value is None else str(value) + text = raw_value.strip() + if text.casefold() in MISSING_MARKERS: + return ParsedClinicalClassification( + raw_value=raw_value, + parse_form="missing_placeholder", + terms=(), + display_terms=(), + duplicate_terms={}, + unparsed_terms=(), + ) + + values, parse_form = _split_values(text) + terms: list[str] = [] + unparsed: list[str] = [] + pending = list(values) + while pending: + item = pending.pop(0) + normalized = _normalized_text(item) + if not normalized: + continue + canonical = _canonical_term(normalized) + if canonical: + terms.append(canonical) + continue + if ";" in normalized: + pending[0:0] = re.split(r"\s*;\s*", normalized) + continue + if "," in normalized: + pending[0:0] = re.split(r"\s*,\s*", normalized) + continue + unparsed.append(normalized) + + counts = Counter(terms) + return ParsedClinicalClassification( + raw_value=raw_value, + parse_form=parse_form, + terms=tuple(terms), + display_terms=tuple(dict.fromkeys(terms)), + duplicate_terms={term: count for term, count in counts.items() if count > 1}, + unparsed_terms=tuple(dict.fromkeys(unparsed)), + ) + + +def normalized_clinical_terms(value: Any, *, strict: bool = True) -> list[str]: + parsed = parse_clinical_classification(value) + if strict and parsed.unparsed_terms: + raise ValueError( + "Unparsed clinical-classification term(s): " + + "; ".join(parsed.unparsed_terms) + + f"; raw value={parsed.raw_value!r}" + ) + return list(parsed.display_terms) + + +def display_group_for_term(term: str) -> str: + try: + return TERM_TO_DISPLAY_GROUP[term] + except KeyError as exc: + raise ValueError(f"Clinical-classification term is not approved: {term!r}") from exc diff --git a/src/datahub/association_evidence_v2/contracts.py b/src/datahub/association_evidence_v2/contracts.py new file mode 100644 index 0000000..25995ed --- /dev/null +++ b/src/datahub/association_evidence_v2/contracts.py @@ -0,0 +1,40 @@ +"""Shared contracts for association evidence model v2.""" + +from __future__ import annotations + +SCHEMA_VERSION = "2.4.0-rc1" + +SOURCE_DETAIL_UNAVAILABLE_FIELDS = ( + "effect_allele", + "effect_size", + "standard_error", + "sample_size", + "ancestry", + "fine_mapping", +) + +PROVIDER_DETAIL_AVAILABLE = "available" +PROVIDER_DETAIL_UNAVAILABLE = "unavailable" +PROVIDER_DETAIL_NOT_APPLICABLE = "not_applicable" + +ASSOCIATION_RECORD_KIND_PROVIDER = "source_association_observation" +ASSOCIATION_RECORD_KIND_SOURCE_SUMMARY = "source_summary_association" +EVIDENCE_GRANULARITY_PROVIDER_RECORD = "provider_record" +EVIDENCE_GRANULARITY_SOURCE_SUMMARY = "source_summary" + +SOURCE_PRIORITY = ( + "legacy_cvd_raw", + "legacy_trait_raw", + "million_veteran_program", +) + +POPULATION_FIELDS = ( + ("gnomad_genome.af.af", "Total"), + ("gnomad_genome.af.af_afr", "African"), + ("gnomad_genome.af.af_amr", "Admixed American"), + ("gnomad_genome.af.af_asj", "Ashkenazi Jewish"), + ("gnomad_genome.af.af_eas", "East Asian"), + ("gnomad_genome.af.af_fin", "Finnish"), + ("gnomad_genome.af.af_nfe", "European"), + ("gnomad_genome.af.af_oth", "Other"), +) diff --git a/src/datahub/association_evidence_v2/normalization.py b/src/datahub/association_evidence_v2/normalization.py new file mode 100644 index 0000000..92abb86 --- /dev/null +++ b/src/datahub/association_evidence_v2/normalization.py @@ -0,0 +1,121 @@ +"""Deterministic normalization helpers for association evidence v2.""" + +from __future__ import annotations + +import ast +import hashlib +import json +import math +import re +from typing import Any + + +MISSING_TOKENS = {"", "na", "nan", "none", "null", "n/a", "."} +SEQUENCE_ALLELE = re.compile(r"^[ACGTN]+$", flags=re.IGNORECASE) + + +def canonical_json(value: Any) -> str: + return json.dumps(value, sort_keys=True, separators=(",", ":"), ensure_ascii=True) + + +def stable_id(prefix: str, value: Any, *, length: int = 24) -> str: + digest = hashlib.sha256(canonical_json(value).encode("utf-8")).hexdigest() + return f"{prefix}:{digest[:length]}" + + +def clean(value: Any) -> str | None: + if value is None: + return None + text = str(value).strip() + return None if text.lower() in MISSING_TOKENS else text + + +def finite_number(value: Any) -> float | None: + text = clean(value) + if text is None: + return None + try: + parsed = float(text) + except (TypeError, ValueError): + return None + return parsed if math.isfinite(parsed) else None + + +def phenotype_slug(value: Any) -> str: + text = clean(value) or "" + return re.sub(r"\s+", "_", text.replace("/", "_").strip().lower()) + + +def normalized_alleles(value: Any) -> list[str]: + text = clean(value) + if text is None: + return [] + try: + parsed = ast.literal_eval(text) + except (ValueError, SyntaxError): + # Some archived association exports use unquoted list notation such as + # ``[A, G]``. Strip only a balanced outer bracket pair before applying + # the delimiter parser; allele roles remain intentionally unordered. + candidate = text[1:-1] if text.startswith("[") and text.endswith("]") else text + parsed = [item.strip() for item in re.split(r"[,/|]", candidate)] + if not isinstance(parsed, (list, tuple, set)): + parsed = [parsed] + result = [] + for item in parsed: + allele = clean(item) + if allele is not None: + result.append(allele.upper()) + return sorted(set(result)) + + +def variation_type_from_alleles(alleles: list[str]) -> dict[str, Any]: + """Derive sequence class without assigning REF, ALT, or effect-allele roles.""" + + if len(alleles) < 2: + return _unresolved_variation("fewer than two normalized alleles are available") + if any(not SEQUENCE_ALLELE.fullmatch(allele) for allele in alleles): + return _unresolved_variation("symbolic or non-sequence allele is present") + lengths = {len(allele) for allele in alleles} + if lengths == {1}: + value = "SNV" + elif len(lengths) > 1: + value = "INDEL" + else: + value = "MNV" + return { + "value": value, + "status": "derived", + "method": "normalized_unordered_source_allele_lengths", + "reason": "REF, ALT, and effect-allele roles remain unresolved", + } + + +def _unresolved_variation(reason: str) -> dict[str, Any]: + return { + "value": None, + "status": "unresolved", + "method": "normalized_unordered_source_allele_lengths", + "reason": reason, + } + + +def assertion_terms(value: Any) -> list[str]: + text = clean(value) + if text is None: + return [] + try: + parsed = ast.literal_eval(text) + except (ValueError, SyntaxError): + parsed = text + values = parsed if isinstance(parsed, (list, tuple, set)) else [parsed] + return sorted({term for item in values if (term := clean(item)) is not None}) + + +def clinical_conflict(values: list[str]) -> bool: + lowered = {item.lower() for item in values} + if any("conflicting" in item for item in lowered): + return True + pathogenic = any("pathogenic" in item and "benign" not in item for item in lowered) + benign = any("benign" in item for item in lowered) + uncertain = any("uncertain" in item for item in lowered) + return pathogenic and (benign or uncertain) diff --git a/src/datahub/association_evidence_v2/schema.py b/src/datahub/association_evidence_v2/schema.py new file mode 100644 index 0000000..91d44c2 --- /dev/null +++ b/src/datahub/association_evidence_v2/schema.py @@ -0,0 +1,340 @@ +"""DuckDB schema for the production-scale association evidence v2 sidecar.""" + +from __future__ import annotations + +from typing import Any + +from .contracts import SCHEMA_VERSION + + +DDL = r""" +CREATE TABLE IF NOT EXISTS build_metadata ( + schema_version VARCHAR NOT NULL, + release_id VARCHAR NOT NULL, + built_at TIMESTAMP, + completed_at TIMESTAMP, + status VARCHAR NOT NULL, + input_root VARCHAR NOT NULL, + variant_index_root VARCHAR, + builder_git_commit VARCHAR, + builder_command VARCHAR, + runtime_json JSON, + counts_json JSON, + limitations_json JSON +); + +CREATE TABLE IF NOT EXISTS phenotype_registry ( + dataset_type VARCHAR NOT NULL, + phenotype_slug VARCHAR NOT NULL, + phenotype_path_json JSON NOT NULL, + phenotype_path_key VARCHAR NOT NULL, + PRIMARY KEY (dataset_type, phenotype_slug) +); + +CREATE TABLE IF NOT EXISTS source_completeness ( + source VARCHAR PRIMARY KEY, + dataset_types_json JSON NOT NULL, + provider_detail_status VARCHAR NOT NULL, + provider_detail_reason VARCHAR, + missing_fields_json JSON NOT NULL, + retained_source_summary_artifact VARCHAR, + input_paths_json JSON NOT NULL, + source_version VARCHAR, + source_release VARCHAR, + notes VARCHAR +); + +CREATE TABLE IF NOT EXISTS source_summary_artifacts ( + source_summary_artifact_id VARCHAR PRIMARY KEY, + source VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + logical_artifact VARCHAR NOT NULL, + file_format VARCHAR NOT NULL, + size_bytes UBIGINT NOT NULL, + modified_at_ns UBIGINT NOT NULL, + provider_detail_status VARCHAR NOT NULL, + provider_detail_reason VARCHAR, + missing_fields_json JSON NOT NULL, + content_contract VARCHAR NOT NULL, + publication_mode VARCHAR NOT NULL, + UNIQUE (source, dataset_type, gene_id, logical_artifact) +); + +CREATE TABLE IF NOT EXISTS provider_records ( + provider_record_id VARCHAR NOT NULL, + association_record_id VARCHAR NOT NULL, + dataset_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + source VARCHAR NOT NULL, + source_file_logical VARCHAR NOT NULL, + source_file_absolute VARCHAR NOT NULL, + source_line UBIGINT NOT NULL, + source_row_sha256 VARCHAR NOT NULL, + marker_id_raw VARCHAR, + p_value_raw VARCHAR, + gwas_summary_path_raw VARCHAR, + study_id_raw VARCHAR, + phenotype_raw VARCHAR, + phenotype_slug VARCHAR, + study_raw VARCHAR, + pmid_raw VARCHAR, + study_genome_build_raw VARCHAR, + variant_id_raw VARCHAR, + dbsnp_build_raw VARCHAR, + allele_string_raw VARCHAR, + chromosome_raw VARCHAR, + hg19_start_raw VARCHAR, + hg19_end_raw VARCHAR, + source_variation_type_raw VARCHAR, + frequency_total_raw VARCHAR, + frequency_afr_raw VARCHAR, + frequency_amr_raw VARCHAR, + frequency_asj_raw VARCHAR, + frequency_eas_raw VARCHAR, + frequency_fin_raw VARCHAR, + frequency_nfe_raw VARCHAR, + frequency_oth_raw VARCHAR, + gene_id_raw VARCHAR, + consequence_raw VARCHAR, + putative_impact_raw VARCHAR, + feature_id_raw VARCHAR, + hgvs_p_raw VARCHAR, + protein_length_raw VARCHAR, + dbnsfp_chrom_raw VARCHAR, + dbnsfp_hg18_start_raw VARCHAR, + dbnsfp_hg18_end_raw VARCHAR, + dbnsfp_hg19_start_raw VARCHAR, + dbnsfp_hg19_end_raw VARCHAR, + dbnsfp_hg38_start_raw VARCHAR, + dbnsfp_hg38_end_raw VARCHAR, + ensembl_protein_id_raw VARCHAR, + ensembl_transcript_id_raw VARCHAR, + clinical_significance_raw VARCHAR, + ingested_at TIMESTAMP NOT NULL +); + +CREATE TABLE IF NOT EXISTS variants ( + variant_id VARCHAR PRIMARY KEY, + normalized_allele_contexts_json JSON NOT NULL, + coordinate_contexts_json JSON NOT NULL, + source_reported_variation_types_json JSON NOT NULL, + variation_type VARCHAR, + variation_type_status VARCHAR NOT NULL, + variation_type_method VARCHAR NOT NULL, + variation_type_reason VARCHAR, + ref_allele_status VARCHAR NOT NULL, + alt_allele_status VARCHAR NOT NULL, + effect_allele_status VARCHAR NOT NULL, + multiallelic_observed BOOLEAN NOT NULL, + build_or_position_conflict BOOLEAN NOT NULL +); + +CREATE TABLE IF NOT EXISTS association_records ( + association_record_id VARCHAR NOT NULL, + record_kind VARCHAR NOT NULL, + evidence_granularity VARCHAR NOT NULL, + source_display_name VARCHAR, + dataset_id VARCHAR, + dataset_type VARCHAR NOT NULL, + source VARCHAR NOT NULL, + gene_id_scope VARCHAR, + variant_id VARCHAR NOT NULL, + phenotype_raw VARCHAR, + phenotype_slug VARCHAR NOT NULL, + phenotype_path_json JSON NOT NULL, + phenotype_path_key VARCHAR NOT NULL, + study_id VARCHAR, + study_title VARCHAR, + pmid VARCHAR, + gwas_summary_path VARCHAR, + reported_genome_build VARCHAR, + reported_p_value DOUBLE, + variation_type VARCHAR, + effect_allele VARCHAR, + effect_allele_status VARCHAR NOT NULL, + effect_size DOUBLE, + effect_size_type VARCHAR, + effect_size_status VARCHAR NOT NULL, + standard_error DOUBLE, + standard_error_status VARCHAR NOT NULL, + sample_size DOUBLE, + sample_size_status VARCHAR NOT NULL, + ancestry_json JSON, + ancestry_status VARCHAR NOT NULL, + fine_mapping_json JSON, + fine_mapping_status VARCHAR NOT NULL, + provider_detail_status VARCHAR NOT NULL, + provider_detail_reason VARCHAR, + missing_fields_json JSON NOT NULL, + retained_source_summary_artifact VARCHAR, + retained_source_summary_json JSON, + provider_record_count UBIGINT NOT NULL, + source_file_count UBIGINT NOT NULL +); + +CREATE TABLE IF NOT EXISTS association_record_genes ( + association_record_id VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + gene_linkage_method VARCHAR NOT NULL, + provider_record_count UBIGINT NOT NULL +); + +CREATE TABLE IF NOT EXISTS association_record_provider_records ( + association_record_id VARCHAR NOT NULL, + provider_record_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS consequence_annotations ( + consequence_annotation_id VARCHAR PRIMARY KEY, + variant_id VARCHAR NOT NULL, + gene_id VARCHAR, + consequence VARCHAR NOT NULL, + putative_impact VARCHAR, + transcript_or_feature_id VARCHAR, + ensembl_transcript_id VARCHAR, + ensembl_protein_id VARCHAR, + hgvs_protein VARCHAR, + protein_length DOUBLE, + annotation_source VARCHAR NOT NULL, + annotation_version VARCHAR, + reported_genome_build VARCHAR, + provider_detail_status VARCHAR NOT NULL, + severity_selection VARCHAR NOT NULL, + provider_record_count UBIGINT NOT NULL +); + +CREATE TABLE IF NOT EXISTS consequence_annotation_provider_records ( + consequence_annotation_id VARCHAR NOT NULL, + provider_record_id VARCHAR NOT NULL, + association_record_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS clinical_assertions ( + clinical_assertion_id VARCHAR PRIMARY KEY, + variant_id VARCHAR NOT NULL, + clinical_significance VARCHAR NOT NULL, + normalized_term VARCHAR NOT NULL, + display_group VARCHAR NOT NULL, + assertion_source VARCHAR NOT NULL, + assertion_version VARCHAR, + condition_name VARCHAR, + condition_status VARCHAR NOT NULL, + provider_detail_status VARCHAR NOT NULL, + evidence_granularity VARCHAR NOT NULL, + provenance_limitation VARCHAR NOT NULL, + raw_source_values_json JSON NOT NULL, + provider_record_count UBIGINT NOT NULL +); + +CREATE TABLE IF NOT EXISTS clinical_assertion_provider_records ( + clinical_assertion_id VARCHAR NOT NULL, + provider_record_id VARCHAR NOT NULL, + association_record_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS population_observations ( + population_observation_id VARCHAR NOT NULL, + variant_id VARCHAR NOT NULL, + frequency_source VARCHAR NOT NULL, + frequency_source_version VARCHAR, + population VARCHAR NOT NULL, + reported_frequency DOUBLE NOT NULL, + allele_context_json JSON NOT NULL, + reported_genome_build VARCHAR, + linkage_method VARCHAR NOT NULL, + association_allele_status VARCHAR NOT NULL, + frequency_allele_status VARCHAR NOT NULL, + provider_detail_status VARCHAR NOT NULL, + provider_record_count UBIGINT NOT NULL +); + +CREATE TABLE IF NOT EXISTS population_observation_provider_records ( + population_observation_id VARCHAR NOT NULL, + provider_record_id VARCHAR NOT NULL, + association_record_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS variant_phenotype_summaries ( + variant_phenotype_summary_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + variant_id VARCHAR NOT NULL, + phenotype_slug VARCHAR NOT NULL, + phenotype_path_json JSON NOT NULL, + phenotype_path_key VARCHAR NOT NULL, + source_observation_count UBIGINT NOT NULL, + association_record_count UBIGINT NOT NULL, + provider_association_record_count UBIGINT NOT NULL, + retained_source_summary_count UBIGINT NOT NULL, + study_count UBIGINT NOT NULL, + population_observation_count UBIGINT NOT NULL, + consequence_annotation_count UBIGINT NOT NULL, + clinical_assertion_count UBIGINT NOT NULL, + minimum_reported_p_value DOUBLE, + minimum_reported_p_value_association_record_id VARCHAR, + minimum_reported_p_value_tied_record_ids_json JSON NOT NULL, + maximum_sample_size DOUBLE, + maximum_sample_size_status VARCHAR NOT NULL, + effect_direction_summary_json JSON, + effect_direction_status VARCHAR NOT NULL, + fine_mapping_summary_json JSON, + fine_mapping_status VARCHAR NOT NULL, + consequence_conflict BOOLEAN NOT NULL, + clinical_significance_conflict BOOLEAN NOT NULL, + variation_type VARCHAR, + variation_type_status VARCHAR NOT NULL, + source_completeness_json JSON NOT NULL +); + +CREATE TABLE IF NOT EXISTS summary_consequence_annotations ( + variant_phenotype_summary_id VARCHAR NOT NULL, + consequence_annotation_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS summary_clinical_assertions ( + variant_phenotype_summary_id VARCHAR NOT NULL, + clinical_assertion_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS summary_association_records ( + variant_phenotype_summary_id VARCHAR NOT NULL, + association_record_id VARCHAR NOT NULL +); + +CREATE TABLE IF NOT EXISTS qa_results ( + check_name VARCHAR PRIMARY KEY, + status VARCHAR NOT NULL, + observed_value VARCHAR, + expected_value VARCHAR, + details_json JSON NOT NULL, + checked_at TIMESTAMP NOT NULL +); +""" + + +INDEX_DDL = r""" +CREATE INDEX IF NOT EXISTS idx_source_summary_gene ON source_summary_artifacts(gene_id, dataset_type); +CREATE INDEX IF NOT EXISTS idx_association_id ON association_records(association_record_id); +CREATE INDEX IF NOT EXISTS idx_summary_gene ON variant_phenotype_summaries(gene_id); +""" + + +def create_schema(connection: Any) -> None: + """Create all v2 tables without replacing existing sidecar contents.""" + + connection.execute(DDL) + + +def create_indexes(connection: Any) -> None: + """Create bounded serving indexes after bulk publication. + + Billion-row provenance links are served from hash-partitioned Parquet sidecars; + building global ART indexes for those tables is intentionally unsupported. + """ + + connection.execute(INDEX_DDL) + + +def schema_version() -> str: + return SCHEMA_VERSION diff --git a/src/datahub/association_evidence_v2/serving.py b/src/datahub/association_evidence_v2/serving.py new file mode 100644 index 0000000..30178cf --- /dev/null +++ b/src/datahub/association_evidence_v2/serving.py @@ -0,0 +1,716 @@ +"""Bounded serving sidecars for the normalized association evidence v2 release.""" + +from __future__ import annotations + +import hashlib +import json +import logging +import os +import platform +import shutil +import threading +import time +from contextlib import contextmanager +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +LOGGER = logging.getLogger(__name__) + +REUSED_TABLES = { + "provider_relations": "provider_relation_partitions", + "provider_variant_projection": "provider_variant_partitions", + "summary_base": "summary_base_partitions", + "consequence_annotations": "published_variant_partitions/consequence_annotations", + "clinical_assertions": "published_variant_partitions/clinical_assertions", + "population_observations": "published_variant_partitions/population_observations", +} + +FINE_VARIANT_TABLES = { + "consequence_annotations": ("consequence_annotations", "variant_id"), + "clinical_assertions": ("clinical_assertions", "variant_id"), + "population_observations": ("population_observations", "variant_id"), +} + +COARSE_PASSTHROUGH_TABLES = ( + "source_summary_associations_by_gene", + "source_summary_association_phenotype_counts_by_gene", + "source_summary_association_base_by_gene", +) + + +def _utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _hardlink_tree(source: Path, destination: Path) -> int: + """Create a self-contained tree without duplicating file data on one volume.""" + + files = 0 + for source_path in sorted(source.rglob("*")): + relative = source_path.relative_to(source) + target = destination / relative + if source_path.is_dir(): + target.mkdir(parents=True, exist_ok=True) + continue + target.parent.mkdir(parents=True, exist_ok=True) + if target.exists(): + source_stat = source_path.stat() + target_stat = target.stat() + if (source_stat.st_ino, source_stat.st_size) != ( + target_stat.st_ino, + target_stat.st_size, + ): + raise RuntimeError(f"Serving hardlink target conflicts: {target}") + else: + os.link(source_path, target) + files += 1 + return files + + +@contextmanager +def _progress_monitor(label: str, root: Path, interval_seconds: int): + """Report bounded filesystem progress while DuckDB executes one atomic COPY.""" + + stopped = threading.Event() + started = time.perf_counter() + + def report() -> None: + while not stopped.wait(interval_seconds): + files = 0 + bytes_written = 0 + if root.exists(): + for path in root.rglob("*.parquet"): + try: + stat = path.stat() + except FileNotFoundError: + continue + files += 1 + bytes_written += stat.st_size + LOGGER.info( + "%s progress elapsed_seconds=%.1f parquet_files=%d bytes=%d root=%s", + label, + time.perf_counter() - started, + files, + bytes_written, + root, + ) + + thread = threading.Thread(target=report, name=f"{label}-progress", daemon=True) + thread.start() + try: + yield + finally: + stopped.set() + thread.join(timeout=max(1, interval_seconds)) + + +class AssociationEvidenceV2ServingBuilder: + """Publish immutable hash partitions used by the backend query service.""" + + def __init__( + self, + *, + source_db: Path, + intermediate_root: Path, + output_root: Path, + memory_limit: str = "8GB", + threads: int = 1, + variant_bucket_characters: int = 3, + progress_interval: int = 60, + coarse_serving_root: Path | None = None, + ) -> None: + if variant_bucket_characters not in {2, 3, 4}: + raise ValueError("variant_bucket_characters must be 2, 3, or 4") + if progress_interval < 1: + raise ValueError("progress_interval must be at least 1 second") + self.source_db = source_db.resolve() + self.intermediate_root = intermediate_root.resolve() + self.output_root = output_root.resolve() + self.memory_limit = memory_limit + self.threads = threads + self.variant_bucket_characters = variant_bucket_characters + self.progress_interval = progress_interval + self.coarse_serving_root = ( + coarse_serving_root.resolve() if coarse_serving_root else None + ) + self.tables_root = self.output_root / "tables" + self.checkpoint_path = self.output_root / "serving-checkpoint.json" + self.manifest_path = self.output_root / "serving-manifest.json" + self.log_path = self.output_root / "serving-build.log" + + def _state(self) -> dict[str, Any]: + if not self.checkpoint_path.exists(): + return {"completed": {}} + return json.loads(self.checkpoint_path.read_text()) + + def _save_state(self, state: dict[str, Any]) -> None: + state["updated_at"] = _utc_now() + temporary = self.checkpoint_path.with_suffix(".json.tmp") + temporary.write_text(json.dumps(state, indent=2, sort_keys=True) + "\n") + temporary.replace(self.checkpoint_path) + + def _partition_source(self, logical_name: str) -> tuple[Path, int, str]: + if self.coarse_serving_root is None: + raise RuntimeError( + "A lower-width --coarse-serving-root is required when publishing " + "more than 256 runtime buckets." + ) + manifest_path = self.coarse_serving_root / "serving-manifest.json" + if not manifest_path.is_file(): + raise FileNotFoundError( + f"Coarse serving manifest is missing: {manifest_path}" + ) + manifest = json.loads(manifest_path.read_text()) + if manifest.get("contract") != "association_evidence_v2_partitioned_serving": + raise RuntimeError(f"Invalid coarse serving contract: {manifest_path}") + table = (manifest.get("tables") or {}).get(logical_name) + if not isinstance(table, dict) or not table.get("path"): + raise RuntimeError( + f"Coarse serving manifest does not register {logical_name}" + ) + characters = int(table.get("bucket_characters") or 2) + layout = str(table.get("directory_layout") or "flat_variant_bucket") + expected_layout = ( + "flat_variant_bucket" if characters == 2 else "nested_coarse_fine" + ) + if ( + characters < 2 + or characters >= self.variant_bucket_characters + or layout != expected_layout + ): + raise RuntimeError( + f"Expected a compatible lower-width partition table for {logical_name}; " + f"target_characters={self.variant_bucket_characters} " + f"found characters={characters} layout={layout}" + ) + root = (self.coarse_serving_root / str(table["path"])).resolve() + if self.coarse_serving_root not in root.parents: + raise RuntimeError(f"Coarse serving table escapes its root: {root}") + if not root.is_dir(): + raise FileNotFoundError(f"Coarse serving table is missing: {root}") + return root, characters, layout + + def _reuse_coarse_generated_table( + self, + logical_name: str, + state: dict[str, Any], + ) -> None: + """Hardlink an unchanged gene projection from the verified coarse package.""" + + if self.coarse_serving_root is None: + return + destination = self.tables_root / logical_name + if destination.exists(): + return + manifest_path = self.coarse_serving_root / "serving-manifest.json" + manifest = json.loads(manifest_path.read_text()) + table = (manifest.get("tables") or {}).get(logical_name) + if not isinstance(table, dict) or not table.get("path"): + raise RuntimeError( + f"Coarse serving manifest does not register {logical_name}" + ) + source = (self.coarse_serving_root / str(table["path"])).resolve() + if not source.is_dir(): + raise FileNotFoundError(f"Reusable coarse serving table is missing: {source}") + started = time.perf_counter() + files = _hardlink_tree(source, destination) + state.setdefault("completed", {})[logical_name] = { + "rows": table.get("rows"), + "files": files, + "complete": True, + "reused_from": str(source), + "elapsed_seconds": round(time.perf_counter() - started, 3), + } + self._save_state(state) + + def _publish_variant_table( + self, + *, + connection: Any, + state: dict[str, Any], + logical_name: str, + source_table: str, + key_column: str, + where_sql: str = "", + ) -> tuple[Path, int]: + """Publish one variant table with table-atomic or coarse-unit checkpoints.""" + + completed = state.setdefault("completed", {}) + target = self.tables_root / logical_name + count = int( + connection.execute( + f"SELECT count(*) FROM {source_table} {where_sql}" + ).fetchone()[0] + ) + previous = completed.get(logical_name) or {} + partition_source = None + if self.variant_bucket_characters > 2: + partition_source = self._partition_source(logical_name) + source_token = ( + "|".join(str(value) for value in partition_source) + if partition_source + else str(self.source_db) + ) + configuration_matches = ( + previous.get("rows") == count + and previous.get("variant_bucket_characters") + == self.variant_bucket_characters + and previous.get("source_token") == source_token + ) + if not configuration_matches: + if target.exists(): + stale = target.with_name( + target.name + + ".stale-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + target.rename(stale) + previous = { + "rows": count, + "variant_bucket_characters": self.variant_bucket_characters, + "source_token": source_token, + "completed_coarse_buckets": [], + "complete": False, + } + completed[logical_name] = previous + self._save_state(state) + if ( + previous.get("complete") + and target.is_dir() + and any(target.rglob("*.parquet")) + ): + return target, count + + if self.variant_bucket_characters == 2: + target.parent.mkdir(parents=True, exist_ok=True) + LOGGER.info("%s partitioning start rows=%d", logical_name, count) + where_clause = f" {where_sql}" if where_sql else "" + with _progress_monitor(logical_name, target, self.progress_interval): + connection.execute( + f""" +COPY ( + SELECT *, substr(sha256({key_column}), 1, 2) AS variant_bucket + FROM {source_table}{where_clause} +) TO '{str(target).replace("'", "''")}' +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + previous.update( + { + "files": len(list(target.rglob("*.parquet"))), + "complete": True, + } + ) + self._save_state(state) + return target, count + + target.mkdir(parents=True, exist_ok=True) + work_root = self.output_root / ".serving-work" / logical_name + work_root.mkdir(parents=True, exist_ok=True) + done = set(previous.get("completed_coarse_buckets") or []) + source_root, source_characters, source_layout = partition_source + source_directory_prefix = ( + "variant_bucket=" + if source_layout == "flat_variant_bucket" + else "coarse_bucket=" + ) + coarse_directories = sorted( + path + for path in source_root.iterdir() + if path.is_dir() + and path.name.startswith(source_directory_prefix) + and len(path.name.rsplit("=", 1)[-1]) == 2 + ) + if not coarse_directories: + raise RuntimeError( + f"No two-character coarse units found in partition source: {source_root}" + ) + for index, source_directory in enumerate(coarse_directories, start=1): + coarse_bucket = source_directory.name.rsplit("=", 1)[-1] + final_unit = target / f"coarse_bucket={coarse_bucket}" + if coarse_bucket in done and final_unit.is_dir(): + continue + if final_unit.exists(): + raise RuntimeError( + f"Uncheckpointed fine partition unit already exists: {final_unit}" + ) + temporary_unit = work_root / f"coarse_bucket={coarse_bucket}.incomplete" + if temporary_unit.exists(): + shutil.rmtree(temporary_unit) + if source_characters == 2: + source_glob_path = source_directory / "*.parquet" + excluded_partition_columns = "variant_bucket" + else: + source_glob_path = source_directory / "variant_bucket=*" / "*.parquet" + excluded_partition_columns = "variant_bucket, coarse_bucket" + source_glob = str(source_glob_path).replace("'", "''") + LOGGER.info( + "%s fine partition start coarse_bucket=%s unit=%d/%d", + logical_name, + coarse_bucket, + index, + len(coarse_directories), + ) + with _progress_monitor( + f"{logical_name} coarse_bucket={coarse_bucket}", + temporary_unit, + self.progress_interval, + ): + connection.execute( + f""" +COPY ( + SELECT * EXCLUDE ({excluded_partition_columns}), + substr(sha256({key_column}), 1, {self.variant_bucket_characters}) AS variant_bucket + FROM read_parquet('{source_glob}', hive_partitioning=true) +) TO '{str(temporary_unit).replace("'", "''")}' +(FORMAT PARQUET, PARTITION_BY (variant_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + temporary_unit.rename(final_unit) + done.add(coarse_bucket) + previous["completed_coarse_buckets"] = sorted(done) + previous["files"] = len(list(target.rglob("*.parquet"))) + self._save_state(state) + LOGGER.info( + "%s fine partition complete coarse_buckets=%d/%d files=%d", + logical_name, + len(done), + len(coarse_directories), + previous["files"], + ) + previous["complete"] = True + previous["files"] = len(list(target.rglob("*.parquet"))) + self._save_state(state) + return target, count + + def build(self) -> dict[str, Any]: + import duckdb + + self.output_root.mkdir(parents=True, exist_ok=True) + self.tables_root.mkdir(parents=True, exist_ok=True) + state = self._state() + completed = state.setdefault("completed", {}) + reused: dict[str, Any] = {} + for logical_name, relative in REUSED_TABLES.items(): + if ( + self.variant_bucket_characters > 2 + and logical_name in FINE_VARIANT_TABLES + ): + continue + source = self.intermediate_root / relative + if not source.is_dir(): + raise FileNotFoundError(f"Required build partition is missing: {source}") + destination = self.tables_root / logical_name + started = time.perf_counter() + file_count = _hardlink_tree(source, destination) + reused[logical_name] = { + "path": str(destination.relative_to(self.output_root)), + "partition_key": ( + "association_record_id" if logical_name == "provider_relations" else "variant_id" + ), + "bucket_function": ( + "association_record_id[12:14]" + if logical_name == "provider_relations" + else "sha256(variant_id)[0:2]" + ), + "files": file_count, + "hardlinked_from": str(source), + } + completed[logical_name] = { + "files": file_count, + "elapsed_seconds": round(time.perf_counter() - started, 3), + } + self._save_state(state) + + connection = duckdb.connect( + str(self.source_db), + read_only=True, + config={"threads": str(self.threads), "memory_limit": self.memory_limit}, + ) + connection.execute("SET preserve_insertion_order=false") + connection.execute("PRAGMA disable_progress_bar") + provider_started = time.perf_counter() + provider_root, provider_count = self._publish_variant_table( + connection=connection, + state=state, + logical_name="provider_records", + source_table="provider_records", + key_column="variant_id_raw", + where_sql="WHERE coalesce(trim(variant_id_raw), '')<>''", + ) + provider_total_count = int( + connection.execute("SELECT count(*) FROM provider_records").fetchone()[0] + ) + provider_unkeyed_count = provider_total_count - provider_count + completed["provider_records"]["elapsed_seconds"] = round( + time.perf_counter() - provider_started, 3 + ) + self._save_state(state) + association_started = time.perf_counter() + association_root, association_count = self._publish_variant_table( + connection=connection, + state=state, + logical_name="association_records", + source_table="association_records", + key_column="variant_id", + ) + completed["association_records"]["elapsed_seconds"] = round( + time.perf_counter() - association_started, 3 + ) + self._save_state(state) + summary_started = time.perf_counter() + summary_root, summary_count = self._publish_variant_table( + connection=connection, + state=state, + logical_name="variant_phenotype_summaries", + source_table="variant_phenotype_summaries", + key_column="variant_id", + ) + completed["variant_phenotype_summaries"]["elapsed_seconds"] = round( + time.perf_counter() - summary_started, 3 + ) + self._save_state(state) + fine_variant_tables: dict[str, Any] = {} + if self.variant_bucket_characters > 2: + for logical_name, (source_table, key_column) in FINE_VARIANT_TABLES.items(): + started = time.perf_counter() + table_root, table_count = self._publish_variant_table( + connection=connection, + state=state, + logical_name=logical_name, + source_table=source_table, + key_column=key_column, + ) + completed[logical_name]["elapsed_seconds"] = round( + time.perf_counter() - started, 3 + ) + self._save_state(state) + fine_variant_tables[logical_name] = { + "path": str(table_root.relative_to(self.output_root)), + "partition_key": key_column, + "bucket_function": ( + f"sha256({key_column})[0:{self.variant_bucket_characters}]" + ), + "bucket_characters": self.variant_bucket_characters, + "coarse_bucket_characters": 2, + "directory_layout": "nested_coarse_fine", + **completed[logical_name], + } + passthrough_tables: dict[str, Any] = {} + if self.coarse_serving_root is not None: + coarse_manifest_path = self.coarse_serving_root / "serving-manifest.json" + coarse_manifest = json.loads(coarse_manifest_path.read_text()) + for logical_name in COARSE_PASSTHROUGH_TABLES: + source_metadata = (coarse_manifest.get("tables") or {}).get(logical_name) + if not isinstance(source_metadata, dict): + continue + self._reuse_coarse_generated_table(logical_name, state) + source_table = dict(source_metadata) + source_table["path"] = str( + (self.tables_root / logical_name).relative_to(self.output_root) + ) + source_table["reused_from_manifest"] = str(coarse_manifest_path) + source_table.update(completed[logical_name]) + passthrough_tables[logical_name] = source_table + self._reuse_coarse_generated_table("summary_base_by_gene", state) + summary_gene_root = self.tables_root / "summary_base_by_gene" + summary_gene_state = completed.get("summary_base_by_gene") or {} + if not ( + summary_gene_state.get("complete") + and summary_gene_root.is_dir() + and any(summary_gene_root.rglob("*.parquet")) + ): + if summary_gene_root.exists(): + stale = summary_gene_root.with_name( + summary_gene_root.name + ".stale-" + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + summary_gene_root.rename(stale) + started = time.perf_counter() + source = str((self.tables_root / "summary_base" / "**" / "*.parquet")).replace("'", "''") + LOGGER.info("Gene summary partitioning start") + copied = connection.execute( + f""" +COPY ( + SELECT * EXCLUDE (variant_bucket), substr(sha256(gene_id), 1, 2) AS gene_bucket + FROM read_parquet('{source}', hive_partitioning=true) +) TO '{str(summary_gene_root).replace("'", "''")}' +(FORMAT PARQUET, PARTITION_BY (gene_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ).fetchone() + completed["summary_base_by_gene"] = { + "rows": int(copied[0]) if copied else None, + "files": len(list(summary_gene_root.rglob("*.parquet"))), + "elapsed_seconds": round(time.perf_counter() - started, 3), + "complete": True, + } + self._save_state(state) + self._reuse_coarse_generated_table("consequence_annotations_by_gene", state) + consequence_gene_root = self.tables_root / "consequence_annotations_by_gene" + consequence_gene_count = int( + connection.execute( + "SELECT count(*) FROM consequence_annotations " + "WHERE coalesce(trim(gene_id), '')<>''" + ).fetchone()[0] + ) + consequence_gene_state = completed.get("consequence_annotations_by_gene") or {} + if not ( + consequence_gene_state.get("complete") + and consequence_gene_state.get("rows") == consequence_gene_count + and consequence_gene_root.is_dir() + and any(consequence_gene_root.rglob("*.parquet")) + ): + if consequence_gene_root.exists(): + stale = consequence_gene_root.with_name( + consequence_gene_root.name + + ".stale-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + consequence_gene_root.rename(stale) + started = time.perf_counter() + LOGGER.info( + "Gene consequence partitioning start rows=%d", consequence_gene_count + ) + connection.execute( + f""" +COPY ( + SELECT *, substr(sha256(upper(gene_id)), 1, 2) AS gene_bucket + FROM consequence_annotations + WHERE coalesce(trim(gene_id), '')<>'' +) TO '{str(consequence_gene_root).replace("'", "''")}' +(FORMAT PARQUET, PARTITION_BY (gene_bucket), COMPRESSION ZSTD, ROW_GROUP_SIZE 100000) +""" + ) + completed["consequence_annotations_by_gene"] = { + "rows": consequence_gene_count, + "files": len(list(consequence_gene_root.rglob("*.parquet"))), + "elapsed_seconds": round(time.perf_counter() - started, 3), + "complete": True, + } + self._save_state(state) + connection.close() + + db_sha_path = self.source_db.with_suffix(self.source_db.suffix + ".sha256") + source_db_sha = None + if db_sha_path.is_file(): + source_db_sha = db_sha_path.read_text().split()[0] + manifest = { + "contract": "association_evidence_v2_partitioned_serving", + "schema_version": "2.6.0-rc1", + "built_at": _utc_now(), + "source_db": str(self.source_db), + "source_db_sha256": source_db_sha, + "coarse_serving_manifest": ( + str(self.coarse_serving_root / "serving-manifest.json") + if self.coarse_serving_root + else None + ), + "coarse_serving_manifest_sha256": ( + _sha256(self.coarse_serving_root / "serving-manifest.json") + if self.coarse_serving_root + else None + ), + "variant_bucket_count": 16 ** self.variant_bucket_characters, + "runtime": { + "python": platform.python_version(), + "duckdb": duckdb.__version__, + "platform": platform.platform(), + "threads": self.threads, + "memory_limit": self.memory_limit, + "progress_interval": self.progress_interval, + }, + "tables": { + **reused, + **passthrough_tables, + **fine_variant_tables, + "provider_records": { + "path": str(provider_root.relative_to(self.output_root)), + "partition_key": "variant_id_raw", + "bucket_function": ( + f"sha256(variant_id_raw)[0:{self.variant_bucket_characters}]" + ), + "bucket_characters": self.variant_bucket_characters, + "normalized_source_rows": provider_total_count, + "excluded_unkeyed_rows": provider_unkeyed_count, + "exclusion_contract": ( + "Rows without variant_id_raw remain in the normalized sidecar " + "but cannot enter variant-keyed serving or drill-down exports." + ), + "coarse_bucket_characters": ( + 2 if self.variant_bucket_characters > 2 else None + ), + "directory_layout": ( + "nested_coarse_fine" + if self.variant_bucket_characters > 2 + else "flat_variant_bucket" + ), + **completed["provider_records"], + }, + "association_records": { + "path": str(association_root.relative_to(self.output_root)), + "partition_key": "variant_id", + "bucket_function": ( + f"sha256(variant_id)[0:{self.variant_bucket_characters}]" + ), + "bucket_characters": self.variant_bucket_characters, + "coarse_bucket_characters": ( + 2 if self.variant_bucket_characters > 2 else None + ), + "directory_layout": ( + "nested_coarse_fine" + if self.variant_bucket_characters > 2 + else "flat_variant_bucket" + ), + **completed["association_records"], + }, + "variant_phenotype_summaries": { + "path": str(summary_root.relative_to(self.output_root)), + "partition_key": "variant_id", + "bucket_function": ( + f"sha256(variant_id)[0:{self.variant_bucket_characters}]" + ), + "bucket_characters": self.variant_bucket_characters, + "coarse_bucket_characters": ( + 2 if self.variant_bucket_characters > 2 else None + ), + "directory_layout": ( + "nested_coarse_fine" + if self.variant_bucket_characters > 2 + else "flat_variant_bucket" + ), + **completed["variant_phenotype_summaries"], + }, + "summary_base_by_gene": { + "path": str(summary_gene_root.relative_to(self.output_root)), + "partition_key": "gene_id", + "bucket_function": "sha256(gene_id)[0:2]", + "bucket_characters": 2, + **completed["summary_base_by_gene"], + }, + "consequence_annotations_by_gene": { + "path": str(consequence_gene_root.relative_to(self.output_root)), + "partition_key": "gene_id", + "bucket_function": "sha256(upper(gene_id))[0:2]", + "bucket_characters": 2, + **completed["consequence_annotations_by_gene"], + }, + }, + "limitations": [ + "MVP source summaries are first-class associations; provider rows are not fabricated.", + "Provider provenance partitions contain only recoverable AWS source rows.", + "Global ART indexes are intentionally not built for hundred-million/billion-row provenance links.", + ], + } + self.manifest_path.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n") + checksum = _sha256(self.manifest_path) + checksum_path = self.manifest_path.with_suffix(".json.sha256") + checksum_path.write_text(f"{checksum} {self.manifest_path.name}\n") + return manifest diff --git a/src/datahub/association_evidence_v2/source_summary_index.py b/src/datahub/association_evidence_v2/source_summary_index.py new file mode 100644 index 0000000..632426f --- /dev/null +++ b/src/datahub/association_evidence_v2/source_summary_index.py @@ -0,0 +1,503 @@ +"""Build a resumable serving index for first-class MVP summary associations. + +The retained compact variant-index artifacts are valid source-summary association +records. They are indexed without fabricating provider, study, allele, effect, +sample-size, ancestry, or fine-mapping detail that the compact artifact omits. +""" + +from __future__ import annotations + +import concurrent.futures +import csv +import gzip +import hashlib +import json +import logging +import os +import time +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Iterator + +import duckdb + +LOGGER = logging.getLogger(__name__) +TABLE_NAME = "source_summary_associations_by_gene" +CONTRACT = "source_summary_association_index_v2" +MVP_DATASET_ID = "hbp_mvp_association" +PRIVATE_SOURCE_PATH_FIELDS = {"source_file", "source_files"} + + +def _without_private_source_paths(value: Any) -> Any: + """Return the public compact-summary fields without host-local paths.""" + + if isinstance(value, dict): + return { + str(key): cleaned + for key, item in value.items() + if str(key) not in PRIVATE_SOURCE_PATH_FIELDS + and (cleaned := _without_private_source_paths(item)) is not None + } + if isinstance(value, list): + return [ + cleaned + for item in value + if (cleaned := _without_private_source_paths(item)) is not None + ] + if isinstance(value, str) and value.startswith("/"): + return None + return value + + +def _utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _iter_json_array(path: Path, *, chunk_size: int = 1024 * 1024) -> Iterator[dict[str, Any]]: + opener = gzip.open if path.name.endswith(".gz") else open + decoder = json.JSONDecoder() + with opener(path, "rt", encoding="utf-8") as stream: + buffer = "" + position = 0 + started = False + eof = False + while True: + if not eof: + chunk = stream.read(chunk_size) + if chunk: + buffer += chunk + else: + eof = True + while True: + while position < len(buffer) and buffer[position].isspace(): + position += 1 + if not started: + if position >= len(buffer): + break + if buffer[position] != "[": + raise ValueError(f"Expected JSON array in {path}") + started = True + position += 1 + continue + while position < len(buffer) and ( + buffer[position].isspace() or buffer[position] == "," + ): + position += 1 + if position < len(buffer) and buffer[position] == "]": + return + if position >= len(buffer): + break + try: + value, end = decoder.raw_decode(buffer, position) + except json.JSONDecodeError: + if eof: + raise + break + if not isinstance(value, dict): + raise ValueError(f"Expected object entries in {path}") + yield value + position = end + if position: + buffer = buffer[position:] + position = 0 + if eof: + if buffer.strip(): + raise ValueError(f"Truncated JSON array in {path}") + return + + +def _source_summary_id( + source: str, + dataset_type: str, + gene: str, + variant_id: str, + phenotype_path: list[str], + logical_artifact: str, +) -> str: + identity = json.dumps( + [source, dataset_type, gene, variant_id, phenotype_path, logical_artifact], + separators=(",", ":"), + ensure_ascii=True, + ) + return "source-summary:" + hashlib.sha256(identity.encode("utf-8")).hexdigest() + + +def _process_artifact(task: tuple[Any, ...]) -> dict[str, Any]: + ( + source, + gene, + dataset_type, + logical_artifact, + size_bytes, + modified_at_ns, + provider_detail_reason, + missing_fields_json, + content_contract, + publication_mode, + path_text, + output_text, + work_text, + ) = task + path = Path(path_text) + output = Path(output_text) + work_root = Path(work_text) + stat = path.stat() + if stat.st_size != int(size_bytes) or stat.st_mtime_ns != int(modified_at_ns): + raise RuntimeError(f"Registered source-summary artifact changed: {path}") + + token = hashlib.sha256(str(logical_artifact).encode("utf-8")).hexdigest() + csv_path = work_root / f"{token}.csv" + parquet_path = work_root / f"{token}.parquet.incomplete" + work_root.mkdir(parents=True, exist_ok=True) + csv_path.unlink(missing_ok=True) + parquet_path.unlink(missing_ok=True) + selected_rows = 0 + seen_summary_ids: set[str] = set() + with csv_path.open("w", newline="", encoding="utf-8") as stream: + writer = csv.writer(stream) + writer.writerow( + ( + "association_record_id", + "source_summary_id", + "record_kind", + "evidence_granularity", + "source", + "source_display_name", + "dataset_id", + "dataset_type", + "gene_id", + "variant_id", + "phenotype_raw", + "phenotype_slug", + "phenotype_kind", + "phenotype_path_json", + "phenotype_path_key", + "reported_p_value", + "variation_type", + "ancestry_json", + "metadata_json", + "representative_source", + "provider_detail_status", + "provider_detail_reason", + "missing_fields_json", + "retained_source_summary_artifact", + "retained_source_summary_json", + "content_contract", + "publication_mode", + "source_row_ordinal", + ) + ) + for ordinal, entry in enumerate(_iter_json_array(path)): + sources = entry.get("sources") + if not isinstance(sources, list): + sources = [entry.get("source")] + if str(source).casefold() not in { + str(value).strip().casefold() for value in sources if value + }: + continue + variant_id = str(entry.get("variant_id") or "").strip() + if not variant_id: + continue + phenotype_path = entry.get("phenotype_path") + if not isinstance(phenotype_path, list) or not phenotype_path: + phenotype = str(entry.get("phenotype") or "").strip() + phenotype_path = [phenotype] if phenotype else [] + phenotype_path = [ + str(value).strip() for value in phenotype_path if str(value).strip() + ] + if not phenotype_path: + continue + representative_source = str(entry.get("source") or "").strip() + p_value = entry.get("p_value") if representative_source.casefold() == str(source).casefold() else None + try: + p_value = float(p_value) if p_value is not None else None + except (TypeError, ValueError): + p_value = None + source_summary_id = _source_summary_id( + str(source), str(dataset_type), str(gene), variant_id, + phenotype_path, str(logical_artifact), + ) + selected_rows += 1 + if source_summary_id in seen_summary_ids: + continue + seen_summary_ids.add(source_summary_id) + row_missing_fields = json.loads(missing_fields_json) + public_metadata = _without_private_source_paths(entry.get("metadata")) + public_entry = _without_private_source_paths( + {key: value for key, value in entry.items() if value is not None} + ) + if entry.get("ancestry") is not None: + row_missing_fields = [ + field for field in row_missing_fields if field != "ancestry" + ] + writer.writerow( + ( + "association-source-summary:" + + source_summary_id.partition(":")[2], + source_summary_id, + "source_summary_association", + "source_summary", + "million_veteran_program", + "MVP", + MVP_DATASET_ID, + dataset_type, + gene, + variant_id, + entry.get("phenotype") or phenotype_path[-1], + phenotype_path[-1], + entry.get("label_key"), + json.dumps(phenotype_path, separators=(",", ":"), ensure_ascii=True), + " > ".join(phenotype_path), + p_value, + entry.get("variation_type"), + json.dumps(entry.get("ancestry"), separators=(",", ":"), ensure_ascii=True) + if entry.get("ancestry") is not None else None, + json.dumps(public_metadata, separators=(",", ":"), ensure_ascii=True) + if public_metadata is not None else None, + representative_source or None, + "not_applicable", + None, + json.dumps(row_missing_fields, separators=(",", ":")), + logical_artifact, + json.dumps( + public_entry, + separators=(",", ":"), + ensure_ascii=True, + ), + content_contract, + publication_mode, + ordinal, + ) + ) + + connection = duckdb.connect(config={"threads": "1", "memory_limit": "1GB"}) + try: + connection.execute("SET preserve_insertion_order=false") + escaped_csv = str(csv_path).replace("'", "''") + escaped_parquet = str(parquet_path).replace("'", "''") + connection.execute( + f""" +COPY ( + SELECT * + FROM read_csv('{escaped_csv}', header=true, all_varchar=true) +) TO '{escaped_parquet}' +(FORMAT PARQUET, COMPRESSION ZSTD, ROW_GROUP_SIZE 50000) +""" + ) + unique_rows = len(seen_summary_ids) + finally: + connection.close() + csv_path.unlink(missing_ok=True) + output.parent.mkdir(parents=True, exist_ok=True) + if output.exists(): + raise RuntimeError(f"Uncheckpointed source-summary target exists: {output}") + parquet_path.replace(output) + return { + "logical_artifact": logical_artifact, + "source": source, + "gene_id": gene, + "dataset_type": dataset_type, + "input_size_bytes": int(size_bytes), + "input_modified_at_ns": int(modified_at_ns), + "selected_source_rows": selected_rows, + "unique_summary_rows": unique_rows, + "output_path": str(output), + "output_size_bytes": output.stat().st_size, + "output_sha256": _sha256(output), + } + + +class SourceSummaryIndexBuilder: + def __init__( + self, + *, + source_db: Path, + artifact_root: Path, + serving_root: Path, + workers: int = 3, + progress_interval: int = 1000, + ) -> None: + self.source_db = source_db.resolve() + self.artifact_root = artifact_root.resolve() + self.serving_root = serving_root.resolve() + self.workers = max(1, int(workers)) + self.progress_interval = max(1, int(progress_interval)) + self.table_root = self.serving_root / "tables" / TABLE_NAME + self.work_root = self.serving_root / ".source-summary-work" + self.checkpoint_db = self.serving_root / "source-summary-index-checkpoint.duckdb" + self.file_manifest = self.serving_root / "source-summary-index-files.jsonl" + self.log_path = self.serving_root / "source-summary-index.log" + + def build(self) -> dict[str, Any]: + self.serving_root.mkdir(parents=True, exist_ok=True) + self.table_root.mkdir(parents=True, exist_ok=True) + checkpoint = duckdb.connect(str(self.checkpoint_db)) + checkpoint.execute( + """CREATE TABLE IF NOT EXISTS processed_artifacts ( + logical_artifact VARCHAR PRIMARY KEY, + source VARCHAR NOT NULL, + gene_id VARCHAR NOT NULL, + dataset_type VARCHAR NOT NULL, + input_size_bytes UBIGINT NOT NULL, + input_modified_at_ns UBIGINT NOT NULL, + selected_source_rows UBIGINT NOT NULL, + unique_summary_rows UBIGINT NOT NULL, + output_path VARCHAR NOT NULL, + output_size_bytes UBIGINT NOT NULL, + output_sha256 VARCHAR NOT NULL, + processed_at TIMESTAMP NOT NULL + )""" + ) + completed = { + row[0] + for row in checkpoint.execute( + "SELECT logical_artifact FROM processed_artifacts" + ).fetchall() + } + source = duckdb.connect(str(self.source_db), read_only=True) + artifacts = source.execute( + """SELECT source, gene_id, dataset_type, logical_artifact, size_bytes, + modified_at_ns, provider_detail_reason, missing_fields_json, + content_contract, publication_mode + FROM source_summary_artifacts + WHERE provider_detail_status IN ('not_applicable', 'unavailable') + AND publication_mode='gene_scoped_on_demand' + ORDER BY dataset_type, gene_id, logical_artifact""" + ).fetchall() + source.close() + tasks: list[tuple[Any, ...]] = [] + for row in artifacts: + source_name, gene, dataset_type, logical = map(str, row[:4]) + if logical in completed: + continue + source_path = (self.artifact_root / logical).resolve() + source_path.relative_to(self.artifact_root) + gene_bucket = hashlib.sha256(gene.encode("utf-8")).hexdigest()[:2] + gene_key = hashlib.sha256(gene.encode("utf-8")).hexdigest()[:16] + file_token = hashlib.sha256(logical.encode("utf-8")).hexdigest()[:24] + output = ( + self.table_root / f"gene_bucket={gene_bucket}" + / f"gene_key={gene_key}" / f"part-{file_token}.parquet" + ) + if output.exists(): + orphaned = output.with_name( + output.name + + ".uncheckpointed-" + + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ") + ) + output.rename(orphaned) + tasks.append( + ( + source_name, gene, dataset_type, logical, int(row[4]), int(row[5]), + row[6], str(row[7]), row[8], row[9], str(source_path), + str(output), str(self.work_root), + ) + ) + started = time.perf_counter() + processed_count = len(completed) + selected_total = int( + checkpoint.execute( + "SELECT coalesce(sum(selected_source_rows),0) FROM processed_artifacts" + ).fetchone()[0] + ) + unique_total = int( + checkpoint.execute( + "SELECT coalesce(sum(unique_summary_rows),0) FROM processed_artifacts" + ).fetchone()[0] + ) + + def save(result: dict[str, Any]) -> None: + nonlocal processed_count, selected_total, unique_total + checkpoint.execute( + """INSERT INTO processed_artifacts VALUES + (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, current_timestamp)""", + [ + result["logical_artifact"], result["source"], result["gene_id"], + result["dataset_type"], result["input_size_bytes"], + result["input_modified_at_ns"], result["selected_source_rows"], + result["unique_summary_rows"], result["output_path"], + result["output_size_bytes"], result["output_sha256"], + ], + ) + processed_count += 1 + selected_total += int(result["selected_source_rows"]) + unique_total += int(result["unique_summary_rows"]) + if processed_count % self.progress_interval == 0: + LOGGER.info( + "source-summary index progress artifacts=%d/%d selected_rows=%d " + "unique_rows=%d elapsed_seconds=%.1f workers=%d", + processed_count, len(artifacts), selected_total, unique_total, + time.perf_counter() - started, self.workers, + ) + + if self.workers == 1: + for task in tasks: + save(_process_artifact(task)) + else: + with concurrent.futures.ProcessPoolExecutor(max_workers=self.workers) as executor: + for result in executor.map(_process_artifact, tasks, chunksize=1): + save(result) + + rows = checkpoint.execute( + """SELECT logical_artifact, source, gene_id, dataset_type, + selected_source_rows, unique_summary_rows, output_path, + output_size_bytes, output_sha256 + FROM processed_artifacts ORDER BY logical_artifact""" + ).fetchall() + with self.file_manifest.open("w", encoding="utf-8") as stream: + for row in rows: + stream.write(json.dumps(dict(zip( + ("logical_artifact", "source", "gene_id", "dataset_type", + "selected_source_rows", "unique_summary_rows", "output_path", + "output_size_bytes", "output_sha256"), row + )), sort_keys=True) + "\n") + checkpoint.close() + manifest_path = self.serving_root / "serving-manifest.json" + manifest = json.loads(manifest_path.read_text()) + manifest.setdefault("tables", {})[TABLE_NAME] = { + "path": str(self.table_root.relative_to(self.serving_root)), + "rows": unique_total, + "selected_source_rows": selected_total, + "files": len(rows), + "source_artifacts": len(artifacts), + "partition_key": "gene_id", + "bucket_function": "sha256(gene_id)[0:2]", + "bucket_characters": 2, + "gene_key_function": "sha256(gene_id)[0:16]", + "directory_layout": "nested_gene_key", + "contract": CONTRACT, + "provider_detail_status": "not_applicable", + "record_kind": "source_summary_association", + "evidence_granularity": "source_summary", + } + manifest["schema_version"] = "2.9.1-rc3" + manifest["source_summary_index_built_at"] = _utc_now() + temporary = manifest_path.with_suffix(".json.tmp") + temporary.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n") + temporary.replace(manifest_path) + checksum = _sha256(manifest_path) + manifest_path.with_suffix(".json.sha256").write_text( + f"{checksum} {manifest_path.name}\n" + ) + return { + "contract": CONTRACT, + "built_at": manifest["source_summary_index_built_at"], + "source_artifacts": len(artifacts), + "selected_source_rows": selected_total, + "unique_summary_rows": unique_total, + "files": len(rows), + "table_root": str(self.table_root), + "file_manifest": str(self.file_manifest), + "checkpoint_db": str(self.checkpoint_db), + "serving_manifest": str(manifest_path), + "serving_manifest_sha256": checksum, + "elapsed_seconds": round(time.perf_counter() - started, 3), + } diff --git a/src/datahub/association_evidence_v2/source_summary_rollup.py b/src/datahub/association_evidence_v2/source_summary_rollup.py new file mode 100644 index 0000000..0211edc --- /dev/null +++ b/src/datahub/association_evidence_v2/source_summary_rollup.py @@ -0,0 +1,270 @@ +"""Build bounded default-summary rollups from the immutable source-summary index.""" + +from __future__ import annotations + +import hashlib +import json +import logging +import platform +import shutil +import time +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import duckdb + +SOURCE_TABLE = "source_summary_associations_by_gene" +BASE_TABLE = "source_summary_association_base_by_gene" +PHENOTYPE_TABLE = "source_summary_association_phenotype_counts_by_gene" +CONTRACT = "source_summary_association_rollup_v2" + + +def _utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _save_json(path: Path, value: Any) -> None: + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text(json.dumps(value, indent=2, sort_keys=True) + "\n") + temporary.replace(path) + + +class SourceSummaryRollupBuilder: + """Create resumable gene-keyed summary projections without altering raw evidence.""" + + def __init__( + self, + *, + serving_root: Path, + threads: int = 2, + memory_limit: str = "4GB", + progress_interval: int = 1, + logger: logging.Logger | None = None, + ) -> None: + self.serving_root = serving_root.resolve() + self.threads = max(1, int(threads)) + self.memory_limit = memory_limit + self.progress_interval = max(1, int(progress_interval)) + self.logger = logger or logging.getLogger(__name__) + self.tables_root = self.serving_root / "tables" + self.source_root = self.tables_root / SOURCE_TABLE + self.base_root = self.tables_root / BASE_TABLE + self.phenotype_root = self.tables_root / PHENOTYPE_TABLE + self.work_root = self.serving_root / ".source-summary-rollup-work" + self.checkpoint_path = self.serving_root / "source-summary-rollup-checkpoint.json" + self.file_manifest = self.serving_root / "source-summary-rollup-files.jsonl" + + def _validate_source(self) -> dict[str, Any]: + manifest_path = self.serving_root / "serving-manifest.json" + checksum_path = self.serving_root / "serving-manifest.json.sha256" + manifest = json.loads(manifest_path.read_text()) + expected = checksum_path.read_text().split()[0] + if _sha256(manifest_path) != expected: + raise RuntimeError("Serving manifest checksum failed before source-summary rollup.") + source = (manifest.get("tables") or {}).get(SOURCE_TABLE) or {} + if manifest.get("schema_version") != "2.9.0-rc2": + raise RuntimeError("Source-summary rollup requires completed schema 2.9.0-rc2.") + if source.get("contract") != "source_summary_association_index_v2": + raise RuntimeError("Unexpected source-summary input contract.") + if int(source.get("files") or 0) != int(source.get("source_artifacts") or -1): + raise RuntimeError("Source-summary files do not reconcile with registered artifacts.") + return manifest + + @staticmethod + def _uncheckpointed(path: Path) -> Path: + return path.with_name(path.name + ".uncheckpointed-" + datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ")) + + def _rollup_bucket(self, bucket: str) -> dict[str, Any]: + source_bucket = self.source_root / f"gene_bucket={bucket}" + final_base = self.base_root / f"gene_bucket={bucket}" + final_phenotype = self.phenotype_root / f"gene_bucket={bucket}" + if final_base.exists() or final_phenotype.exists(): + raise RuntimeError(f"Uncheckpointed rollup output exists for bucket {bucket}.") + work = self.work_root / f"bucket={bucket}.incomplete" + if work.exists(): + work.rename(self._uncheckpointed(work)) + base_work = work / "base" + phenotype_work = work / "phenotype" + base_work.mkdir(parents=True) + phenotype_work.mkdir(parents=True) + source_glob = str(source_bucket / "gene_key=*" / "*.parquet").replace("'", "''") + base_target = str(base_work).replace("'", "''") + phenotype_target = str(phenotype_work).replace("'", "''") + connection = duckdb.connect(config={ + "threads": str(self.threads), + "memory_limit": self.memory_limit, + }) + try: + connection.execute("SET preserve_insertion_order=false") + connection.execute(f"SET temp_directory='{str(work / 'duckdb-tmp').replace(chr(39), chr(39)*2)}'") + connection.execute(f""" +COPY ( + SELECT gene_id, dataset_type, variant_id, + min(try_cast(reported_p_value AS DOUBLE)) AS p_value, + to_json(list_sort(list(DISTINCT phenotype_path_json))) AS phenotype_paths_json, + count(*)::UBIGINT AS association_record_count, + count(*)::UBIGINT AS source_summary_association_count, + count(*)::UBIGINT AS retained_source_summary_count, + list(association_record_id ORDER BY + try_cast(reported_p_value AS DOUBLE) NULLS LAST, + association_record_id)[1] + AS minimum_reported_p_value_association_record_id, + first(source ORDER BY source) AS source, + first(source_display_name ORDER BY source) AS source_display_name, + first(evidence_granularity ORDER BY source) AS evidence_granularity, + first(provider_detail_status ORDER BY source) AS provider_detail_status, + first(provider_detail_reason ORDER BY source) AS provider_detail_reason, + first(missing_fields_json ORDER BY source) AS missing_fields_json, + first(retained_source_summary_artifact ORDER BY source) AS retained_source_summary_artifact, + first(publication_mode ORDER BY source) AS publication_mode, + substr(sha256(gene_id), 1, 16) AS gene_key + FROM read_parquet('{source_glob}', hive_partitioning=true, union_by_name=true) + GROUP BY gene_id, dataset_type, variant_id +) TO '{base_target}' (FORMAT PARQUET, COMPRESSION ZSTD, PARTITION_BY (gene_key), ROW_GROUP_SIZE 50000) +""") + connection.execute(f""" +COPY ( + SELECT gene_id, dataset_type, phenotype_path_json, phenotype_path_key, + count(DISTINCT variant_id)::UBIGINT AS distinct_variant_count, + to_json(list_sort(list(DISTINCT variant_id))) AS variant_ids_json, + min(try_cast(reported_p_value AS DOUBLE)) AS minimum_reported_p_value, + list(association_record_id ORDER BY + try_cast(reported_p_value AS DOUBLE) NULLS LAST, + association_record_id)[1] + AS minimum_reported_p_value_association_record_id, + substr(sha256(gene_id), 1, 16) AS gene_key + FROM read_parquet('{source_glob}', hive_partitioning=true, union_by_name=true) + GROUP BY gene_id, dataset_type, phenotype_path_json, phenotype_path_key +) TO '{phenotype_target}' (FORMAT PARQUET, COMPRESSION ZSTD, PARTITION_BY (gene_key), ROW_GROUP_SIZE 50000) +""") + finally: + connection.close() + base_files = sorted(base_work.rglob("*.parquet")) + phenotype_files = sorted(phenotype_work.rglob("*.parquet")) + base_rows = phenotype_rows = 0 + check = duckdb.connect() + try: + if base_files: + base_rows = int(check.execute( + "SELECT count(*) FROM read_parquet(?, hive_partitioning=true)", + [str(base_work / "gene_key=*" / "*.parquet")], + ).fetchone()[0]) + if phenotype_files: + phenotype_rows = int(check.execute( + "SELECT count(*) FROM read_parquet(?, hive_partitioning=true)", + [str(phenotype_work / "gene_key=*" / "*.parquet")], + ).fetchone()[0]) + finally: + check.close() + final_base.parent.mkdir(parents=True, exist_ok=True) + final_phenotype.parent.mkdir(parents=True, exist_ok=True) + base_work.rename(final_base) + phenotype_work.rename(final_phenotype) + shutil.rmtree(work, ignore_errors=True) + return { + "bucket": bucket, + "base_rows": base_rows, + "phenotype_rows": phenotype_rows, + "base_files": len(base_files), + "phenotype_files": len(phenotype_files), + } + + def build(self) -> dict[str, Any]: + started = time.perf_counter() + manifest = self._validate_source() + checkpoint = json.loads(self.checkpoint_path.read_text()) if self.checkpoint_path.exists() else { + "contract": CONTRACT, + "completed_buckets": {}, + "started_at": _utc_now(), + } + completed = checkpoint.setdefault("completed_buckets", {}) + buckets = [f"{value:02x}" for value in range(256)] + self.work_root.mkdir(parents=True, exist_ok=True) + for index, bucket in enumerate(buckets, start=1): + if bucket in completed: + continue + result = self._rollup_bucket(bucket) + completed[bucket] = result + checkpoint["updated_at"] = _utc_now() + _save_json(self.checkpoint_path, checkpoint) + if index % self.progress_interval == 0: + self.logger.info( + "source-summary rollup progress buckets=%d/256 base_rows=%d phenotype_rows=%d elapsed=%.1fs", + len(completed), + sum(int(row["base_rows"]) for row in completed.values()), + sum(int(row["phenotype_rows"]) for row in completed.values()), + time.perf_counter() - started, + ) + file_rows = [] + for table_name, table_root in ((BASE_TABLE, self.base_root), (PHENOTYPE_TABLE, self.phenotype_root)): + for path in sorted(table_root.rglob("*.parquet")): + file_rows.append({ + "table": table_name, + "path": str(path.relative_to(self.serving_root)), + "size_bytes": path.stat().st_size, + "sha256": _sha256(path), + }) + with self.file_manifest.open("w", encoding="utf-8") as stream: + for row in file_rows: + stream.write(json.dumps(row, sort_keys=True) + "\n") + base_rows = sum(int(row["base_rows"]) for row in completed.values()) + phenotype_rows = sum(int(row["phenotype_rows"]) for row in completed.values()) + tables = manifest.setdefault("tables", {}) + common = { + "contract": CONTRACT, + "partition_key": "gene_id", + "bucket_function": "sha256(gene_id)[0:2]", + "bucket_characters": 2, + "gene_key_function": "sha256(gene_id)[0:16]", + "directory_layout": "nested_gene_key", + } + tables[BASE_TABLE] = { + **common, + "path": str(self.base_root.relative_to(self.serving_root)), + "rows": base_rows, + "files": sum(int(row["base_files"]) for row in completed.values()), + "unit": "gene_dataset_variant", + } + tables[PHENOTYPE_TABLE] = { + **common, + "path": str(self.phenotype_root.relative_to(self.serving_root)), + "rows": phenotype_rows, + "files": sum(int(row["phenotype_files"]) for row in completed.values()), + "unit": "gene_dataset_exact_phenotype_path", + } + manifest["schema_version"] = "2.10.0-rc2" + manifest["source_summary_rollup_built_at"] = _utc_now() + manifest["source_summary_rollup_runtime"] = { + "python": platform.python_version(), + "duckdb": duckdb.__version__, + "threads": self.threads, + "memory_limit": self.memory_limit, + } + manifest_path = self.serving_root / "serving-manifest.json" + _save_json(manifest_path, manifest) + checksum = _sha256(manifest_path) + (self.serving_root / "serving-manifest.json.sha256").write_text( + f"{checksum} serving-manifest.json\n" + ) + result = { + "contract": CONTRACT, + "built_at": manifest["source_summary_rollup_built_at"], + "base_rows": base_rows, + "phenotype_rows": phenotype_rows, + "files": len(file_rows), + "file_manifest": str(self.file_manifest), + "checkpoint": str(self.checkpoint_path), + "serving_manifest": str(manifest_path), + "serving_manifest_sha256": checksum, + "elapsed_seconds": round(time.perf_counter() - started, 3), + } + return result diff --git a/tests/test_association_evidence_v2.py b/tests/test_association_evidence_v2.py new file mode 100644 index 0000000..267c7e2 --- /dev/null +++ b/tests/test_association_evidence_v2.py @@ -0,0 +1,525 @@ +from __future__ import annotations + +import csv +import gzip +import importlib.util +import json +import sys +from pathlib import Path + +import duckdb + +from datahub.association_evidence_v2.builder import EvidenceV2Builder +from datahub.association_evidence_v2.normalization import variation_type_from_alleles +from datahub.association_evidence_v2.serving import AssociationEvidenceV2ServingBuilder + + +HEADER = [ + "MarkerID", + "pval", + "Phenotype", + "Study", + "PMID", + "StudyGenomeBuild", + "dbsnp.rsid", + "dbsnp.dbsnp_build", + "dbsnp.alleles.allele", + "dbsnp.chrom", + "dbsnp.hg19.start", + "dbsnp.hg19.end", + "dbsnp.vartype", + "gnomad_genome.af.af", + "gnomad_genome.af.af_afr", + "snpeff.ann.gene_id", + "snpeff.ann.effect", + "snpeff.ann.putative_impact", + "snpeff.ann.feature_id", + "snpeff.ann.hgvs_p", + "snpeff.ann.protein.length", + "clinvar.rcv.clinical_significance", +] + + +def _write_source(path: Path, delimiter: str, rows: list[list[str]]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with path.open("w", newline="") as stream: + writer = csv.writer(stream, delimiter=delimiter) + writer.writerow(HEADER) + writer.writerows(rows) + + +def _write_variant_index(path: Path, payload: list[dict]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with gzip.open(path, "wt", encoding="utf-8") as stream: + json.dump(payload, stream) + + +def _load_build_script(): + repo_root = Path(__file__).resolve().parents[1] + path = ( + repo_root + / "scripts" + / "dataset_specific_scripts" + / "unified" + / "build_association_evidence_v2.py" + ) + spec = importlib.util.spec_from_file_location("build_association_evidence_v2", path) + assert spec is not None + assert spec.loader is not None + module = importlib.util.module_from_spec(spec) + sys.modules[spec.name] = module + spec.loader.exec_module(module) + return module + + +def test_build_script_writes_serving_checksum_contract(tmp_path: Path) -> None: + module = _load_build_script() + database = tmp_path / "association-evidence-v2.duckdb" + database.write_bytes(b"normalized-sidecar") + digest = module.sha256_file(database) + + checksum_path = module.write_sha256_sidecar(database, digest) + + assert checksum_path.name == "association-evidence-v2.duckdb.sha256" + assert checksum_path.read_text() == f"{digest} {database.name}\n" + + +def _builder(tmp_path: Path) -> EvidenceV2Builder: + cvd_root = tmp_path / "raw" / "cvd" + trait_root = tmp_path / "raw" / "trait" + row_base = [ + "rs1", + "1e-8", + "hypertension", + "Study A", + "PMID: 1", + "b37", + "rs1", + "155", + "['A', 'G']", + "1", + "100", + "100", + "snp", + "0.2", + "0.3", + "GENE1", + "missense_variant", + "MODERATE", + "ENST1", + "p.A1G", + "100", + "Pathogenic", + ] + synonymous = row_base.copy() + synonymous[16] = "synonymous_variant" + synonymous[17] = "LOW" + synonymous[18] = "ENST2" + synonymous[19] = "p.A1A" + synonymous[21] = "Benign" + alternate_protein_length = row_base.copy() + alternate_protein_length[20] = "101" + _write_source( + cvd_root / "hypertension.txt", + ",", + [row_base, synonymous, alternate_protein_length], + ) + _write_source(trait_root / "rr_interval.txt", "\t", [row_base]) + variant_root = tmp_path / "variant_index" + _write_variant_index( + variant_root / "CVD" / "TTN.json.gz", + [ + { + "variant_id": "rsMVP", + "gene_id": "TTN", + "dataset_type": "CVD", + "phenotype": "cardiomyopathy", + "phenotype_path": ["cardiomyopathies", "cardiomyopathy"], + "source": "million_veteran_program", + "sources": ["million_veteran_program"], + "p_value": 1e-9, + "variation_type": "SNP", + } + ], + ) + tree = tmp_path / "phenotype_tree.json" + tree.write_text( + json.dumps( + { + "CVD": {"vascular_diseases": {"hypertensive_disorders": ["hypertension"]}}, + "TRAIT": {"ecg_measures": ["hypertension"]}, + } + ) + ) + return EvidenceV2Builder( + cvd_root=cvd_root, + trait_root=trait_root, + variant_index_root=variant_root, + output_db=tmp_path / "evidence.duckdb", + checkpoint_path=tmp_path / "checkpoint.json", + phenotype_tree=tree, + release_id="test-v2", + command="pytest fixture", + threads=1, + memory_limit="1GB", + smoke_files=1, + progress_interval=1, + ) + + +def test_allele_derived_variation_type_never_assigns_roles() -> None: + result = variation_type_from_alleles(["A", "G"]) + assert result["value"] == "SNV" + assert result["status"] == "derived" + assert "roles remain unresolved" in result["reason"] + assert variation_type_from_alleles(["A"])["status"] == "unresolved" + + +def test_full_separation_and_provenance_contract(tmp_path: Path) -> None: + builder = _builder(tmp_path) + result = builder.run() + assert result["counts"]["provider_records"] == 4 + connection = duckdb.connect(str(tmp_path / "evidence.duckdb"), read_only=True) + + provider_ids = connection.execute( + "SELECT provider_record_id, source_line, source_row_sha256 FROM provider_records ORDER BY provider_record_id" + ).fetchall() + assert len(provider_ids) == len({row[0] for row in provider_ids}) == 4 + assert {row[1] for row in provider_ids} == {2, 3, 4} + assert all(len(row[2]) == 64 for row in provider_ids) + + consequences = connection.execute( + "SELECT consequence, severity_selection, provider_record_count FROM consequence_annotations ORDER BY consequence" + ).fetchall() + assert {row[0] for row in consequences} == {"missense_variant", "synonymous_variant"} + assert sum(row[0] == "missense_variant" for row in consequences) == 2 + assert all(row[1] == "not_performed" for row in consequences) + + assertions = connection.execute( + """SELECT clinical_significance, normalized_term, display_group, + assertion_source, condition_status, evidence_granularity, + provenance_limitation, raw_source_values_json, + provider_record_count + FROM clinical_assertions + ORDER BY normalized_term, assertion_source""" + ).fetchall() + assert [(row[0], row[3]) for row in assertions] == [ + ("Benign", "legacy_cvd_raw"), + ("Pathogenic", "legacy_cvd_raw"), + ("Pathogenic", "legacy_trait_raw"), + ] + assert all(row[0] == row[1] for row in assertions) + assert {row[2] for row in assertions} == {"Benign spectrum", "Pathogenic spectrum"} + assert all(row[4] == "unavailable" for row in assertions) + assert all(row[5] == "legacy_variant_level_annotation" for row in assertions) + assert all("may not refer to the selected HBP phenotype" in row[6] for row in assertions) + assert all(json.loads(row[7]) for row in assertions) + assert all(row[8] > 0 for row in assertions) + assert connection.execute( + """SELECT count(*)=count(DISTINCT variant_id || chr(31) || normalized_term || chr(31) || assertion_source) + FROM clinical_assertions""" + ).fetchone() == (True,) + + mvp = connection.execute( + """ + SELECT provider_detail_status, provider_detail_reason, missing_fields_json, + logical_artifact, publication_mode + FROM source_summary_artifacts WHERE source='million_veteran_program' + """ + ).fetchone() + assert mvp[0] == "not_applicable" + assert mvp[1] is None + assert "effect_allele" in json.loads(mvp[2]) + assert mvp[3].endswith("TTN.json.gz") + assert mvp[4] == "gene_scoped_on_demand" + assert connection.execute( + "SELECT count(*) FROM association_records WHERE source='million_veteran_program'" + ).fetchone() == (0,) + + recoverable = connection.execute( + """ + SELECT effect_allele_status, effect_size_status, standard_error_status, + sample_size_status, ancestry_status, fine_mapping_status, + missing_fields_json + FROM association_records + WHERE source='legacy_cvd_raw' + ORDER BY association_record_id + LIMIT 1 + """ + ).fetchone() + assert recoverable[:6] == ("unavailable",) * 6 + assert set(json.loads(recoverable[6])) == { + "effect_allele", + "effect_size", + "standard_error", + "sample_size", + "ancestry", + "fine_mapping", + } + + summary = connection.execute( + """ + SELECT minimum_reported_p_value, minimum_reported_p_value_association_record_id, + consequence_conflict, clinical_significance_conflict + FROM variant_phenotype_summaries + WHERE gene_id='GENE1' AND dataset_type='CVD' AND variant_id='rs1' + """ + ).fetchone() + assert summary[0] == 1e-8 + supplier = connection.execute( + "SELECT reported_p_value FROM association_records WHERE association_record_id=?", + [summary[1]], + ).fetchone() + assert supplier == (1e-8,) + assert summary[2] is True + assert summary[3] is True + + qa = connection.execute("SELECT status FROM qa_results").fetchall() + assert qa and {row[0] for row in qa} == {"passed"} + connection.close() + + +def test_resume_keeps_provider_ids_and_counts_stable(tmp_path: Path) -> None: + builder = _builder(tmp_path) + first = builder.run() + connection = duckdb.connect(str(tmp_path / "evidence.duckdb"), read_only=True) + first_ids = connection.execute( + "SELECT provider_record_id FROM provider_records ORDER BY provider_record_id" + ).fetchall() + connection.close() + + second = _builder(tmp_path).run() + connection = duckdb.connect(str(tmp_path / "evidence.duckdb"), read_only=True) + second_ids = connection.execute( + "SELECT provider_record_id FROM provider_records ORDER BY provider_record_id" + ).fetchall() + connection.close() + assert first["counts"] == second["counts"] + assert first_ids == second_ids + + +def test_clean_rebuild_keeps_provider_ids_and_source_lines_stable(tmp_path: Path) -> None: + first_root = tmp_path / "first" + second_root = tmp_path / "second" + first = _builder(first_root) + second = _builder(second_root) + + # The fixture factory writes identical logical source files below each root. + first.run() + second.run() + first_connection = duckdb.connect(str(first_root / "evidence.duckdb"), read_only=True) + second_connection = duckdb.connect(str(second_root / "evidence.duckdb"), read_only=True) + query = """ + SELECT provider_record_id, source_file_logical, source_line, source_row_sha256 + FROM provider_records ORDER BY provider_record_id + """ + assert first_connection.execute(query).fetchall() == second_connection.execute(query).fetchall() + first_connection.close() + second_connection.close() + + +def test_variant_builder_streams_past_insert_batch_boundary(tmp_path: Path) -> None: + cvd_root = tmp_path / "raw" / "cvd" + trait_root = tmp_path / "raw" / "trait" + rows = [] + for index in range(10005): + rows.append( + [ + f"rs{index}", + "1e-8", + "hypertension", + "Study A", + "1", + "b37", + f"rs{index}", + "155", + "['A', 'G']", + "1", + str(index), + str(index), + "snp", + "", + "", + "GENE1", + "missense_variant", + "MODERATE", + "ENST1", + "p.A1G", + "100", + "", + ] + ) + _write_source(cvd_root / "hypertension.txt", ",", rows) + _write_source(trait_root / "rr_interval.txt", "\t", [rows[0]]) + variant_root = tmp_path / "variant_index" + (variant_root / "CVD").mkdir(parents=True) + (variant_root / "TRAIT").mkdir(parents=True) + tree = tmp_path / "phenotype_tree.json" + tree.write_text( + json.dumps( + { + "CVD": {"vascular_diseases": {"hypertensive_disorders": ["hypertension"]}}, + "TRAIT": {"ecg_measures": ["hypertension"]}, + } + ) + ) + builder = EvidenceV2Builder( + cvd_root=cvd_root, + trait_root=trait_root, + variant_index_root=variant_root, + output_db=tmp_path / "evidence.duckdb", + checkpoint_path=tmp_path / "checkpoint.json", + phenotype_tree=tree, + release_id="batch-boundary-test", + command="pytest fixture", + threads=1, + memory_limit="1GB", + progress_interval=1, + ) + result = builder.run() + assert result["counts"]["variants"] == 10005 + connection = duckdb.connect(str(tmp_path / "evidence.duckdb"), read_only=True) + assert connection.execute( + "SELECT count(*) FROM variants WHERE variation_type='SNV'" + ).fetchone() == (10005,) + connection.close() + + + +def test_partitioned_serving_preserves_provider_rows_and_associations(tmp_path: Path) -> None: + builder = _builder(tmp_path) + result = builder.run() + build_script = _load_build_script() + source_db = tmp_path / "evidence.duckdb" + source_digest = build_script.sha256_file(source_db) + build_script.write_sha256_sidecar(source_db, source_digest) + coarse_root = tmp_path / "serving-coarse" + coarse_manifest = AssociationEvidenceV2ServingBuilder( + source_db=tmp_path / "evidence.duckdb", + intermediate_root=tmp_path / "duckdb_tmp", + output_root=coarse_root, + memory_limit="1GB", + threads=1, + variant_bucket_characters=2, + progress_interval=1, + ).build() + passthrough_names = ( + "source_summary_associations_by_gene", + "source_summary_association_phenotype_counts_by_gene", + "source_summary_association_base_by_gene", + ) + for logical_name in passthrough_names: + table_root = coarse_root / "tables" / logical_name / "gene_bucket=00" + table_root.mkdir(parents=True) + (table_root / "part-0.parquet").write_bytes(b"fixture") + coarse_manifest["tables"][logical_name] = { + "path": f"tables/{logical_name}", + "partition_key": "gene_id", + "bucket_function": "sha256(gene_id)[0:2]", + "bucket_characters": 2, + "rows": 1, + "files": 1, + } + (coarse_root / "serving-manifest.json").write_text( + json.dumps(coarse_manifest, indent=2, sort_keys=True) + "\n" + ) + serving_root = tmp_path / "serving" + manifest = AssociationEvidenceV2ServingBuilder( + source_db=tmp_path / "evidence.duckdb", + intermediate_root=tmp_path / "duckdb_tmp", + output_root=serving_root, + memory_limit="1GB", + threads=1, + variant_bucket_characters=3, + progress_interval=1, + coarse_serving_root=coarse_root, + ).build() + finer_root = tmp_path / "serving-finer" + finer_manifest = AssociationEvidenceV2ServingBuilder( + source_db=tmp_path / "evidence.duckdb", + intermediate_root=tmp_path / "duckdb_tmp", + output_root=finer_root, + memory_limit="1GB", + threads=1, + variant_bucket_characters=4, + progress_interval=1, + coarse_serving_root=serving_root, + ).build() + for logical_name in passthrough_names: + assert manifest["tables"][logical_name]["rows"] == 1 + assert finer_manifest["tables"][logical_name]["rows"] == 1 + coarse_file = next((coarse_root / "tables" / logical_name).rglob("*.parquet")) + serving_file = next((serving_root / "tables" / logical_name).rglob("*.parquet")) + finer_file = next((finer_root / "tables" / logical_name).rglob("*.parquet")) + assert coarse_file.stat().st_ino == serving_file.stat().st_ino + assert serving_file.stat().st_ino == finer_file.stat().st_ino + assert manifest["contract"] == "association_evidence_v2_partitioned_serving" + assert manifest["source_db_sha256"] == source_digest + assert manifest["tables"]["provider_records"]["normalized_source_rows"] == result["counts"]["provider_records"] + assert manifest["tables"]["provider_records"]["excluded_unkeyed_rows"] == 0 + assert manifest["tables"]["provider_records"]["rows"] == result["counts"]["provider_records"] + assert manifest["tables"]["association_records"]["rows"] == result["counts"]["association_records"] + assert ( + manifest["tables"]["variant_phenotype_summaries"]["rows"] + == result["counts"]["variant_phenotype_summaries"] + ) + assert manifest["schema_version"] == "2.6.0-rc1" + assert manifest["variant_bucket_count"] == 4096 + assert finer_manifest["variant_bucket_count"] == 65536 + assert manifest["tables"]["provider_records"]["bucket_characters"] == 3 + assert manifest["tables"]["association_records"]["bucket_characters"] == 3 + assert manifest["tables"]["variant_phenotype_summaries"]["bucket_characters"] == 3 + for logical_name in ( + "consequence_annotations", + "clinical_assertions", + "population_observations", + ): + assert manifest["tables"][logical_name]["bucket_characters"] == 3 + assert manifest["tables"][logical_name]["directory_layout"] == "nested_coarse_fine" + assert manifest["tables"]["association_records"]["directory_layout"] == "nested_coarse_fine" + assert any( + path.name.startswith("coarse_bucket=") + for path in (serving_root / "tables" / "association_records").iterdir() + ) + for logical_name in ( + "provider_records", + "association_records", + "variant_phenotype_summaries", + "consequence_annotations", + "clinical_assertions", + "population_observations", + ): + assert finer_manifest["tables"][logical_name]["bucket_characters"] == 4 + assert finer_manifest["tables"][logical_name]["directory_layout"] == "nested_coarse_fine" + assert ( + manifest["tables"]["consequence_annotations_by_gene"]["rows"] + == result["counts"]["consequence_annotations"] + ) + connection = duckdb.connect() + provider_rows = connection.execute( + f"SELECT provider_record_id, source_line, source_row_sha256 " + f"FROM read_parquet('{serving_root / 'tables/provider_records/**/*.parquet'}') " + "ORDER BY provider_record_id" + ).fetchall() + association_rows = connection.execute( + f"SELECT association_record_id, variant_id, reported_p_value " + f"FROM read_parquet('{serving_root / 'tables/association_records/**/*.parquet'}') " + "ORDER BY association_record_id" + ).fetchall() + finer_association_rows = connection.execute( + f"SELECT association_record_id, variant_id, reported_p_value " + f"FROM read_parquet('{finer_root / 'tables/association_records/**/*.parquet'}') " + "ORDER BY association_record_id" + ).fetchall() + gene_consequences = connection.execute( + f"SELECT gene_id, variant_id, consequence " + f"FROM read_parquet('{serving_root / 'tables/consequence_annotations_by_gene/**/*.parquet'}') " + "ORDER BY gene_id, variant_id, consequence" + ).fetchall() + assert len(provider_rows) == result["counts"]["provider_records"] + assert all(row[1] >= 2 and len(row[2]) == 64 for row in provider_rows) + assert len(association_rows) == result["counts"]["association_records"] + assert finer_association_rows == association_rows + assert len(gene_consequences) == result["counts"]["consequence_annotations"] + assert {row[0] for row in gene_consequences} == {"GENE1"} + connection.close() diff --git a/tests/test_association_evidence_v2_audit.py b/tests/test_association_evidence_v2_audit.py new file mode 100644 index 0000000..f5c4327 --- /dev/null +++ b/tests/test_association_evidence_v2_audit.py @@ -0,0 +1,250 @@ +"""Contract tests for the resumable association evidence v2 release audit.""" + +from __future__ import annotations + +import gzip +import hashlib +import importlib.util +import json +import logging +from pathlib import Path + +import duckdb + + +def _module(): + path = Path(__file__).parents[1] / "scripts" / "dataset_specific_scripts" / "unified" / "audit_association_evidence_v2_release.py" + spec = importlib.util.spec_from_file_location("evidence_v2_audit", path) + module = importlib.util.module_from_spec(spec) + assert spec.loader is not None + spec.loader.exec_module(module) + return module + + +def test_count_intermediate_round_trip() -> None: + module = _module() + category = {("HMGCR", "CVD", "vc", "snv"): 379} + totals = {("HMGCR", "CVD"): {"v2_distinct_variants": 379, "association_record_count": 400}} + encoded = module.encode_count_bundle(category, totals) + observed_category, observed_totals = module.decode_count_bundle(encoded, rich_totals=True) + assert observed_category == category + assert observed_totals == totals + + +def test_checkpoint_is_configuration_bound(tmp_path: Path) -> None: + module = _module() + checkpoint = tmp_path / "checkpoint.json" + module.save_checkpoint(checkpoint, "config-a", "v2_counts_complete", {"rows": 3}) + assert module.load_checkpoint(checkpoint, "config-a")["details"]["rows"] == 3 + try: + module.load_checkpoint(checkpoint, "config-b") + except RuntimeError as exc: + assert "does not match" in str(exc) + else: + raise AssertionError("mismatched audit configuration was accepted") + + +def test_category_difference_reason_distinguishes_missing_gene_and_new_term() -> None: + module = _module() + totals = {("HMGCR", "CVD"): {"v2_distinct_variants": 379}} + missing = module.category_difference_reason( + gene="5S_rRNA", dataset_type="CVD", dimension="vc", + v1_count=40, v2_count=0, v2_totals=totals, + ) + assert "No v2 gene/dataset summary" in missing + added = module.category_difference_reason( + gene="HMGCR", dataset_type="CVD", dimension="msc", + v1_count=0, v2_count=12, v2_totals=totals, + ) + assert "not exposed by the v1 representative" in added + + +def test_streaming_array_and_resumable_membership_registry(tmp_path: Path) -> None: + module = _module() + artifact_root = tmp_path / "artifacts" + artifact = artifact_root / "variant_index" / "TRAIT" / "HMGCR.json.gz" + artifact.parent.mkdir(parents=True) + with gzip.open(artifact, "wt", encoding="utf-8") as stream: + json.dump([ + { + "variant_id": "rs1", + "source": "million_veteran_program", + "sources": ["million_veteran_program"], + }, + { + "variant_id": "rs1", + "source": "million_veteran_program", + "sources": ["million_veteran_program"], + }, + { + "variant_id": "rs2", + "source": "other", + "sources": ["other"], + }, + ], stream) + assert [row["variant_id"] for row in module.iter_json_array(artifact, chunk_size=7)] == [ + "rs1", "rs1", "rs2" + ] + + registry_path = tmp_path / "registry.duckdb" + registry = duckdb.connect(str(registry_path)) + registry.execute( + """CREATE TABLE source_summary_artifacts ( + source VARCHAR, gene_id VARCHAR, dataset_type VARCHAR, + logical_artifact VARCHAR, size_bytes UBIGINT, modified_at_ns UBIGINT, + provider_detail_status VARCHAR, publication_mode VARCHAR); + CREATE TABLE variant_phenotype_summaries ( + gene_id VARCHAR, dataset_type VARCHAR, variant_id VARCHAR, + source_observation_count UBIGINT, association_record_count UBIGINT, + provider_association_record_count UBIGINT); + CREATE TABLE variants (variant_id VARCHAR, variation_type VARCHAR); + CREATE TABLE consequence_annotations ( + variant_id VARCHAR, gene_id VARCHAR, consequence VARCHAR); + CREATE TABLE clinical_assertions ( + variant_id VARCHAR, clinical_significance VARCHAR); + INSERT INTO variants VALUES ('rs1', 'SNV'); + INSERT INTO consequence_annotations VALUES ('rs1', 'HMGCR', 'missense_variant'); + INSERT INTO clinical_assertions VALUES ('rs1', 'Pathogenic');""" + ) + stat = artifact.stat() + registry.execute( + "INSERT INTO source_summary_artifacts VALUES (?, ?, ?, ?, ?, ?, ?, ?)", + [ + "million_veteran_program", "HMGCR", "TRAIT", + "variant_index/TRAIT/HMGCR.json.gz", stat.st_size, stat.st_mtime_ns, + "unavailable", "gene_scoped_on_demand", + ], + ) + runtime_path = tmp_path / "runtime.duckdb" + logger = logging.getLogger("test.association_evidence_v2.audit") + first = module.materialize_compact_membership( + registry, + artifact_root=artifact_root, + runtime_db_path=runtime_path, + logger=logger, + progress_interval=1, + batch_rows=1, + ) + second = module.materialize_compact_membership( + registry, + artifact_root=artifact_root, + runtime_db_path=runtime_path, + logger=logger, + progress_interval=1, + batch_rows=1, + ) + registry.close() + read_only_registry = duckdb.connect(str(registry_path), read_only=True) + category, totals = module.v2_counts( + read_only_registry, + runtime_db_path=runtime_path, + v2_db_path=registry_path, + ) + read_only_registry.close() + assert first == second + assert first["artifacts"] == 1 + assert first["selected_source_rows"] == 2 + assert first["distinct_gene_dataset_variants"] == 1 + assert category[("HMGCR", "TRAIT", "vc", "snv")] == 1 + assert category[("HMGCR", "TRAIT", "msc", "missense variant")] == 1 + assert category[("HMGCR", "TRAIT", "cs", "pathogenic")] == 1 + assert totals[("HMGCR", "TRAIT")]["v2_distinct_variants"] == 1 + assert totals[("HMGCR", "TRAIT")]["association_record_count"] == 0 + + +def test_artifact_reader_preserves_selected_count_and_deduplicates_membership( + tmp_path: Path, +) -> None: + module = _module() + artifact = tmp_path / "HMGCR.json.gz" + with gzip.open(artifact, "wt", encoding="utf-8") as stream: + json.dump([ + {"variant_id": "rs1", "sources": ["million_veteran_program"]}, + {"variant_id": "rs1", "sources": ["million_veteran_program"]}, + {"variant_id": "rs2", "source": "other"}, + ], stream) + stat = artifact.stat() + result = module._read_membership_artifact(( + "million_veteran_program", "HMGCR", "TRAIT", "logical/HMGCR.json.gz", + stat.st_size, stat.st_mtime_ns, str(artifact), + )) + assert result["selected_rows"] == 2 + assert result["membership_rows"] == [("HMGCR", "TRAIT", "rs1")] + + + +def test_indexed_membership_reuses_completed_serving_contract(tmp_path: Path) -> None: + module = _module() + serving = tmp_path / "serving" + table_root = serving / "tables" / "source_summary_associations_by_gene" + connection = duckdb.connect() + manifest_rows = [] + for bucket in range(256): + target = table_root / f"gene_bucket={bucket:02x}" / "gene_key=fixture" + target.mkdir(parents=True, exist_ok=True) + connection.execute( + """COPY (SELECT 'source-summary:1'::VARCHAR AS source_summary_id, + 'million_veteran_program'::VARCHAR AS source, + 'CVD'::VARCHAR AS dataset_type, + 'HMGCR'::VARCHAR AS gene_id, + 'rs1'::VARCHAR AS variant_id) + TO ? (FORMAT PARQUET)""", + [str(target / f"part-{bucket:02x}.parquet")], + ) + manifest_rows.append(json.dumps({ + "logical_artifact": f"variant_index/CVD/HMGCR-{bucket:02x}.json.gz", + "source": "million_veteran_program", + "gene_id": "HMGCR", + "dataset_type": "CVD", + "selected_source_rows": 1, + })) + connection.close() + (serving / "source-summary-index-files.jsonl").write_text( + "\n".join(manifest_rows) + "\n" + ) + manifest = { + "schema_version": "2.9.0-rc2", + "tables": {"source_summary_associations_by_gene": { + "path": "tables/source_summary_associations_by_gene", + "contract": "source_summary_association_index_v2", + "files": 256, + "source_artifacts": 256, + "selected_source_rows": 256, + }}, + } + manifest_path = serving / "serving-manifest.json" + manifest_path.write_text(json.dumps(manifest)) + checksum = hashlib.sha256(manifest_path.read_bytes()).hexdigest() + (serving / "serving-manifest.json.sha256").write_text( + f"{checksum} serving-manifest.json\n" + ) + for schema_version in ("2.9.0-rc2", "2.10.0-rc2"): + manifest["schema_version"] = schema_version + if schema_version == "2.10.0-rc2": + manifest["tables"]["source_summary_association_base_by_gene"] = { + "path": "tables/source_summary_associations_by_gene", + "contract": "source_summary_association_rollup_v2", + "files": 256, + "rows": 256, + } + (serving / "source-summary-rollup-files.jsonl").write_text("{}\n") + else: + manifest["tables"].pop( + "source_summary_association_base_by_gene", None + ) + manifest_path.write_text(json.dumps(manifest)) + checksum = hashlib.sha256(manifest_path.read_bytes()).hexdigest() + (serving / "serving-manifest.json.sha256").write_text( + f"{checksum} serving-manifest.json\n" + ) + stats = module.materialize_indexed_compact_membership( + serving_root=serving, + runtime_db_path=tmp_path / f"runtime-indexed-{schema_version}.duckdb", + logger=logging.getLogger("test.association_evidence_v2.indexed-audit"), + ) + assert stats == { + "artifacts": 256, + "distinct_gene_dataset_variants": 1, + "selected_source_rows": 256, + "indexed_gene_buckets": 256, + } diff --git a/tests/test_association_evidence_v2_source_summary_index.py b/tests/test_association_evidence_v2_source_summary_index.py new file mode 100644 index 0000000..17b8664 --- /dev/null +++ b/tests/test_association_evidence_v2_source_summary_index.py @@ -0,0 +1,100 @@ +from __future__ import annotations + +import gzip +import json +from pathlib import Path + +import duckdb + +from datahub.association_evidence_v2.source_summary_index import SourceSummaryIndexBuilder + + +def test_index_is_gene_keyed_and_does_not_create_provider_rows(tmp_path: Path) -> None: + artifact = tmp_path / "artifacts/variant_index/CVD/TTN.json.gz" + artifact.parent.mkdir(parents=True) + payload = { + "variant_id": "rsMVP", + "phenotype": "cardiomyopathy", + "phenotype_path": ["cardiomyopathies", "cardiomyopathy"], + "source": "million_veteran_program", + "sources": ["million_veteran_program"], + "p_value": 1e-9, + "ancestry": {"European": 1}, + "label_key": "disease", + "metadata": { + "phenotype_keys": ["Phe_425_11"], + "source_files": [ + "/N/scratch/kvand/aws_import/aggregated_phenotypes/cardiomyopathy.csv.gz" + ], + }, + } + with gzip.open(artifact, "wt", encoding="utf-8") as stream: + json.dump([payload, payload, {**payload, "sources": ["legacy_cvd_raw"]}], stream) + stat = artifact.stat() + source_db = tmp_path / "source.duckdb" + connection = duckdb.connect(str(source_db)) + connection.execute( + """CREATE TABLE source_summary_artifacts ( + source_summary_artifact_id VARCHAR, source VARCHAR, dataset_type VARCHAR, + gene_id VARCHAR, logical_artifact VARCHAR, file_format VARCHAR, + size_bytes UBIGINT, modified_at_ns UBIGINT, + provider_detail_status VARCHAR, provider_detail_reason VARCHAR, + missing_fields_json JSON, content_contract VARCHAR, publication_mode VARCHAR + )""" + ) + connection.execute( + "INSERT INTO source_summary_artifacts VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?)", + [ + "artifact-1", "million_veteran_program", "CVD", "TTN", + "variant_index/CVD/TTN.json.gz", "json.gz", stat.st_size, + stat.st_mtime_ns, "not_applicable", None, + json.dumps(["effect_allele", "sample_size", "ancestry"]), + "retained_compact_variant_index_v1", "gene_scoped_on_demand", + ], + ) + connection.close() + serving_root = tmp_path / "serving" + serving_root.mkdir() + (serving_root / "serving-manifest.json").write_text( + json.dumps({"contract": "association_evidence_v2_partitioned_serving", "tables": {}}) + ) + + result = SourceSummaryIndexBuilder( + source_db=source_db, + artifact_root=tmp_path / "artifacts", + serving_root=serving_root, + workers=1, + progress_interval=1, + ).build() + + assert result["selected_source_rows"] == 2 + assert result["unique_summary_rows"] == 1 + manifest = json.loads((serving_root / "serving-manifest.json").read_text()) + table = manifest["tables"]["source_summary_associations_by_gene"] + assert table["provider_detail_status"] == "not_applicable" + assert table["record_kind"] == "source_summary_association" + assert table["evidence_granularity"] == "source_summary" + rows = duckdb.connect().execute( + f"SELECT source, source_display_name, dataset_id, variant_id, " + f"phenotype_slug, phenotype_kind, record_kind, evidence_granularity, " + f"provider_detail_status, reported_p_value, retained_source_summary_artifact, " + f"ancestry_json, missing_fields_json, metadata_json, " + f"retained_source_summary_json FROM read_parquet(" + f"'{serving_root / 'tables/source_summary_associations_by_gene/**/*.parquet'}')" + ).fetchall() + assert rows == [ + ( + "million_veteran_program", "MVP", "hbp_mvp_association", "rsMVP", + "cardiomyopathy", "disease", "source_summary_association", + "source_summary", "not_applicable", "1e-09", + "variant_index/CVD/TTN.json.gz", '{"European":1}', + '["effect_allele","sample_size"]', + '{"phenotype_keys":["Phe_425_11"]}', + '{"variant_id":"rsMVP","phenotype":"cardiomyopathy",' + '"phenotype_path":["cardiomyopathies","cardiomyopathy"],' + '"source":"million_veteran_program",' + '"sources":["million_veteran_program"],"p_value":1e-09,' + '"ancestry":{"European":1},"label_key":"disease",' + '"metadata":{"phenotype_keys":["Phe_425_11"]}}', + ) + ] diff --git a/tests/test_association_evidence_v2_source_summary_rollup.py b/tests/test_association_evidence_v2_source_summary_rollup.py new file mode 100644 index 0000000..029037c --- /dev/null +++ b/tests/test_association_evidence_v2_source_summary_rollup.py @@ -0,0 +1,64 @@ +"""Tests for evidence-v2 source-summary default rollups.""" + +from __future__ import annotations + +import hashlib +import json +from pathlib import Path + +import duckdb + +from datahub.association_evidence_v2.source_summary_rollup import ( + BASE_TABLE, + PHENOTYPE_TABLE, + SourceSummaryRollupBuilder, +) + + +def test_rollup_preserves_variant_and_exact_phenotype_counts(tmp_path: Path) -> None: + serving = tmp_path / "serving" + source = serving / "tables" / "source_summary_associations_by_gene" / "gene_bucket=00" / "gene_key=fixture" + source.mkdir(parents=True) + connection = duckdb.connect() + connection.execute( + """COPY (SELECT * FROM (VALUES + ('a1','mvp','MVP','source_summary','CVD','HMGCR','rs1','[\"vascular\",\"hypertension\"]','vascular > hypertension','1e-6','not_applicable',NULL,'[\"effect_allele\"]','variant_index/CVD/HMGCR.json.gz','gene_scoped_on_demand'), + ('a2','mvp','MVP','source_summary','CVD','HMGCR','rs1','[\"vascular\",\"cad\"]','vascular > cad','1e-9','not_applicable',NULL,'[\"effect_allele\"]','variant_index/CVD/HMGCR.json.gz','gene_scoped_on_demand'), + ('a3','mvp','MVP','source_summary','CVD','HMGCR','rs2','[\"vascular\",\"hypertension\"]','vascular > hypertension','0.2','not_applicable',NULL,'[\"effect_allele\"]','variant_index/CVD/HMGCR.json.gz','gene_scoped_on_demand') + ) t(association_record_id,source,source_display_name,evidence_granularity,dataset_type,gene_id,variant_id,phenotype_path_json,phenotype_path_key,reported_p_value,provider_detail_status,provider_detail_reason,missing_fields_json,retained_source_summary_artifact,publication_mode)) + TO ? (FORMAT PARQUET)""", + [str(source / "part.parquet")], + ) + connection.close() + manifest = { + "contract": "association_evidence_v2_partitioned_serving", + "schema_version": "2.9.0-rc2", + "tables": {"source_summary_associations_by_gene": { + "contract": "source_summary_association_index_v2", + "files": 1, + "source_artifacts": 1, + }}, + } + manifest_path = serving / "serving-manifest.json" + manifest_path.write_text(json.dumps(manifest)) + checksum = hashlib.sha256(manifest_path.read_bytes()).hexdigest() + (serving / "serving-manifest.json.sha256").write_text(f"{checksum} serving-manifest.json\n") + builder = SourceSummaryRollupBuilder(serving_root=serving) + result = builder._rollup_bucket("00") + assert result["base_rows"] == 2 + assert result["phenotype_rows"] == 2 + check = duckdb.connect() + base = check.execute( + "SELECT variant_id,p_value,source_summary_association_count,minimum_reported_p_value_association_record_id FROM read_parquet(?) ORDER BY variant_id", + [str(serving / "tables" / BASE_TABLE / "gene_bucket=00" / "gene_key=*" / "*.parquet")], + ).fetchall() + phenotypes = check.execute( + "SELECT phenotype_path_key,distinct_variant_count,variant_ids_json FROM read_parquet(?) ORDER BY phenotype_path_key", + [str(serving / "tables" / PHENOTYPE_TABLE / "gene_bucket=00" / "gene_key=*" / "*.parquet")], + ).fetchall() + check.close() + assert base == [("rs1", 1e-9, 2, "a2"), ("rs2", 0.2, 1, "a3")] + assert phenotypes == [ + ("vascular > cad", 1, '["rs1"]'), + ("vascular > hypertension", 2, '["rs1","rs2"]'), + ] diff --git a/tests/test_association_evidence_v2_variant_lookup.py b/tests/test_association_evidence_v2_variant_lookup.py new file mode 100644 index 0000000..6ecca7c --- /dev/null +++ b/tests/test_association_evidence_v2_variant_lookup.py @@ -0,0 +1,168 @@ +from __future__ import annotations + +import importlib.util +import json +import sys +from argparse import Namespace +from pathlib import Path + +import duckdb + + +def load_builder(): + path = ( + Path(__file__).resolve().parents[1] + / "scripts" + / "dataset_specific_scripts" + / "unified" + / "build_association_evidence_v2_variant_lookup.py" + ) + spec = importlib.util.spec_from_file_location("build_v2_variant_lookup", path) + assert spec and spec.loader + module = importlib.util.module_from_spec(spec) + sys.modules[spec.name] = module + spec.loader.exec_module(module) + return module + + +def load_duckdb_builder(): + path = ( + Path(__file__).resolve().parents[1] + / "scripts" + / "dataset_specific_scripts" + / "unified" + / "build_association_evidence_v2_variant_lookup_duckdb.py" + ) + spec = importlib.util.spec_from_file_location("build_v2_variant_lookup_duckdb", path) + assert spec and spec.loader + module = importlib.util.module_from_spec(spec) + sys.modules[spec.name] = module + spec.loader.exec_module(module) + return module + + +def create_source(path: Path) -> None: + connection = duckdb.connect(str(path)) + connection.execute( + """ + CREATE TABLE variants ( + variant_id VARCHAR, + normalized_allele_contexts_json JSON, + coordinate_contexts_json JSON, + source_reported_variation_types_json JSON, + variation_type VARCHAR, + variation_type_status VARCHAR, + variation_type_method VARCHAR, + variation_type_reason VARCHAR, + ref_allele_status VARCHAR, + alt_allele_status VARCHAR, + effect_allele_status VARCHAR, + multiallelic_observed BOOLEAN, + build_or_position_conflict BOOLEAN + ) + """ + ) + connection.executemany( + "INSERT INTO variants VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?)", + [ + ("rs1", "[]", "[]", "[]", "SNV", "derived", "test", None, + "unresolved", "unresolved", "unresolved", False, False), + ("rs2", "[]", "[]", "[]", None, "unresolved", "test", "missing", + "unresolved", "unresolved", "unresolved", False, False), + ], + ) + connection.close() + + +def arguments(source: Path, output: Path) -> Namespace: + return Namespace( + source_db=source, + source_sha256="source-digest", + output_root=output, + bucket_characters=3, + threads=1, + memory_limit="1GB", + progress_interval=1, + smoke_limit=None, + reset=False, + verbose=False, + ) + + +def test_variant_lookup_is_partitioned_checksummed_and_resumable(tmp_path: Path) -> None: + module = load_builder() + source = tmp_path / "source.duckdb" + output = tmp_path / "lookup" + create_source(source) + + first = module.build(arguments(source, output)) + second = module.build(arguments(source, output)) + + assert first == second + assert first["tables"]["variants"]["rows"] == 2 + assert first["tables"]["variants"]["files"] == 2 + assert len((output / "checksums.sha256").read_text().splitlines()) == 3 + checkpoint = json.loads((output / "variant-lookup-checkpoint.json").read_text()) + assert checkpoint["complete"] is True + + rows = duckdb.connect().execute( + "SELECT variant_id, variation_type FROM read_parquet(?, hive_partitioning=true)", + [str(output / "tables" / "variants" / "**" / "*.parquet")], + ).fetchall() + assert sorted(rows) == [("rs1", "SNV"), ("rs2", None)] + + +def test_variant_lookup_rejects_checkpoint_configuration_change(tmp_path: Path) -> None: + module = load_builder() + source = tmp_path / "source.duckdb" + output = tmp_path / "lookup" + create_source(source) + module.build(arguments(source, output)) + + changed = arguments(source, output) + changed.source_sha256 = "different" + try: + module.build(changed) + except RuntimeError as error: + assert "checkpoint configuration" in str(error) + else: + raise AssertionError("Configuration mismatch was accepted") + + +def test_indexed_duckdb_lookup_is_unique_checksummed_and_resumable(tmp_path: Path) -> None: + partition_builder = load_builder() + duckdb_builder = load_duckdb_builder() + source = tmp_path / "source.duckdb" + parquet_output = tmp_path / "lookup-parquet" + create_source(source) + partition_builder.build(arguments(source, parquet_output)) + output_db = tmp_path / "lookup.duckdb" + args = Namespace( + parquet_root=parquet_output / "tables" / "variants", + parquet_manifest_sha256="parquet-manifest-digest", + output_db=output_db, + expected_rows=2, + threads=1, + memory_limit="1GB", + progress_interval=1, + smoke_limit=None, + reset=False, + verbose=False, + ) + + first = duckdb_builder.build(args) + second = duckdb_builder.build(args) + + assert first == second + assert first["rows"] == 2 + assert first["duplicate_variant_ids"] == 0 + assert first["database_sha256"] == duckdb_builder.sha256_file(output_db) + connection = duckdb.connect(str(output_db), read_only=True) + assert connection.execute( + "SELECT variation_type FROM variants WHERE variant_id='rs1'" + ).fetchone() == ("SNV",) + indexes = connection.execute( + "SELECT index_name FROM duckdb_indexes() WHERE table_name='variants'" + ).fetchall() + connection.close() + assert indexes == [("idx_variants_variant_id",)] diff --git a/tests/test_clinical_classification.py b/tests/test_clinical_classification.py new file mode 100644 index 0000000..597c485 --- /dev/null +++ b/tests/test_clinical_classification.py @@ -0,0 +1,59 @@ +from __future__ import annotations + +import pytest + +from datahub.association_evidence_v2.clinical_classification import ( + APPROVED_TERMS, + display_group_for_term, + normalized_clinical_terms, + parse_clinical_classification, +) + + +@pytest.mark.parametrize( + ("raw_value", "expected"), + [ + ("Pathogenic", ["Pathogenic"]), + ("['Benign', 'Likely benign']", ["Benign", "Likely benign"]), + ('["Benign", "association"]', ["Benign", "association"]), + ("[Benign, Likely benign, Benign]", ["Benign", "Likely benign"]), + ("Pathogenic, other", ["Pathogenic", "other"]), + ("['Pathogenic, other', 'Benign']", ["Pathogenic", "other", "Benign"]), + ( + "['Conflicting interpretations of pathogenicity, Affects, other']", + ["Conflicting interpretations of pathogenicity", "Affects", "other"], + ), + ("association; protective", ["association", "protective"]), + ("Benign / Likely benign", ["Benign/Likely benign"]), + (" affects. ", ["Affects"]), + ("likely_benign", ["Likely benign"]), + ], +) +def test_audited_clinical_classification_parser(raw_value, expected): + assert normalized_clinical_terms(raw_value) == expected + + +def test_repeated_terms_are_removed_only_from_display_projection(): + parsed = parse_clinical_classification("[Benign, Benign, Likely benign]") + assert parsed.terms == ("Benign", "Benign", "Likely benign") + assert parsed.display_terms == ("Benign", "Likely benign") + assert parsed.duplicate_terms == {"Benign": 2} + assert parsed.raw_value == "[Benign, Benign, Likely benign]" + + +def test_unknown_meaningful_term_fails_closed(): + with pytest.raises(ValueError, match="Unparsed clinical-classification"): + normalized_clinical_terms("new unsupported term") + + +def test_exact_audited_vocabulary_and_groups(): + assert len(APPROVED_TERMS) == 15 + assert set(APPROVED_TERMS) == { + "Benign", "Likely benign", "Benign/Likely benign", "Pathogenic", + "Likely pathogenic", "Pathogenic/Likely pathogenic", + "Uncertain significance", "Conflicting interpretations of pathogenicity", + "association", "risk factor", "protective", "drug response", "Affects", + "other", "not provided", + } + assert display_group_for_term("association") == "Association and risk context" + assert display_group_for_term("not provided") == "Details only"