fix: restore gnomAD and expression evidence layers for complete 6-layer scoring

Three bugs prevented gnomAD and expression data from contributing to scores:
1. gnomAD COLUMN_VARIANTS mapped "gene" (HGNC symbol) to gene_id instead of
   gene_symbol, causing JOIN miss with gene_universe (Ensembl IDs)
2. Expression HPA data was fetched but never merged (lf_hpa unused)
3. GTEx versioned Ensembl IDs (ENSG*.5) didn't match gene_universe

Results: gnomAD 78.5% coverage, expression 87.4%, 19946 genes with ≥4 layers.
HIGH tier refined from 44 → 18 candidates. Validation PASSED (CDH23 96.5th pctl).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit is contained in:
2026-02-16 05:25:37 +08:00
parent b63251a996
commit fe8e13c1a1
13 changed files with 18927 additions and 17136 deletions

Binary file not shown.

View File

@@ -1,12 +1,12 @@
generated_at: '2026-02-15T20:47:53.707513+00:00' generated_at: '2026-02-15T21:13:11.954116+00:00'
output_files: output_files:
- candidates.tsv - candidates.tsv
- candidates.parquet - candidates.parquet
statistics: statistics:
total_candidates: 19342 total_candidates: 21103
high_count: 44 high_count: 18
medium_count: 7268 medium_count: 9577
low_count: 12030 low_count: 11508
column_count: 22 column_count: 22
column_names: column_names:
- gene_id - gene_id

File diff suppressed because it is too large Load Diff

Binary file not shown.

Before

Width:  |  Height:  |  Size: 116 KiB

After

Width:  |  Height:  |  Size: 116 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 87 KiB

After

Width:  |  Height:  |  Size: 90 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 94 KiB

After

Width:  |  Height:  |  Size: 88 KiB

View File

@@ -1,6 +1,6 @@
{ {
"run_id": "0bcdae80-84c2-486b-809a-36040ce4821d", "run_id": "e7486ff1-f9be-403b-a68d-115fc845f4a1",
"timestamp": "2026-02-15T20:47:54.482074+00:00", "timestamp": "2026-02-15T21:13:12.288563+00:00",
"pipeline_version": "0.1.0", "pipeline_version": "0.1.0",
"parameters": { "parameters": {
"gnomad": 0.2, "gnomad": 0.2,
@@ -49,9 +49,9 @@
], ],
"validation_metrics": {}, "validation_metrics": {},
"tier_statistics": { "tier_statistics": {
"total": 19342, "total": 21103,
"high": 44, "high": 18,
"medium": 7268, "medium": 9577,
"low": 12030 "low": 11508
} }
} }

View File

@@ -1,7 +1,7 @@
# Pipeline Reproducibility Report # Pipeline Reproducibility Report
**Run ID:** `0bcdae80-84c2-486b-809a-36040ce4821d` **Run ID:** `e7486ff1-f9be-403b-a68d-115fc845f4a1`
**Timestamp:** 2026-02-15T20:47:54.482074+00:00 **Timestamp:** 2026-02-15T21:13:12.288563+00:00
**Pipeline Version:** 0.1.0 **Pipeline Version:** 0.1.0
## Parameters ## Parameters
@@ -39,7 +39,7 @@
## Tier Statistics ## Tier Statistics
- **Total Candidates:** 19342 - **Total Candidates:** 21103
- **HIGH:** 44 - **HIGH:** 18
- **MEDIUM:** 7268 - **MEDIUM:** 9577
- **LOW:** 12030 - **LOW:** 11508

View File

