usher-exploring/.planning/phases/04-scoring-integration/04-RESEARCH.md

Phase 4: Scoring & Integration - Research

Researched: 2026-02-11 Domain: Multi-Evidence Weighted Scoring Systems for Gene Prioritization Confidence: MEDIUM-HIGH

Summary

Phase 4 implements a multi-evidence weighted scoring system that integrates six heterogeneous evidence layers (genetic constraint, expression, annotation, localization, animal models, literature) into a composite gene prioritization score. This research identifies best practices for weighted scoring with missing data, positive control validation, quality control checks, and DuckDB-based integration.

The core technical challenges are: (1) handling missing data as "unknown" rather than penalizing genes (NULL preservation pattern established in prior phases), (2) configurable per-layer weights that sum to 1.0, (3) validating scoring logic using known cilia/Usher genes as positive controls before exclusion, and (4) quality control checks for distribution anomalies, outliers, and missing data rates per layer.

The existing codebase already has the necessary foundation: Pydantic ScoringWeights schema with validation, DuckDB storage of normalized 0-1 scores per evidence layer, polars for efficient DataFrame operations, and provenance tracking. The scoring phase extends these patterns to read all evidence tables, apply configurable weights, handle NULLs explicitly via SQL COALESCE with "unknown" flags, and compute composite scores with quality metrics.

Primary recommendation: Use DuckDB SQL to join all evidence layer tables, apply weighted scoring with COALESCE for NULL handling, compute composite scores and evidence breadth metrics, validate against known gene lists (CiliaCarta SCGS v2, SYSCILIA, OMIM Usher genes), and implement distribution-based QC checks (missing data rates, score percentiles, outlier detection via MAD). Store results in DuckDB with per-gene quality flags and per-layer contribution breakdowns.

Standard Stack

Core

Library	Version	Purpose	Why Standard
DuckDB	1.4+	Multi-table joins, weighted aggregation	Native polars integration, columnar aggregation optimized for scoring queries, SQL COALESCE for NULL handling
polars	1.38+	DataFrame operations, score computation	Already in codebase, 5-10x faster than pandas, lazy evaluation for large gene sets
pydantic	2.12+	Config validation for weights	Already in codebase, validates weight constraints (0-1 range, sum to 1.0)
scipy	1.14+	Statistical QC (MAD, percentiles, distributions)	Industry standard for scientific computing, robust outlier detection methods
structlog	25.5+	Structured logging for QC warnings	Already in codebase, essential for tracking missing data rates and anomalies

Supporting

Library	Version	Purpose	When to Use
numpy	2.0+	Vectorized score computation	If polars expressions insufficient for complex formulas
pandas	2.2+	Interop with scipy for QC checks	scipy functions expect pandas/numpy, not polars

Alternatives Considered

Instead of	Could Use	Tradeoff
DuckDB SQL	Pure polars joins	DuckDB optimized for multi-table analytics, SQL COALESCE cleaner than polars.when chains
scipy.stats	Custom percentile/MAD	scipy provides robust outlier detection (MAD-based), well-tested statistical functions
Weighted average	Manual score loops	DuckDB/polars vectorized operations 100x faster than Python loops

Installation:

# scipy is the only new dependency (others already in codebase)
pip install 'scipy>=1.14'

Architecture Patterns

Recommended Project Structure

src/
├── scoring/
│   ├── __init__.py
│   ├── integration.py      # Multi-evidence weighted scoring
│   ├── known_genes.py      # CiliaCarta, SYSCILIA, OMIM gene lists
│   ├── quality_control.py  # QC checks (missing data, distributions, outliers)
│   └── validation.py       # Positive control validation (known gene ranking)
├── cli/
│   └── score.py            # CLI command for scoring phase
└── config/
    └── schema.py           # ScoringWeights (already exists)

Pattern 1: Multi-Table Join with NULL Preservation

What: Join all evidence layer tables on gene_id, preserving NULLs to distinguish "missing data" from "zero score" When to use: Integrating 6+ evidence tables with heterogeneous coverage Example:

# Source: DuckDB NULL handling best practices
import duckdb

def join_evidence_layers(con: duckdb.DuckDBPyConnection) -> duckdb.DuckDBPyRelation:
    """Join all evidence layers preserving NULL for missing data.

    Returns relation with columns:
    - gene_id
    - gene_symbol
    - gnomad_score (NULL if no constraint data)
    - expression_score (NULL if no expression data)
    - annotation_score (NULL if no annotation data)
    - localization_score (NULL if no localization data)
    - animal_model_score (NULL if no phenotype data)
    - literature_score (NULL if no literature data)
    - evidence_count (number of non-NULL scores)
    """
    query = """
    SELECT
        g.gene_id,
        g.gene_symbol,
        c.loeuf_normalized AS gnomad_score,
        e.composite_expression_score AS expression_score,
        a.composite_annotation_score AS annotation_score,
        l.localization_score,
        am.phenotype_score AS animal_model_score,
        lit.literature_score,
        -- Count non-NULL evidence layers
        (
            CASE WHEN c.loeuf_normalized IS NOT NULL THEN 1 ELSE 0 END +
            CASE WHEN e.composite_expression_score IS NOT NULL THEN 1 ELSE 0 END +
            CASE WHEN a.composite_annotation_score IS NOT NULL THEN 1 ELSE 0 END +
            CASE WHEN l.localization_score IS NOT NULL THEN 1 ELSE 0 END +
            CASE WHEN am.phenotype_score IS NOT NULL THEN 1 ELSE 0 END +
            CASE WHEN lit.literature_score IS NOT NULL THEN 1 ELSE 0 END
        ) AS evidence_count
    FROM genes g
    LEFT JOIN gnomad_constraint c ON g.gene_id = c.gene_id
    LEFT JOIN expression_evidence e ON g.gene_id = e.gene_id
    LEFT JOIN annotation_evidence a ON g.gene_id = a.gene_id
    LEFT JOIN localization_evidence l ON g.gene_id = l.gene_id
    LEFT JOIN animal_model_evidence am ON g.gene_id = am.gene_id
    LEFT JOIN literature_evidence lit ON g.gene_id = lit.gene_id
    """
    return con.sql(query)

