Skip to content

just_dna_enricher.gwas

just_dna_enricher.gwas

The NHGRI-EBI GWAS Catalog pass — published effect sizes as a derived fact table (0.6, RM90).

Fills gwas_effects.csv for the variants variants.csv names, from the Catalog's REST API. What it does not do is fill weight: a consumer asked for exactly that and it is barred — MODULE_LIFECYCLE § Stage 3 names weight/direction/effect_size in the cells no tool fills, and every check in this tier reports rather than repairs. The effect lands in its own table beside the authored column, and a consumer chooses between them wholesale.

Everything below was established by probing the real API on 2026-08-17, not from its documentation, and three findings shaped the pass:

  • The association payload is thin. It carries the effect, the p-value and the risk allele, and nothing else — pmid, study_accession, ancestry, trait and trait_efo_id all sit behind _links.study and _links.efoTraits. So a complete row costs two extra requests, and the pass costs 1 + 2N per variant with N associations.

_LinkCache memoizes by resolved URL, and the measurement corrected the prediction that motivated it. The expectation was that associations share studies heavily, so caching would collapse most of the cost. Run against reference_examples/hfe_hemochromatosis, it saved nothing: rs1800562 alone carries 189 associations and each names its own study, so the real figure was 382 requests and 0 cache hits. The cache stays — it costs a dict and it does pay on a module whose variants share literature — but the honest budget is 1 + 2N, and --no-study-facts exists because of that measurement rather than in anticipation of it. The pass reports both halves so an operator sees the real number instead of this docstring's guess. * riskAlleleName is rs4149056-C, or rs4149056-? when the study never established which allele carries the effect — 2 of the first 3 associations for that variant. -? becomes None, never a guess, and the row is still written: an effect relative to an unknown allele is real evidence that cannot be used as a weight, and dropping it would hide that from a consumer who needs to know. * betaUnit is free text and frequently uninterpretable — umol/l on one association and unit on two others for the same variant. It is stored verbatim including the useless values, because "these betas are on unknown and possibly different scales" is the fact a consumer must have.

Rate limits are NOT established. EBI publishes no numeric budget that could be found, which makes this pass unlike gnomad.py, whose 10-requests-per-60-seconds is real, documented and load-bearing. The gate here is a conservative default rather than a transcribed limit, and it is spelled that way so nobody later "corrects" it against a number that does not exist.

GwasError

Bases: RuntimeError

A GWAS Catalog fetch or parse failed in a way the caller must see.

Every transport and parse failure is translated into this before it leaves the module. A client that leaks httpx.HTTPError has no contract — the caller would have to import the transport library to catch it, and a swap of transport would silently become a breaking change.

GwasNotFound

Bases: GwasError

The Catalog answered, and it has no record of this variant at all (HTTP 404).

A subclass rather than a flag, because this is an answer and not a failure. Found by running the pass against a real module: the Catalog holds only variants with a published genome-wide association, so rs111033563 — a rare clinical HFE variant — 404s rather than returning an empty list. The first version of this pass treated that as an outage and died on the first rare variant in the module, which is every clinically-authored module.

It stays inside GwasError so a caller that does not care about the distinction still catches one type; associations_for catches it and reports the empty answer, which the pass records as not_found. What must never happen is the reverse — a transport failure read as "no associations", which would write a confident negative about a variant nobody could ask about.

GwasResult dataclass

GwasResult(
    rows: list[GwasEffectRow] = list(),
    covered: list[str] = list(),
    missing: list[str] = list(),
    requests_made: int = 0,
    requests_saved: int = 0,
    p_value_underflows: int = 0,
    unusable: int = 0,
    skipped_offline: bool = False,
)

What one pass did, including what it spent.

GwasCatalogClient dataclass

GwasCatalogClient(
    endpoint: str = DEFAULT_GWAS_ENDPOINT,
    timeout: float = 30.0,
    gate: PacingGate = (
        lambda: PacingGate(DEFAULT_REQUEST_INTERVAL)
    )(),
    _client: Client | None = None,
)

