just_dna_enricher.expression

just_dna_enricher.expression

Fill expression_effects.csv from AlphaGenome's Atlas API (0.7, RM194 + RM200).

One row per (variant, gene): which way a variant moves that gene's predicted expression, how many of the scorer's 371 tissue tracks agree, and how far the variant sits from the gene. A recording pass, not a check — it is the source of the rows rather than a judge of them, so it writes no verification record and adds no member to VALID_VERIFICATION_CHECKS, exactly as gwas.py does not.

RNA_SEQ is the only scorer here, and RM200 measured all twenty-two to say so. It is the only one with a gene axis, the only one whose disagreement across tracks is meaningful rather than flat, and the only one that attributes its own claim to a gene — which is what @gene-map-is-another-sources-attribution requires. The rest either restate AVI_SCORE or rank noise: CAGE's top-5 of 546 tracks carry 2-7% of the effect, running backwards to effect size.

Two ways to name the interval, and the gene filter is mandatory in both. --chrom/--start/--end wins when given; otherwise the span comes from the MANE lane widened by the model's measured ±512 kb horizon (gene_spans). But gene is required either way, because the server-side gene filter is not an optimisation: unfiltered, a 32 bp interval answers with 43 MB against grpc's 4 MB receive limit, and atlas_client.score_interval refuses such a request before sending it.

Distance is recorded because a threshold without it is wrong. Scores 100-500 kb out run about an order of magnitude lower than scores at the gene, so a flat --min-score silently keeps only the proximal variants — which is the failure RM194 exists to prevent, not a tuning detail. Distance needs the gene's span, so the MANE lane is consulted even when the interval was supplied by hand; when it is absent, distance_to_gene is null and the pass says so rather than substituting an interval edge.

Non-commercial, and the gate runs before the RPC. Atlas Output is commercial_use=False — the Output Terms say so in their opening sentence — so check_declared_use gates the fetch (@acquisition-gate-is-not-a-read-gate). With no --use declared that returns a skip, which means a no-flag run writes nothing and explains why; every documented invocation carries --use non-commercial for that reason.

ExpressionError

Bases: RuntimeError

This pass could not do its job. The local half — bad input, an unreadable sidecar, a refusal.

ExpressionUnavailable

Bases: ExpressionError

The Atlas did not answer. Retryable, and its own type rather than a bare ExpressionError.

@client-exception-contract: a pass owes its own type too, and translating to an *Unavailable subclass is what lets a caller tell "the service is down, try later" from "your input is wrong". gwas.py is exempt from that guard only because its client and its pass share one error type declared in one module; AtlasError lives in atlas_client, a foreign module, so this pass is not exempt and does not claim to be.

ExpressionResult dataclass

ExpressionResult(
    rows: list[ExpressionEffectRow] = list(),
    gene: str | None = None,
    interval: tuple[str, int, int] | None = None,
    span: GeneSpan | None = None,
    candidates: int = 0,
    written: int = 0,
    withheld: dict[str, int] = dict(),
    warnings: list[str] = list(),
    dataset: str | None = None,
    skipped: bool = False,
)

What one run did, with every admitted candidate accounted for.

accounts_for_every_candidate

accounts_for_every_candidate() -> bool

Every admitted candidate was either written or counted under a named reason.

Source code in enricher/src/just_dna_enricher/expression.py

def accounts_for_every_candidate(self) -> bool:
    """Every admitted candidate was either written or counted under a named reason."""
    return self.candidates == self.written + sum(self.withheld.values())

enrich_expression

enrich_expression(
    spec_dir: Path,
    gene: str,
    *,
    chrom: str | None = None,
    start: int | None = None,
    end: int | None = None,
    client=None,
    mane_cache: Path | None = None,
    min_score: float | None = None,
    max_rows: int = DEFAULT_MAX_ROWS,
    declared_use: str = "unstated",
    dataset: str | None = None,
    offline: bool = False,
    write: bool = True,
) -> ExpressionResult

Query one gene's interval and merge the answers into expression_effects.csv.