Pattern 2: Weighted Scoring with COALESCE and Unknown Flags

What: Compute weighted composite score treating NULLs as "unknown" rather than zero penalty When to use: Scoring with missing data that should not penalize genes Example:

# Source: DuckDB COALESCE pattern + weighted scoring best practices
def compute_weighted_scores(
    con: duckdb.DuckDBPyConnection,
    weights: dict[str, float]
) -> duckdb.DuckDBPyRelation:
    """Compute weighted composite scores handling missing data explicitly.

    Scoring strategy:
    1. Composite score = weighted average of AVAILABLE evidence only
    2. NULL evidence layers don't contribute (neither positive nor negative)
    3. Quality flag marks genes with <3 evidence layers as "sparse_evidence"

    Args:
        weights: dict with keys matching evidence layers (gnomad, expression, etc.)

    Returns:
        Relation with composite_score, available_weight, quality_flag columns
    """
    query = f"""
    WITH evidence AS (
        -- Join all layers (from Pattern 1)
        {join_evidence_layers(con).sql_query()}
    ),
    weighted AS (
        SELECT
            gene_id,
            gene_symbol,
            evidence_count,
            -- Calculate sum of weights for AVAILABLE evidence
            (
                CASE WHEN gnomad_score IS NOT NULL THEN {weights['gnomad']} ELSE 0 END +
                CASE WHEN expression_score IS NOT NULL THEN {weights['expression']} ELSE 0 END +
                CASE WHEN annotation_score IS NOT NULL THEN {weights['annotation']} ELSE 0 END +
                CASE WHEN localization_score IS NOT NULL THEN {weights['localization']} ELSE 0 END +
                CASE WHEN animal_model_score IS NOT NULL THEN {weights['animal_model']} ELSE 0 END +
                CASE WHEN literature_score IS NOT NULL THEN {weights['literature']} ELSE 0 END
            ) AS available_weight,
            -- Calculate weighted sum (only non-NULL layers)
            (
                COALESCE(gnomad_score * {weights['gnomad']}, 0) +
                COALESCE(expression_score * {weights['expression']}, 0) +
                COALESCE(annotation_score * {weights['annotation']}, 0) +
                COALESCE(localization_score * {weights['localization']}, 0) +
                COALESCE(animal_model_score * {weights['animal_model']}, 0) +
                COALESCE(literature_score * {weights['literature']}, 0)
            ) AS weighted_sum,
            -- Individual layer contributions for explainability
            gnomad_score,
            expression_score,
            annotation_score,
            localization_score,
            animal_model_score,
            literature_score
        FROM evidence
    )
    SELECT
        gene_id,
        gene_symbol,
        evidence_count,
        available_weight,
        -- Composite score = weighted sum / available weight
        -- NULL if no evidence layers available
        CASE
            WHEN available_weight > 0 THEN weighted_sum / available_weight
            ELSE NULL
        END AS composite_score,
        -- Quality flag based on evidence breadth
        CASE
            WHEN evidence_count >= 4 THEN 'sufficient_evidence'
            WHEN evidence_count >= 2 THEN 'moderate_evidence'
            WHEN evidence_count >= 1 THEN 'sparse_evidence'
            ELSE 'no_evidence'
        END AS quality_flag,
        -- Individual scores for explainability
        gnomad_score,
        expression_score,
        annotation_score,
        localization_score,
        animal_model_score,
        literature_score
    FROM weighted
    """
    return con.sql(query)

Pattern 3: Known Gene Lists for Positive Control Validation

What: Compile known cilia/Usher genes from authoritative sources, validate scoring by checking if they rank highly When to use: Before finalizing scoring weights, to validate scoring logic works Example:

# Source: Gene prioritization validation best practices
import httpx
import polars as pl
from pathlib import Path

def fetch_known_cilia_genes() -> pl.DataFrame:
    """Fetch known ciliary genes from CiliaCarta and SYSCILIA gold standard.

    Returns DataFrame with columns:
    - gene_symbol
    - source (ciliacarta_scgs_v2, syscilia, omim_usher)
    - confidence (HIGH for SCGS, MEDIUM for predicted, HIGH for OMIM)
    """
    known_genes = []

    # SYSCILIA Gold Standard v2 (SCGSv2) - 686 high-confidence ciliary genes
    # Source: https://www.molbiolcell.org/doi/10.1091/mbc.E21-05-0226
    # Note: Requires manual download or scraping; no direct API
    # Placeholder - implement actual fetch during development
    scgs_genes = fetch_scgs_v2()  # Returns list of gene symbols
    for gene in scgs_genes:
        known_genes.append({
            "gene_symbol": gene,
            "source": "syscilia_scgs_v2",
            "confidence": "HIGH"
        })

    # OMIM Usher Syndrome genes
    # Source: OMIM entries for Usher syndrome (USH1A-USH3A)
    # Known genes: MYO7A, USH1C, CDH23, PCDH15, USH1G, CIB2, USH2A, ADGRV1, WHRN, CLRN1
    usher_genes = [
        "MYO7A", "USH1C", "CDH23", "PCDH15", "USH1G", "CIB2",
        "USH2A", "ADGRV1", "WHRN", "CLRN1"
    ]
    for gene in usher_genes:
        known_genes.append({
            "gene_symbol": gene,
            "source": "omim_usher",
            "confidence": "HIGH"
        })

    return pl.DataFrame(known_genes)