Thin REST client for the Catalog, paced and retried.

endpoint and gate are injectable so a test can drive the real parsing code against recorded payloads without a network or a real sleep — the same shape EnsemblResolver and the gnomAD client use.

associations_for

associations_for(rsid: str) -> list[dict]

Every association the Catalog holds for one rsID.

Three outcomes, not two — the S20 shape. A non-empty list is an answer; [] is also an answer (the Catalog was reached and holds nothing for this variant, which the caller records as not_found); and a raised GwasError means it could not be asked, which is unchecked rather than empty. Fusing the last two would turn a failed request into a definite negative.

Source code in enricher/src/just_dna_enricher/gwas.py
def associations_for(self, rsid: str) -> list[dict]:
    """Every association the Catalog holds for one rsID.

    **Three outcomes, not two** — the S20 shape. A non-empty list is an answer; `[]` is *also* an
    answer (the Catalog was reached and holds nothing for this variant, which the caller records
    as `not_found`); and a raised `GwasError` means it could not be asked, which is unchecked
    rather than empty. Fusing the last two would turn a failed request into a definite negative.
    """
    url = f"{self.endpoint}/singleNucleotidePolymorphisms/{rsid}/associations"
    try:
        payload = self._get(url)
    except GwasNotFound:
        # The Catalog holds only variants with a published association, so it 404s on a rare
        # clinical variant rather than returning an empty list. That is the empty ANSWER.
        return []
    embedded = payload.get("_embedded") or {}
    associations = embedded.get("associations")
    if associations is None:
        # A 200 with neither `_embedded.associations` nor an empty list is a shape we do not
        # model. Raising beats returning `[]`, which would record a definite "no associations".
        raise GwasError(
            f"GWAS Catalog response for {rsid} carries no `_embedded.associations` — the API "
            f"shape has changed and this pass would otherwise record a false negative"
        )
    return list(associations)

follow

follow(url: str) -> dict

A linked sub-resource, or {} when the Catalog has none.

A 404 here withholds the study/trait facts for one association rather than sinking the pass: the association itself is still a real published effect, and dropping it because its study record moved would lose evidence over metadata.

Source code in enricher/src/just_dna_enricher/gwas.py
def follow(self, url: str) -> dict:
    """A linked sub-resource, or `{}` when the Catalog has none.

    A 404 here withholds the study/trait facts for one association rather than sinking the pass:
    the association itself is still a real published effect, and dropping it because its study
    record moved would lose evidence over metadata.
    """
    try:
        return self._get(url)
    except GwasNotFound:
        logger.warning("GWAS Catalog has no record at %s; that association keeps null study facts", url)
        return {}

enrich_gwas

enrich_gwas(
    spec_dir: Path,
    *,
    mode: str = "best_effort",
    offline: bool = False,
    write: bool = True,
    client: GwasCatalogClient | None = None,
    dataset: str | None = None,
    declared_use: str = "unstated",
    study_facts: bool = True,
) -> GwasResult

Record the Catalog's published effect sizes for this module's variants.

Existing rows are authoritative and merged, never clobbered — the standing rule for every pass, with the standing consequence: to regenerate after a machinery change, delete the file first.

offline makes the pass a no-op with a warning rather than a failure: the Catalog publishes a bulk download, but this pass reads the REST API, and there is no snapshot for it to fall back on. An injected client still wins, because handing over a transport you already hold is not egress.

A variant the Catalog holds nothing for gets a not_found row, unlike enrich_gene_validity's silence, and the difference is real rather than stylistic: a curating body's silence means nobody has assessed the gene, whereas the Catalog's empty answer means no genome-wide association has been published for this variant — which is a fact about it, and one a consumer weighing an authored weight against the literature wants to see.