@@ -1360,11 +1360,15 @@ def expression_cmd(ctx, force, skip_cellxgene):
click.echo(" Skipping CellxGene (--skip-cellxgene flag)") click.echo(" Skipping CellxGene (--skip-cellxgene flag)")
try: try:
# Build gene_symbol_map for HPA merge (HPA uses gene_symbol, not gene_id)
gene_symbol_map = gene_universe.select(["gene_id", "gene_symbol"])
df = process_expression_evidence( df = process_expression_evidence(
gene_ids=gene_ids, gene_ids=gene_ids,
cache_dir=expression_dir, cache_dir=expression_dir,
force=force, force=force,
skip_cellxgene=skip_cellxgene, skip_cellxgene=skip_cellxgene,
gene_symbol_map=gene_symbol_map,
) )
click.echo(click.style( click.echo(click.style(
f" Processed {len(df)} genes", f" Processed {len(df)} genes",

View File

@@ -389,6 +389,12 @@ def fetch_gtex_expression(
lf = lf.select(select_cols).rename(rename_map) lf = lf.select(select_cols).rename(rename_map)
# Strip version suffix from Ensembl gene IDs (e.g., ENSG00000223972.5 → ENSG00000223972)
# GTEx uses versioned IDs but gene_universe uses unversioned
lf = lf.with_columns(
pl.col("gene_id").str.replace(r"\.\d+$", "").alias("gene_id")
)
# Add NULL columns for missing tissues # Add NULL columns for missing tissues
for our_key, gtex_tissue in target_tissue_cols.items(): for our_key, gtex_tissue in target_tissue_cols.items():
col_name = f"gtex_{our_key}_tpm" col_name = f"gtex_{our_key}_tpm"

View File

@@ -187,6 +187,7 @@ def process_expression_evidence(
cache_dir: Optional[Path] = None, cache_dir: Optional[Path] = None,
force: bool = False, force: bool = False,
skip_cellxgene: bool = False, skip_cellxgene: bool = False,
gene_symbol_map: Optional[pl.DataFrame] = None,
) -> pl.DataFrame: ) -> pl.DataFrame:
"""End-to-end expression evidence processing pipeline. """End-to-end expression evidence processing pipeline.
@@ -217,17 +218,23 @@ def process_expression_evidence(
gene_universe = pl.LazyFrame({"gene_id": gene_ids}) gene_universe = pl.LazyFrame({"gene_id": gene_ids})
# Merge GTEx with gene universe (left join to preserve all genes) # Merge GTEx with gene universe (left join to preserve all genes)
# GTEx has gene_id, HPA has gene_symbol - need to handle join carefully
lf_merged = gene_universe.join(lf_gtex, on="gene_id", how="left") lf_merged = gene_universe.join(lf_gtex, on="gene_id", how="left")
# For HPA, we need gene_symbol mapping # Merge HPA data via gene_symbol mapping
# We'll need to load gene universe with gene_symbol from DuckDB or pass it in # HPA returns gene_symbol as key; we need gene_symbol_map to bridge to gene_id
# For now, we'll fetch HPA separately and join on gene_symbol later if gene_symbol_map is not None:
# This requires gene_symbol in our gene_ids input or from gene universe logger.info("merging_hpa_via_symbol_map")
# lf_hpa has: gene_symbol, hpa_retina_tpm, hpa_cerebellum_tpm, ...
# DEVIATION: HPA uses gene_symbol, but we're working with gene_ids # gene_symbol_map has: gene_id, gene_symbol
# We need gene_symbol mapping. For simplicity, we'll collect HPA separately # Join HPA → symbol_map to get gene_id, then join into merged
# and merge in load.py after enriching with gene_symbol from gene universe lf_hpa_with_id = lf_hpa.join(
gene_symbol_map.select(["gene_id", "gene_symbol"]).lazy(),
on="gene_symbol",
how="inner",
).drop("gene_symbol")
lf_merged = lf_merged.join(lf_hpa_with_id, on="gene_id", how="left")
else:
logger.warning("hpa_skipped_no_symbol_map", msg="gene_symbol_map not provided; HPA data will be NULL")
# Fetch CellxGene if not skipped # Fetch CellxGene if not skipped
if not skip_cellxgene: if not skip_cellxgene:

View File

@@ -29,6 +29,19 @@ def load_to_duckdb(
""" """
logger.info("gnomad_load_start", row_count=len(df)) logger.info("gnomad_load_start", row_count=len(df))
# Enrich with Ensembl gene_id from gene_universe if missing
# gnomAD data only has gene_symbol (HGNC); we need Ensembl gene_id for scoring JOINs
if "gene_id" not in df.columns or df["gene_id"].null_count() == len(df):
logger.info("gnomad_enriching_gene_ids", msg="Mapping gene_symbol to Ensembl gene_id via gene_universe")
gene_map = store.conn.execute(
"SELECT gene_id, gene_symbol FROM gene_universe"
).pl()
if "gene_id" in df.columns:
df = df.drop("gene_id")
df = df.join(gene_map, on="gene_symbol", how="left")
matched = df.filter(pl.col("gene_id").is_not_null()).height
logger.info("gnomad_gene_id_enrichment", matched=matched, total=len(df))
# Calculate summary statistics for provenance # Calculate summary statistics for provenance
measured_count = df.filter(pl.col("quality_flag") == "measured").height measured_count = df.filter(pl.col("quality_flag") == "measured").height
incomplete_count = df.filter(pl.col("quality_flag") == "incomplete_coverage").height incomplete_count = df.filter(pl.col("quality_flag") == "incomplete_coverage").height

View File

@@ -13,8 +13,8 @@ GNOMAD_CONSTRAINT_URL = (
# v4.x uses: gene, transcript, mane_select, lof.pLI, lof.oe_ci.upper (LOEUF), mean_proportion_covered # v4.x uses: gene, transcript, mane_select, lof.pLI, lof.oe_ci.upper (LOEUF), mean_proportion_covered
# NOTE: In gnomAD data, what's called "upper" is actually the LOEUF value we want (observed/expected upper bound) # NOTE: In gnomAD data, what's called "upper" is actually the LOEUF value we want (observed/expected upper bound)
COLUMN_VARIANTS = { COLUMN_VARIANTS = {
"gene_id": ["gene", "gene_id"], "gene_id": ["gene_id"], # gnomAD doesn't provide Ensembl IDs; enriched from gene_universe in load
"gene_symbol": ["gene_symbol", "gene"], "gene_symbol": ["gene", "gene_symbol"], # gnomAD "gene" column = HGNC symbol
"transcript": ["transcript", "canonical_transcript", "mane_select"], "transcript": ["transcript", "canonical_transcript", "mane_select"],
"pli": ["pLI", "lof.pLI", "pli"], "pli": ["pLI", "lof.pLI", "pli"],
# LOEUF is the "upper" column in gnomAD (oe_lof_upper = observed/expected upper bound) # LOEUF is the "upper" column in gnomAD (oe_lof_upper = observed/expected upper bound)