def validate_known_gene_ranking(
    con: duckdb.DuckDBPyConnection,
    known_genes: pl.DataFrame,
    percentile_threshold: float = 0.75
) -> dict:
    """Validate scoring by checking if known genes rank in top percentile.

    Returns validation metrics:
    - median_rank_percentile: median percentile rank of known genes (should be >0.75)
    - top_quartile_count: number of known genes in top 25%
    - total_known_genes: total known genes in dataset
    """
    # Join scored genes with known gene list
    query = """
    SELECT
        s.gene_symbol,
        s.composite_score,
        k.source AS known_source,
        k.confidence AS known_confidence,
        -- Calculate percentile rank (0 = lowest score, 1 = highest score)
        PERCENT_RANK() OVER (ORDER BY s.composite_score) AS percentile_rank
    FROM scored_genes s
    INNER JOIN known_genes k ON s.gene_symbol = k.gene_symbol
    WHERE s.composite_score IS NOT NULL
    """
    result = con.execute(query).pl()

    median_percentile = result["percentile_rank"].median()
    top_quartile = (result["percentile_rank"] >= 0.75).sum()
    total = len(result)

    return {
        "median_rank_percentile": median_percentile,
        "top_quartile_count": top_quartile,
        "total_known_genes": total,
        "validation_passed": median_percentile >= percentile_threshold
    }

Pattern 4: Distribution-Based Quality Control

What: Detect anomalies in score distributions, outliers, and missing data rates per layer When to use: After scoring, before outputting results Example:

# Source: Quality control for genomic scoring systems
import scipy.stats as stats
import numpy as np

def compute_qc_metrics(con: duckdb.DuckDBPyConnection) -> dict:
    """Compute QC metrics for score distributions and missing data.

    Returns dict with:
    - missing_data_rates: dict per layer
    - score_distributions: dict per layer (mean, median, std, skew)
    - outliers: dict per layer (count, gene_symbols)
    - composite_score_stats: overall distribution stats
    """
    qc_metrics = {}

    # Missing data rates per layer
    missing_query = """
    SELECT
        COUNT(*) AS total_genes,
        SUM(CASE WHEN gnomad_score IS NULL THEN 1 ELSE 0 END) AS gnomad_missing,
        SUM(CASE WHEN expression_score IS NULL THEN 1 ELSE 0 END) AS expression_missing,
        SUM(CASE WHEN annotation_score IS NULL THEN 1 ELSE 0 END) AS annotation_missing,
        SUM(CASE WHEN localization_score IS NULL THEN 1 ELSE 0 END) AS localization_missing,
        SUM(CASE WHEN animal_model_score IS NULL THEN 1 ELSE 0 END) AS animal_model_missing,
        SUM(CASE WHEN literature_score IS NULL THEN 1 ELSE 0 END) AS literature_missing
    FROM scored_genes
    """
    missing_result = con.execute(missing_query).fetchone()
    total = missing_result[0]

    qc_metrics["missing_data_rates"] = {
        "gnomad": missing_result[1] / total,
        "expression": missing_result[2] / total,
        "annotation": missing_result[3] / total,
        "localization": missing_result[4] / total,
        "animal_model": missing_result[5] / total,
        "literature": missing_result[6] / total,
    }

    # Score distributions per layer (using polars for statistical functions)
    layers = ["gnomad", "expression", "annotation", "localization", "animal_model", "literature"]
    qc_metrics["score_distributions"] = {}

    for layer in layers:
        layer_scores = con.execute(
            f"SELECT {layer}_score FROM scored_genes WHERE {layer}_score IS NOT NULL"
        ).pl()[f"{layer}_score"].to_numpy()

        if len(layer_scores) > 0:
            qc_metrics["score_distributions"][layer] = {
                "mean": float(np.mean(layer_scores)),
                "median": float(np.median(layer_scores)),
                "std": float(np.std(layer_scores)),
                "skew": float(stats.skew(layer_scores)),
                "kurtosis": float(stats.kurtosis(layer_scores)),
            }

    # Outlier detection using MAD (Median Absolute Deviation)
    # Source: https://academic.oup.com/nargab/article/3/1/lqab005/6155871
    qc_metrics["outliers"] = {}

    for layer in layers:
        layer_df = con.execute(
            f"SELECT gene_symbol, {layer}_score FROM scored_genes WHERE {layer}_score IS NOT NULL"
        ).pl()

        scores = layer_df[f"{layer}_score"].to_numpy()
        median = np.median(scores)
        mad = np.median(np.abs(scores - median))

        # Outliers: >3 MAD from median
        threshold = 3
        outlier_mask = np.abs(scores - median) > threshold * mad
        outlier_genes = layer_df.filter(pl.Series(outlier_mask))["gene_symbol"].to_list()

        qc_metrics["outliers"][layer] = {
            "count": len(outlier_genes),
            "genes": outlier_genes[:10] if len(outlier_genes) > 10 else outlier_genes  # Limit to 10 for logging
        }

    # Composite score distribution
    composite_scores = con.execute(
        "SELECT composite_score FROM scored_genes WHERE composite_score IS NOT NULL"
    ).pl()["composite_score"].to_numpy()

    qc_metrics["composite_score_stats"] = {
        "mean": float(np.mean(composite_scores)),
        "median": float(np.median(composite_scores)),
        "std": float(np.std(composite_scores)),
        "percentiles": {
            "p10": float(np.percentile(composite_scores, 10)),
            "p25": float(np.percentile(composite_scores, 25)),
            "p50": float(np.percentile(composite_scores, 50)),
            "p75": float(np.percentile(composite_scores, 75)),
            "p90": float(np.percentile(composite_scores, 90)),
        }
    }

    return qc_metrics