mode is the severity ladder the CLI's --strict sets, and it is deliberately not about missing: see the block above the raise at the end of this function for why the Catalog's empty answer is a fact rather than a shortfall, and what strict reads instead.

Source code in enricher/src/just_dna_enricher/gwas.py
def enrich_gwas(
    spec_dir: Path,
    *,
    mode: str = "best_effort",
    offline: bool = False,
    write: bool = True,
    client: GwasCatalogClient | None = None,
    dataset: str | None = None,
    declared_use: str = "unstated",
    study_facts: bool = True,
) -> GwasResult:
    """Record the Catalog's published effect sizes for this module's variants.

    Existing rows are authoritative and merged, never clobbered — the standing rule for every pass,
    with the standing consequence: to regenerate after a machinery change, delete the file first.

    `offline` makes the pass a no-op with a warning rather than a failure: the Catalog publishes a
    bulk download, but this pass reads the REST API, and there is no snapshot for it to fall back on.
    An injected `client` still wins, because handing over a transport you already hold is not egress.

    A variant the Catalog holds nothing for gets a **`not_found` row**, unlike
    `enrich_gene_validity`'s silence, and the difference is real rather than stylistic: a curating
    body's silence means nobody has assessed the gene, whereas the Catalog's empty answer means no
    genome-wide association has been published for this variant — which *is* a fact about it, and one
    a consumer weighing an authored weight against the literature wants to see.

    `mode` is the severity ladder the CLI's `--strict` sets, and it is deliberately **not** about
    `missing`: see the block above the raise at the end of this function for why the Catalog's empty
    answer is a fact rather than a shortfall, and what `strict` reads instead.
    """
    spec_dir = Path(spec_dir)
    output_path = sidecar_path(spec_dir, "gwas_effects.csv", error=GwasError)
    if write:
        # Fail on a placeholder or half-edited licence table now, before the fetch (S98, RM231).
        require_sources_file(spec_dir, error=GwasError)

    if offline and client is None:
        logger.warning(
            "GWAS pass skipped: --offline. This pass reads the Catalog's REST API and has no "
            "snapshot to fall back on, so it is a no-op offline rather than a failure."
        )
        return GwasResult(skipped_offline=True)

    existing_rows: list[GwasEffectRow] = []
    if output_path.exists():
        parsed, errors, _ = load_csv_rows(output_path, GwasEffectRow, output_path.name)
        if errors:
            raise GwasError(f"existing {output_path.name} is invalid: {errors[0]}")
        existing_rows = parsed

    subjects = _module_subjects(spec_dir)
    if not subjects:
        logger.warning(
            "GWAS pass has no subjects: no row in variants.csv carries an rsID, and the Catalog is "
            "queried by rsID. Nothing was fetched and nothing was written."
        )
        return GwasResult(rows=existing_rows)

    owned = client is None
    # `try/finally`, like every sibling pass (RM100). The close used to be a bare
    # `if client is None: catalog.close()` after ~80 lines of fetching and writing, so any
    # exception in between -- a `GwasError`, a sidecar collision, a write failure -- leaked the
    # httpx client. `frequencies`, `gene_metrics` and `enrich` all had this shape already.
    try:
        catalog = client or GwasCatalogClient()
        cache = _LinkCache(fetch=catalog.follow)
        release = dataset or f"gwas_catalog_{now_utc_iso()[:10]}"
        fetched_at = now_utc_iso()

        seen = {_merge_key(row) for row in existing_rows}
        out: list[GwasEffectRow] = list(existing_rows)
        covered: list[str] = []
        missing: list[str] = []
        unusable = 0
        direct_requests = 0
        underflows: list[str] = []

        for rsid, variant_key in subjects:
            associations = catalog.associations_for(rsid)
            direct_requests += 1
            if not associations:
                missing.append(rsid)
                row = GwasEffectRow(
                    association_id=f"{rsid}:not_found",
                    variant_key=variant_key,
                    rsid=rsid,
                    dataset=release,
                    source=GWAS_SOURCE,
                    status="not_found",
                    fetched_at=fetched_at,
                )
                if _merge_key(row) not in seen:
                    seen.add(_merge_key(row))
                    out.append(row)
                continue
            covered.append(rsid)
            for association in associations:
                built = _build_row(
                    association,
                    rsid=rsid,
                    variant_key=variant_key,
                    cache=cache,
                    release=release,
                    fetched_at=fetched_at,
                    underflows=underflows,
                    study_facts=study_facts,
                )
                if built is None:
                    unusable += 1
                    continue
                key = _merge_key(built)
                if key in seen:
                    continue
                seen.add(key)
                out.append(built)

        if underflows:
            # Aggregated by reason, once. On a well-studied variant this fires dozens of times, and the
            # rows are all still there — only the queryable number is withheld.
            logger.warning(
                "GWAS pass withheld p_value_num on %d association(s) whose p-value the Catalog reports "
                "below float64's range (it publishes 0.0); the verbatim p_value string is kept and the "
                "associations are recorded in full.",
                len(underflows),
            )
        if unusable:
            # Aggregated by reason, once — a per-row warning over a well-studied variant's dozens of
            # associations is a wall nobody reads.
            logger.warning(
                "GWAS pass skipped %d association(s) the Catalog published without an id this pass "
                "could key on; every usable association was still recorded.",
                unusable,
            )

        out.sort(key=_sort_key)
        result = GwasResult(
            rows=out,
            covered=covered,
            missing=missing,
            requests_made=direct_requests + cache.misses,
            requests_saved=cache.hits,
            p_value_underflows=len(underflows),
            unusable=unusable,
        )

        if write:
            # The licence row lands inside the table's commit (S98, RM231).
            _write_gwas_csv(
                out,
                output_path,
                before_commit=lambda: merge_sources_file(
                    [GWAS_CATALOG_TERMS.row("gwas_effect", declared_use=declared_use, dataset=release)],
                    spec_dir,
                    error=GwasError,
                ),
            )
    finally:
        if owned:
            catalog.close()

    # **What `mode` means here, and why it is not `missing`** (RM100). The parameter was
    # accepted and never read while the CLI advertised `--strict` as a severity ladder, so the
    # flag was inert. Every sibling pass escalates on `result.missing`, and that would be wrong
    # here for the reason this pass's own docstring gives: the Catalog holding nothing for a
    # variant is a FACT about the variant, recorded as a `not_found` row, and true of most
    # variants. Escalating it would refuse nearly every module and say nothing.
    #
    # What `strict` does escalate is the pass failing to hold what the source DID say: an
    # association served without an id to key on, and a p-value below float64 whose queryable
    # number is therefore withheld. Both are already warned about in `best_effort`, aggregated
    # by reason; `strict` means a reproducible artifact, and both of these are places where the
    # artifact is missing something the Catalog published.
    if mode == "strict" and (result.unusable or result.p_value_underflows):
        raise GwasError(
            f"strict GWAS enrichment: {result.unusable} association(s) were served without an "
            f"id this pass can key on and {result.p_value_underflows} carried a p-value below "
            f"float64's range, so the artifact does not hold everything the Catalog published. "
            f"Both are the Catalog's shape rather than an authoring mistake -- use "
            f"mode='best_effort' to record what is holdable and read the warnings."
        )
    return result

gwas_source_row

gwas_source_row(
    *,
    declared_use: str = "unstated",
    dataset: str | None = None,
) -> SourceRow

The licence row this pass writes. Exposed so a caller can record the terms without fetching.

Source code in enricher/src/just_dna_enricher/gwas.py
def gwas_source_row(*, declared_use: str = "unstated", dataset: str | None = None) -> SourceRow:
    """The licence row this pass writes. Exposed so a caller can record the terms without fetching."""
    return GWAS_CATALOG_TERMS.row("gwas_effect", declared_use=declared_use, dataset=dataset)