Merge-not-clobber on (variant_key, gene): an existing row is authoritative and a re-run gap-fills. Delete the file to re-derive it, which costs nothing because no authored judgement lives in it.

Source code in enricher/src/just_dna_enricher/expression.py

def enrich_expression(
    spec_dir: Path,
    gene: str,
    *,
    chrom: str | None = None,
    start: int | None = None,
    end: int | None = None,
    client=None,
    mane_cache: Path | None = None,
    min_score: float | None = None,
    max_rows: int = DEFAULT_MAX_ROWS,
    declared_use: str = "unstated",
    dataset: str | None = None,
    offline: bool = False,
    write: bool = True,
) -> ExpressionResult:
    """Query one gene's interval and merge the answers into `expression_effects.csv`.

    Merge-not-clobber on `(variant_key, gene)`: an existing row is authoritative and a re-run
    gap-fills. Delete the file to re-derive it, which costs nothing because no authored judgement
    lives in it.
    """
    spec_dir = Path(spec_dir)
    result = ExpressionResult(gene=gene)
    output_path = sidecar_path(spec_dir, SIDECAR_NAME, error=ExpressionError)
    if write:
        # Fail on a placeholder or half-edited licence table now, before the fetch (S98, RM231).
        require_sources_file(spec_dir, error=ExpressionError)

    # **`offline` outranks an injected client** (RM220), which is `pgx`'s reading and not `gwas`'s.
    # The two shapes coexist in this tier and the difference is the source's licence, not a
    # preference: the GWAS Catalog is ungated, so handing over a transport you already hold really is
    # not egress. The Atlas is **not** — its Additional Terms bar classes of holder outright — so a
    # live client under a flag documented as making no egress is exactly the loophole RM38 closed for
    # the PGx sources, and `test_pgx_licensing.py` already asserts the strict reading by name. This
    # gated `offline and client is None`, so an injected client fetched from a gated source under
    # `--offline`. `@flag-means-same`.
    if offline:
        note = (
            "expression pass skipped: --offline. This pass reads the Atlas API and has no snapshot "
            "to fall back on, so it is a no-op offline rather than a failure. An injected client "
            "does not override it: the Atlas is licence-gated, so --offline means no egress."
        )
        result.warnings.append(note)
        result.skipped = True
        return result

    # The gate runs before anything is fetched (`@acquisition-gate-is-not-a-read-gate`). With
    # `commercial_use=False` an undeclared run is a SKIP, not permission — the tool must not assert a
    # purpose on the operator's behalf.
    declared_use, declared_from = effective_declared_use(
        spec_dir, ALPHAGENOME_ATLAS_TERMS, declared_use
    )  # S105
    refusal = check_declared_use(ALPHAGENOME_ATLAS_TERMS, declared_use)
    if declared_from is not None:
        result.warnings.append(
            f"alphagenome: use {declared_use!r} read from {declared_from}, recorded by an earlier run; pass "
            f"--use to declare otherwise."
        )
    if refusal is not None:
        result.warnings.append(refusal)
        result.skipped = True
        return result

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

    interval = _resolve_interval(gene, chrom, start, end, mane_cache, result)
    if interval is None:
        return result
    result.interval = interval
    logger.warning("expression pass: %s", _cost_note(interval))

    # Every coordinate writer owes this (`@restamp-for-build` one layer up): AlphaGenome publishes
    # GRCh38 only, and a GRCh37 module would silently receive GRCh38 positions. Reported, never
    # repaired and never refused.
    mismatch = source_build_mismatch(spec_dir, SOURCE_NAME, source_build=SOURCE_BUILD)
    if mismatch:
        result.warnings.append(mismatch)

    release = dataset or f"alphagenome_atlas_{now_utc_iso()[:10]}"
    fetched_at = now_utc_iso()
    result.dataset = release

    if client is None and not ATLAS_CLIENT_AVAILABLE:
        # The ONE fix that applies, never both (`@specific-rejection`). This used to name the extra
        # and the generator in one sentence, which sent a wheel user missing `[atlas]` off to
        # `atlas generate` and into "grpcio-tools is not installed" — a third error about a fourth
        # thing, when their fix was one `pip install`.
        raise ExpressionError(f"no Atlas client, so nothing was queried: {client_absence()}")
    try:
        scores = (client or _connect()).score_interval(
            interval[0], interval[1], interval[2], scorers=(SCORER,), gene_names=(gene,)
        )
    except AtlasUnavailable as exc:
        raise ExpressionUnavailable(
            f"the Atlas did not answer for {gene} at {interval[0]}:{interval[1]}-{interval[2]}: {exc}"
        ) from exc
    except AtlasError as exc:
        raise ExpressionError(
            f"the Atlas refused the query for {gene} at {interval[0]}:{interval[1]}-{interval[2]}: {exc}"
        ) from exc

    if not scores:
        # `@unreachable-not-absent`: write no row, name it separately. An interval beyond the model's
        # ±512 kb horizon and a symbol the service does not attribute to anything both look like this
        # from here, so the warning names both rather than picking one.
        note = (
            f"{gene}: the Atlas attributed no scored variant in "
            f"{interval[0]}:{interval[1]}-{interval[2]}. Either the interval lies beyond the model's "
            f"±512 kb attribution horizon for this gene, or the service does not know the symbol — "
            f"this surface cannot tell those apart, and neither is written as a row."
        )
        result.warnings.append(note)
        return result

    seen = {merge_key(row) for row in existing}
    out = list(existing)

    for score in scores:
        result.candidates += 1
        values, refused, gene_id = _gene_block(score)
        if refused is not None:
            result.withhold(refused)
            continue
        effect_size, direction, agreeing, total = _summarise(values)
        if min_score is not None and (effect_size is None or abs(effect_size) < min_score):
            result.withhold("below_min_score")
            continue

        # The artifact ships `chr22`; a `VariantRow` stores `22`. RM193 found this exact join
        # matching nothing while looking like an uncovered region, so the strip happens here.
        contig = str(score.chrom).removeprefix("chr")
        row = ExpressionEffectRow(
            variant_key=derive_variant_key(
                None, contig, score.position, score.ref, score.alt, build=SOURCE_BUILD
            ),
            chrom=contig,
            start=score.position,
            ref=score.ref,
            alt=score.alt,
            gene=gene,
            gene_id=gene_id,
            effect_size=effect_size,
            effect_measure=SCORER,
            effect_direction=direction,
            tracks_agreeing=agreeing,
            tracks_total=total,
            distance_to_gene=_distance(result.span, score.position),
            dataset=release,
            source=SOURCE_NAME,
            status="resolved",
            fetched_at=fetched_at,
        )
        key = merge_key(row)
        if key in seen:
            result.withhold("already_recorded")
            continue
        seen.add(key)
        out.append(row)
        result.written += 1

    if len(out) > max_rows:
        raise ExpressionError(
            f"{len(out):,} rows would be written to {SIDECAR_NAME}, over the {max_rows:,} cap. "
            f"A CSV sidecar of that size is not a sidecar, and truncating silently would be worse "
            f"than refusing. Raise the bar with --min-score (distal scores run ~10x lower than "
            f"scores at the gene, so a distance-aware bar keeps more of what matters), narrow the "
            f"interval, or raise --max-rows deliberately."
        )

    out.sort(key=_sort_key)
    if write and result.written:
        # The licence row lands inside the table's commit (S98, RM231): 12,003 non-commercial rows
        # once sat on disk under `FAILED` with no licence record, because the row was merged after
        # the write and a scaffold's placeholder row made the merge refuse.
        _write_csv(
            out,
            output_path,
            before_commit=lambda: merge_sources_file(
                [ALPHAGENOME_ATLAS_TERMS.row(SOURCE_LAYER, declared_use=declared_use, dataset=release)],
                spec_dir,
                error=ExpressionError,
            ),
        )
    result.rows = out
    return result