Pattern 5: Scoring Results Persistence with Explainability

What: Store composite scores with per-layer contributions and quality flags in DuckDB When to use: Final step of scoring phase Example:

# Source: DuckDB analytics table best practices
def persist_scoring_results(
    con: duckdb.DuckDBPyConnection,
    weights: dict[str, float],
    config_hash: str
) -> None:
    """Persist scoring results to DuckDB with full explainability.

    Creates table: scored_genes with columns:
    - gene_id, gene_symbol
    - composite_score (weighted average of available evidence)
    - evidence_count (number of non-NULL layers)
    - available_weight (sum of weights for non-NULL layers)
    - quality_flag (sufficient/moderate/sparse/no evidence)
    - Individual layer scores (for explainability)
    - Individual layer contributions (score * weight, for debugging)
    - scoring_metadata (weights used, config hash, timestamp)
    """
    query = f"""
    CREATE OR REPLACE TABLE scored_genes AS
    WITH weighted_scores AS (
        -- Pattern 2 query here
        {compute_weighted_scores(con, weights).sql_query()}
    )
    SELECT
        gene_id,
        gene_symbol,
        composite_score,
        evidence_count,
        available_weight,
        quality_flag,
        -- Individual layer scores
        gnomad_score,
        expression_score,
        annotation_score,
        localization_score,
        animal_model_score,
        literature_score,
        -- Individual layer contributions (for explainability)
        gnomad_score * {weights['gnomad']} AS gnomad_contribution,
        expression_score * {weights['expression']} AS expression_contribution,
        annotation_score * {weights['annotation']} AS annotation_contribution,
        localization_score * {weights['localization']} AS localization_contribution,
        animal_model_score * {weights['animal_model']} AS animal_model_contribution,
        literature_score * {weights['literature']} AS literature_contribution,
        -- Metadata
        '{config_hash}' AS config_hash,
        CURRENT_TIMESTAMP AS scored_at
    FROM weighted_scores
    ORDER BY composite_score DESC NULLS LAST
    """
    con.execute(query)

Anti-Patterns to Avoid

• Don't replace NULL with zero in weighted averaging: This penalizes genes for missing data rather than treating it as "unknown"
• Don't normalize composite scores across all genes: Keep scores on 0-1 scale based on evidence layer normalization; percentile ranks can be added separately
• Don't validate scoring after excluding known genes: Validation must happen BEFORE exclusion to confirm known genes rank highly
• Don't ignore QC warnings: High missing data rates or distribution anomalies indicate upstream evidence layer failures
• Don't hardcode weights: Use Pydantic config to make weights configurable and validate they sum to 1.0

Don't Hand-Roll

Problem	Don't Build	Use Instead	Why
Outlier detection	Manual z-score thresholds	scipy MAD-based detection	MAD is robust to non-normal distributions common in genomic data; z-scores fail with skewed distributions
Percentile ranking	Manual sorting and indexing	DuckDB PERCENT_RANK() window function	Vectorized, handles ties correctly, optimized for large datasets
Statistical distributions	Custom mean/std/skew calculations	scipy.stats module	Handles edge cases (empty arrays, NaNs), numerically stable algorithms
Weighted averages with NULLs	Python loops over genes	DuckDB COALESCE in SQL query	100-1000x faster, vectorized operations, cleaner NULL semantics
Known gene list fetching	Manual web scraping	Structured data downloads + provenance	Scraping is brittle; use official TSV/CSV downloads with version tracking

Key insight: Scoring with missing data has subtle edge cases (divide by zero when all layers NULL, weight normalization, NULL vs zero distinction). DuckDB SQL and scipy provide battle-tested solutions that handle these correctly.

Common Pitfalls

Pitfall 1: Replacing NULL with Zero Penalizes Poorly-Studied Genes

What goes wrong: Using COALESCE(score, 0) in weighted average treats "no data" as "score of zero", systematically lowering scores for genes with missing evidence. This creates bias against novel/poorly-studied genes—exactly the candidates the pipeline aims to discover. Why it happens: Standard weighted average formulas assume all values are present; developers default to zero-filling for simplicity. How to avoid: Use weighted average of AVAILABLE evidence only: weighted_sum / available_weight where available_weight excludes NULL layers. Track evidence_count and flag genes with sparse evidence rather than penalizing them. Warning signs: Known well-studied genes rank higher than expected; genes with sparse annotation cluster at low scores; correlation between composite score and total publication count.

Pitfall 2: Weights Don't Sum to 1.0 Due to Configuration Errors

What goes wrong: User modifies weight config file but weights sum to 0.95 or 1.15 instead of 1.0. Composite scores become non-comparable across configurations. Results appear valid but are systematically biased. Why it happens: No validation on weight sum; YAML editing errors; copy-paste mistakes. How to avoid: Add Pydantic validator to ScoringWeights that asserts sum equals 1.0 (within floating point tolerance). Reject invalid configs at pipeline startup. Warning signs: Composite scores outside expected 0-1 range; scores don't match manual calculations; weight changes have disproportionate impact.

Pitfall 3: Known Gene Validation Happens After Exclusion

What goes wrong: Known cilia/Usher genes are excluded from results, then validation tries to check their ranking. Validation fails because genes aren't in scored dataset. Developer concludes scoring is broken when it's actually working correctly. Why it happens: Misunderstanding of validation purpose—it's to validate scoring logic, not final results. How to avoid: Run validation on FULL scored gene set BEFORE applying any filters or exclusions. Log validation metrics (median percentile, top quartile count) to confirm known genes rank highly. Then exclude known genes for final output. Warning signs: Validation reports zero known genes found; known genes absent from top-ranked results; validation fails every time.

Pitfall 4: Missing Data Rates Not Monitored, Upstream Failures Silent

What goes wrong: Expression evidence layer fails to retrieve data for 90% of genes (API error, identifier mismatch, etc.). Scoring runs successfully but expression weight is essentially ignored. Results are dominated by other layers without user awareness. Why it happens: No QC checks on missing data rates; scoring accepts any evidence table structure. How to avoid: Compute missing data rate per layer (% of genes with NULL). Log warnings if rate exceeds thresholds (>50% = warning, >80% = error). Include missing data rates in provenance/QC report. Warning signs: Composite scores barely change when expression weight is modified; evidence_count clusters at low values; layer-specific scores are NULL for most genes.

Pitfall 5: Score Distributions Not Inspected, Anomalies Undetected

What goes wrong: All genes in localization layer score exactly 1.0 (normalization bug), or annotation scores cluster at 0.0 and 1.0 with nothing in between (bad binning). Composite scores are distorted but no one notices until manual inspection of top candidates. Why it happens: No automated distribution checks; developers assume normalization worked correctly. How to avoid: Compute distribution stats per layer (mean, median, std, skew). Flag anomalies: std < 0.01 (no variation), skew > 2.0 (extreme skew), > 50% genes at single value. Log distribution stats to QC report. Warning signs: All top-ranked genes have same composite score; weights have no effect on results; manual inspection reveals nonsensical patterns.

Pitfall 6: Outlier Genes Not Flagged, Distort Results

What goes wrong: Single gene has animal_model_score = 50.0 due to upstream normalization bug (should be 0-1). Weighted average is dominated by this outlier. Gene ranks #1 despite weak evidence in other layers. Why it happens: No outlier detection; assumption that all evidence layers are correctly normalized. How to avoid: Use MAD-based outlier detection per layer (>3 MAD from median = outlier). Flag outlier genes in quality_flag column. Log outlier counts to QC report. Consider capping scores at 1.0 as sanity check. Warning signs: Top-ranked genes have suspiciously high composite scores (>1.0); single layer dominates composite score; manual inspection reveals impossible values.

Pitfall 7: Evidence Count Used as Scoring Factor Instead of Quality Flag

What goes wrong: Developer modifies scoring to boost genes with high evidence_count (e.g., composite_score * log(evidence_count)). This creates bias toward well-studied genes that have data in all layers, defeating the purpose of NULL preservation. Why it happens: Desire to reward "comprehensive evidence"; misunderstanding that evidence breadth ≠ evidence quality. How to avoid: Use evidence_count ONLY as quality flag for downstream filtering (e.g., "high confidence" tier requires ≥4 layers). Never multiply composite score by evidence count. Trust weighted average of available evidence. Warning signs: Known under-studied genes rank lower than expected; correlation between composite score and total publications; genes with sparse evidence systematically deprioritized.

Code Examples

Verified patterns from official sources:

DuckDB Multi-Table Join with NULL Preservation

-- Source: DuckDB NULL handling documentation
-- Join all evidence layers, preserving NULL for missing data
CREATE OR REPLACE VIEW evidence_integrated AS
SELECT
    g.gene_id,
    g.gene_symbol,
    -- Evidence layer scores (NULL if layer has no data for gene)
    c.loeuf_normalized AS gnomad_score,
    e.composite_expression_score AS expression_score,
    a.composite_annotation_score AS annotation_score,
    l.localization_score,
    am.phenotype_score AS animal_model_score,
    lit.literature_score,
    -- Count non-NULL layers
    (
        (CASE WHEN c.loeuf_normalized IS NOT NULL THEN 1 ELSE 0 END) +
        (CASE WHEN e.composite_expression_score IS NOT NULL THEN 1 ELSE 0 END) +
        (CASE WHEN a.composite_annotation_score IS NOT NULL THEN 1 ELSE 0 END) +
        (CASE WHEN l.localization_score IS NOT NULL THEN 1 ELSE 0 END) +
        (CASE WHEN am.phenotype_score IS NOT NULL THEN 1 ELSE 0 END) +
        (CASE WHEN lit.literature_score IS NOT NULL THEN 1 ELSE 0 END)
    ) AS evidence_count
FROM genes g
LEFT JOIN gnomad_constraint c ON g.gene_id = c.gene_id
LEFT JOIN expression_evidence e ON g.gene_id = e.gene_id
LEFT JOIN annotation_evidence a ON g.gene_id = a.gene_id
LEFT JOIN localization_evidence l ON g.gene_id = l.gene_id
LEFT JOIN animal_model_evidence am ON g.gene_id = am.gene_id
LEFT JOIN literature_evidence lit ON g.gene_id = lit.gene_id;

Weighted Scoring with Configurable Weights

# Source: Pydantic config + DuckDB weighted aggregation
from pydantic import BaseModel, field_validator
import duckdb

class ScoringWeights(BaseModel):
    """Weights for evidence layers (must sum to 1.0)."""
    gnomad: float = 0.20
    expression: float = 0.20
    annotation: float = 0.15
    localization: float = 0.15
    animal_model: float = 0.15
    literature: float = 0.15

    @field_validator("*")
    @classmethod
    def weights_in_range(cls, v: float) -> float:
        if not 0.0 <= v <= 1.0:
            raise ValueError(f"Weight must be in [0.0, 1.0], got {v}")
        return v

    def validate_sum(self) -> None:
        """Validate weights sum to 1.0."""
        total = sum([
            self.gnomad, self.expression, self.annotation,
            self.localization, self.animal_model, self.literature
        ])
        if not abs(total - 1.0) < 1e-6:
            raise ValueError(f"Weights must sum to 1.0, got {total}")

def compute_composite_scores(con: duckdb.DuckDBPyConnection, weights: ScoringWeights) -> None:
    """Compute weighted composite scores with NULL handling."""
    weights.validate_sum()  # Assert weights sum to 1.0

    query = f"""
    CREATE OR REPLACE TABLE scored_genes AS
    SELECT
        gene_id,
        gene_symbol,
        evidence_count,
        -- Calculate available weight (sum of weights for non-NULL layers)
        (
            (CASE WHEN gnomad_score IS NOT NULL THEN {weights.gnomad} ELSE 0 END) +
            (CASE WHEN expression_score IS NOT NULL THEN {weights.expression} ELSE 0 END) +
            (CASE WHEN annotation_score IS NOT NULL THEN {weights.annotation} ELSE 0 END) +
            (CASE WHEN localization_score IS NOT NULL THEN {weights.localization} ELSE 0 END) +
            (CASE WHEN animal_model_score IS NOT NULL THEN {weights.animal_model} ELSE 0 END) +
            (CASE WHEN literature_score IS NOT NULL THEN {weights.literature} ELSE 0 END)
        ) AS available_weight,
        -- Weighted sum (COALESCE NULL to 0 for addition, not for penalization)
        (
            COALESCE(gnomad_score * {weights.gnomad}, 0) +
            COALESCE(expression_score * {weights.expression}, 0) +
            COALESCE(annotation_score * {weights.annotation}, 0) +
            COALESCE(localization_score * {weights.localization}, 0) +
            COALESCE(animal_model_score * {weights.animal_model}, 0) +
            COALESCE(literature_score * {weights.literature}, 0)
        ) AS weighted_sum,
        -- Composite score = weighted average of AVAILABLE evidence
        CASE
            WHEN available_weight > 0 THEN weighted_sum / available_weight
            ELSE NULL
        END AS composite_score,
        -- Quality flag
        CASE
            WHEN evidence_count >= 4 THEN 'sufficient_evidence'
            WHEN evidence_count >= 2 THEN 'moderate_evidence'
            WHEN evidence_count >= 1 THEN 'sparse_evidence'
            ELSE 'no_evidence'
        END AS quality_flag,
        -- Individual scores for explainability
        gnomad_score,
        expression_score,
        annotation_score,
        localization_score,
        animal_model_score,
        literature_score
    FROM evidence_integrated
    ORDER BY composite_score DESC NULLS LAST
    """
    con.execute(query)

Known Gene Validation (Positive Controls)

# Source: Gene prioritization validation best practices
import polars as pl
import structlog

logger = structlog.get_logger()

def load_known_genes() -> pl.DataFrame:
    """Load known cilia/Usher genes for validation.

    Returns DataFrame with gene_symbol, source, confidence columns.
    """
    # OMIM Usher syndrome genes (HIGH confidence)
    usher_genes = [
        "MYO7A", "USH1C", "CDH23", "PCDH15", "USH1G", "CIB2",
        "USH2A", "ADGRV1", "WHRN", "CLRN1"
    ]

    # SYSCILIA Gold Standard v2 genes (HIGH confidence)
    # Note: Full list requires downloading from publication supplementary data
    # Placeholder for now - replace with actual fetch during implementation
    scgs_v2_genes = [
        "IFT88", "IFT140", "BBS1", "BBS2", "BBS4", "BBS5",
        "RPGRIP1L", "CEP290", "ARL13B", "INPP5E"
        # ... full SCGS v2 list has 686 genes
    ]

    known_genes = []
    for gene in usher_genes:
        known_genes.append({"gene_symbol": gene, "source": "omim_usher", "confidence": "HIGH"})
    for gene in scgs_v2_genes:
        known_genes.append({"gene_symbol": gene, "source": "syscilia_scgs_v2", "confidence": "HIGH"})

    return pl.DataFrame(known_genes)

def validate_known_gene_ranking(con: duckdb.DuckDBPyConnection) -> dict:
    """Validate scoring by checking known gene rankings.

    Returns validation metrics dict.
    """
    known_genes = load_known_genes()

    # Create temp table for known genes
    con.execute("CREATE TEMP TABLE known_genes AS SELECT * FROM known_genes_df")

    # Query percentile ranks for known genes
    query = """
    SELECT
        s.gene_symbol,
        s.composite_score,
        k.source,
        PERCENT_RANK() OVER (ORDER BY s.composite_score) AS percentile_rank
    FROM scored_genes s
    INNER JOIN known_genes k ON s.gene_symbol = k.gene_symbol
    WHERE s.composite_score IS NOT NULL
    """
    result = con.execute(query).pl()

    if len(result) == 0:
        logger.error("validation_failed", reason="No known genes found in scored dataset")
        return {"validation_passed": False, "reason": "no_known_genes_found"}

    median_percentile = result["percentile_rank"].median()
    top_quartile_count = (result["percentile_rank"] >= 0.75).sum()
    total = len(result)

    metrics = {
        "median_percentile": float(median_percentile),
        "top_quartile_count": int(top_quartile_count),
        "total_known_genes": total,
        "top_quartile_fraction": top_quartile_count / total,
        "validation_passed": median_percentile >= 0.75  # Known genes should rank in top quartile
    }

    if metrics["validation_passed"]:
        logger.info("validation_passed", **metrics)
    else:
        logger.warning("validation_failed", **metrics)

    return metrics

Quality Control Checks

# Source: Genomic data QC best practices
import numpy as np
import scipy.stats as stats
import structlog

logger = structlog.get_logger()

def run_qc_checks(con: duckdb.DuckDBPyConnection) -> dict:
    """Run quality control checks on scoring results.

    Returns QC metrics dict with warnings/errors.
    """
    qc_metrics = {"warnings": [], "errors": []}

    # 1. Missing data rates per layer
    missing_query = """
    SELECT
        COUNT(*) AS total,
        AVG(CASE WHEN gnomad_score IS NULL THEN 1.0 ELSE 0.0 END) AS gnomad_missing,
        AVG(CASE WHEN expression_score IS NULL THEN 1.0 ELSE 0.0 END) AS expression_missing,
        AVG(CASE WHEN annotation_score IS NULL THEN 1.0 ELSE 0.0 END) AS annotation_missing,
        AVG(CASE WHEN localization_score IS NULL THEN 1.0 ELSE 0.0 END) AS localization_missing,
        AVG(CASE WHEN animal_model_score IS NULL THEN 1.0 ELSE 0.0 END) AS animal_model_missing,
        AVG(CASE WHEN literature_score IS NULL THEN 1.0 ELSE 0.0 END) AS literature_missing
    FROM scored_genes
    """
    missing = con.execute(missing_query).fetchone()

    qc_metrics["missing_data_rates"] = {
        "gnomad": float(missing[1]),
        "expression": float(missing[2]),
        "annotation": float(missing[3]),
        "localization": float(missing[4]),
        "animal_model": float(missing[5]),
        "literature": float(missing[6]),
    }

    # Check for high missing rates
    for layer, rate in qc_metrics["missing_data_rates"].items():
        if rate > 0.8:
            qc_metrics["errors"].append(f"{layer}_missing_rate_critical: {rate:.1%}")
            logger.error("missing_data_critical", layer=layer, rate=rate)
        elif rate > 0.5:
            qc_metrics["warnings"].append(f"{layer}_missing_rate_high: {rate:.1%}")
            logger.warning("missing_data_high", layer=layer, rate=rate)

    # 2. Score distributions per layer
    layers = ["gnomad", "expression", "annotation", "localization", "animal_model", "literature"]
    qc_metrics["distributions"] = {}

    for layer in layers:
        scores = con.execute(
            f"SELECT {layer}_score FROM scored_genes WHERE {layer}_score IS NOT NULL"
        ).pl()[f"{layer}_score"].to_numpy()

        if len(scores) < 10:
            qc_metrics["errors"].append(f"{layer}_insufficient_data: only {len(scores)} genes")
            continue

        dist_stats = {
            "mean": float(np.mean(scores)),
            "median": float(np.median(scores)),
            "std": float(np.std(scores)),
            "min": float(np.min(scores)),
            "max": float(np.max(scores)),
        }
        qc_metrics["distributions"][layer] = dist_stats

        # Check for anomalies
        if dist_stats["std"] < 0.01:
            qc_metrics["warnings"].append(f"{layer}_no_variation: std={dist_stats['std']:.4f}")
        if dist_stats["min"] < 0.0 or dist_stats["max"] > 1.0:
            qc_metrics["errors"].append(f"{layer}_out_of_range: [{dist_stats['min']:.2f}, {dist_stats['max']:.2f}]")

    # 3. Outlier detection (MAD-based)
    qc_metrics["outliers"] = {}

    for layer in layers:
        scores_df = con.execute(
            f"SELECT gene_symbol, {layer}_score FROM scored_genes WHERE {layer}_score IS NOT NULL"
        ).pl()

        if len(scores_df) < 10:
            continue

        scores = scores_df[f"{layer}_score"].to_numpy()
        median = np.median(scores)
        mad = np.median(np.abs(scores - median))

        if mad == 0:
            continue  # No variation, skip outlier detection

        # Outliers: >3 MAD from median
        outlier_mask = np.abs(scores - median) > 3 * mad
        outlier_count = np.sum(outlier_mask)

        if outlier_count > 0:
            outlier_genes = scores_df.filter(pl.Series(outlier_mask))["gene_symbol"].to_list()[:5]
            qc_metrics["outliers"][layer] = {
                "count": int(outlier_count),
                "example_genes": outlier_genes
            }
            logger.warning("outliers_detected", layer=layer, count=outlier_count, examples=outlier_genes)

    # 4. Composite score distribution
    composite = con.execute(
        "SELECT composite_score FROM scored_genes WHERE composite_score IS NOT NULL"
    ).pl()["composite_score"].to_numpy()

    qc_metrics["composite_stats"] = {
        "mean": float(np.mean(composite)),
        "median": float(np.median(composite)),
        "std": float(np.std(composite)),
        "percentiles": {
            "p10": float(np.percentile(composite, 10)),
            "p25": float(np.percentile(composite, 25)),
            "p50": float(np.percentile(composite, 50)),
            "p75": float(np.percentile(composite, 75)),
            "p90": float(np.percentile(composite, 90)),
        }
    }

    return qc_metrics

State of the Art

Old Approach	Current Approach	When Changed	Impact
Impute missing data with zeros	NULL preservation with explicit flags	2020s genomics best practices	Avoids bias against poorly-studied genes
Equal weights for all layers	Configurable validated weights	Always preferred	Allows domain expertise to guide prioritization
Validation after filtering	Positive controls validated before exclusion	Gene prioritization methodology	Ensures scoring logic works correctly
Manual QC inspection	Automated distribution checks + MAD outliers	2020s data pipelines	Catches anomalies before results published
Z-score outlier detection	MAD-based robust outliers	2020s genomic QC	Handles non-normal distributions correctly
Python loops for scoring	DuckDB/polars vectorized operations	2020s performance	100-1000x faster, handles large datasets

Deprecated/outdated:

• Imputation for genomic missing data: NULL preservation is now standard in bioinformatics
• Z-score outlier detection: MAD is robust to non-normal distributions, preferred for genomic data
• Hardcoded scoring weights: Configurable weights with validation are standard practice
• Post-hoc validation: Positive control validation should happen during development, not after publication

Open Questions

1. Exact SYSCILIA Gold Standard v2 gene list source

   • What we know: SCGS v2 published in Molecular Biology of the Cell (2021) with 686 genes
   • What's unclear: Direct download URL for full gene list (supplementary data access)
   • Recommendation: Download supplementary data from publication DOI 10.1091/mbc.E21-05-0226; parse Excel/TSV; store in pipeline data directory with provenance

2. Weight configuration optimization strategy

   • What we know: Weights are configurable, default to equal (0.15-0.20 per layer)
   • What's unclear: Whether to optimize weights based on known gene validation or use domain expertise
   • Recommendation: Start with default weights, validate against known genes, document weight choice rationale; consider v2 feature for sensitivity analysis across weight configs

3. Evidence count threshold for quality tiers

   • What we know: Genes with sparse evidence should be flagged, not penalized
   • What's unclear: Optimal threshold for "sufficient" vs "moderate" vs "sparse" evidence (3+ layers? 4+?)
   • Recommendation: Use 4+ layers = sufficient, 2-3 = moderate, 1 = sparse based on 6 total layers (2/3 majority); document in code

4. Outlier handling strategy

   • What we know: MAD-based detection identifies outliers (>3 MAD from median)
   • What's unclear: Should outliers be capped at 1.0, flagged but kept, or excluded?
   • Recommendation: Flag outliers in QC report but keep in dataset; cap composite_score at 1.0 as sanity check; investigate upstream normalization bugs

5. Missing data rate thresholds for errors vs warnings

   • What we know: High missing rates indicate upstream failures
   • What's unclear: Exact thresholds for QC warnings (50%? 60%?) and errors (80%? 90%?)
   • Recommendation: Warn at 50%, error at 80% missing rate; document in QC module; adjust based on real data during implementation

Sources

Primary (HIGH confidence)

• DuckDB NULL Handling - Official NULL semantics and COALESCE
• DuckDB Window Functions - PERCENT_RANK() for validation
• DuckDB Polars Integration - Zero-copy Arrow conversion
• scipy.stats Documentation - Statistical functions for QC
• SYSCILIA Gold Standard v2 - 686 high-confidence ciliary genes
• CiliaCarta Publication - Integrated ciliary gene compendium
• OMIM Usher Syndrome Entries - Known Usher genes

Secondary (MEDIUM confidence)

• Weighted Scoring Model Best Practices - Weight validation, calibration
• Gene Prioritization Validation - Positive control methodology
• Probabilistic Outlier Detection for RNA-seq - MAD-based robust outliers
• Missing Data in Large-Scale Assessments - Missing data handling strategies
• Quality Control for GWAS - Distribution checks, anomaly detection

Tertiary (LOW confidence - needs validation)

• Composite Scoring vs Imputation - Item-level vs composite-level imputation (medical context, not genomics)
• SYSCILIA website (https://www.syscilia.org/goldstandard.shtml) - May have direct download links (not verified in 2026)
• CiliaCarta web tool download - Requires verification of current URL and data format

Metadata

Confidence breakdown:

• Standard stack: HIGH - DuckDB, polars, scipy all verified from official docs and existing codebase
• Architecture patterns: HIGH - NULL preservation, weighted averaging, validation established in prior phases and literature
• Known gene sources: MEDIUM-HIGH - Publications verified, exact download URLs require implementation phase verification
• QC thresholds: MEDIUM - Based on genomic QC literature, but specific thresholds (50%/80% missing) are recommendations to be tuned
• Code examples: HIGH - Adapted from official DuckDB/scipy docs and verified patterns from prior phases

Research date: 2026-02-11 Valid until: 2026-04-11 (60 days) - Stable domain. Scoring methodology is well-established; DuckDB/scipy APIs stable; known gene lists update infrequently (SCGS v2 published 2021, no v3 announced). Re-verify OMIM Usher genes if new subtypes discovered.