just_dna_enricher.alphagenome_avi_build

just_dna_enricher.alphagenome_avi_build

The AlphaGenome AVI artifact, re-encoded as an operator-built snapshot (RM191).

Nothing here fetches, and that is not the usual inject-only rule — this is the network tier, so it is allowed to. The AVI artifact is 88.5 GB behind a sign-in whose eligibility clause bars whole classes of holder (ALPHAGENOME_ATLAS.md § 2.1), so acquisition is the operator's act and this builder reads the file they already have. --input is required for that reason and has no default URL to fall back on.

What it writes. <out>/data/alphagenome_avi-<contig>.parquet — one row per (contig, position, ref, alt) with raw_score as Int32 at a scale of 10⁵ — plus avi_knots.parquet, release.json and LICENSE.txt.

Four shape decisions, each measured rather than argued (§§ 1.4, 4.4–4.9 of the probe):

• Int32×10⁵, not Float32. Both published columns print at most 5 decimals, so an integer scale is exactly lossless while Float32 silently rounds the fifth — and is larger, 3.154 bytes/row against 2.413. Wherever a source publishes fixed decimals a float is the wrong container: its low mantissa bits are noise the source never had, and noise does not compress. The builder does not take that on trust — _scaled_scores re-derives the scaling and refuses a value that does not land on the grid, so losslessness is checked over every row written rather than sampled.

Read the integer, do not divide it back. The exactness is about the decimal: raw_score_e5 is the printed value shifted five places, and nothing is lost. Recovering a float with raw_score_e5 / 1e5 rounds a second time and lands one ulp off float(printed) on 53% of rows — measured, not feared. Compare thresholds in the integer domain (score >= 0.1 becomes raw_score_e5 >= 10_000) and the question never arises. * PHRED is not stored. Measured over all 8,812,917,339 rows it is an exact within-corpus rank (PHRED ≥ p keeps 10^(-p/10) of the corpus, to four significant figures across four decades), so it is 24.7 GB of a number that is a function of raw_score. The knot table carries the curve instead, in 466 KB. * The knot table is rebuilt here, not copied. docs/probes/alphagenome_knots/avi_knots.parquet is evidence; the lane's copy is the artifact, and sum(n) over its knots must equal the rows this build wrote or the two halves describe different data. * No threshold. The whole corpus is 34.2 GB measured keeping raw_score alone — inside any stated budget, with the sign intact. 49.30% of rows are negative and a negative AVI is low conservation, evidence against impact, not down-regulation (§ 4.7): abs() would discard what half the corpus says. A threshold is a consumer's slice, not this artifact's shape.

Absence is row-absence. AVI covers about 95% of the assembly and writes 672,931 genuine zeros, so a position with no row is unscored and a row with raw_score_e5 == 0 is scored zero. Collapsing the two would be @unreachable-not-absent at nine-billion-row scale.

AlphaGenomeBuildError

Bases: RuntimeError

The build could not produce a snapshot. Never raised for a legitimately empty contig.

ContigResult dataclass

ContigResult(contig: str, rows: int, path: Path)

One contig's parquet, and what went into it.

AviBuildResult dataclass

AviBuildResult(
    contigs: tuple[ContigResult, ...],
    rows: int,
    knots: int,
    negative_rows: int,
    zero_rows: int,
    source_path: Path,
    source_sha256: str | None,
    source_mtime: str | None,
    dataset: str | None,
)

What the build produced, with every number it computed rather than only the ones it used.

A count computed and discarded is a count every consumer recomputes (@dont-discard-computed), and here the recomputation costs a pass over 88.5 GB.

alts_for_ref

alts_for_ref(ref: str) -> tuple[str, str, str]

The three ALT bases a locus with this REF carries, in the order the columns hold them.

{A,C,G,T} − ref, ascending. This is the whole reason the wide layout costs nothing to read: a caller with (pos, ref, alt) finds its column as alts_for_ref(ref).index(alt), and a reader reconstructing long rows walks the tuple. No lookup table travels with the artifact.

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

def alts_for_ref(ref: str) -> tuple[str, str, str]:
    """The three ALT bases a locus with this REF carries, in the order the columns hold them.

    `{A,C,G,T} − ref`, ascending. This is the whole reason the wide layout costs nothing to read:
    a caller with `(pos, ref, alt)` finds its column as `alts_for_ref(ref).index(alt)`, and a reader
    reconstructing long rows walks the tuple. No lookup table travels with the artifact.
    """
    rest = tuple(b for b in BASES if b != ref.upper())
    if len(rest) != 3:
        raise AlphaGenomeBuildError(
            f"ref {ref!r} is not one of {BASES}, so it has no three-ALT complement. AVI is SNV-only "
            "and every REF it publishes is a single standard base."
        )
    return rest

check_input_is_the_avi_artifact

check_input_is_the_avi_artifact(source: Path) -> None

Refuse a file from the other two bulk artifacts, by name, before anything is read.

Not a courtesy check. The three artifacts are not one licence: only AVI is in the Permissive class, and the merged splicing and SHAP files are non-commercial-only. A lane whose SourceRow describes AVI, silently filled from the splicing file, would publish terms that do not govern its bytes — which is the failure @a-hosts-terms-are-not-its-contents-terms names from the other direction. The name is what an operator actually has; the schema check below catches the rest.

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

def check_input_is_the_avi_artifact(source: Path) -> None:
    """Refuse a file from the other two bulk artifacts, by name, before anything is read.

    Not a courtesy check. The three artifacts are **not one licence**: only AVI is in the Permissive
    class, and the merged splicing and SHAP files are non-commercial-only. A lane whose `SourceRow`
    describes AVI, silently filled from the splicing file, would publish terms that do not govern
    its bytes — which is the failure `@a-hosts-terms-are-not-its-contents-terms` names from the
    other direction. The name is what an operator actually has; the schema check below catches the
    rest.
    """
    stem = source.name
    for foreign in FOREIGN_ARTIFACTS:
        if foreign in stem:
            raise AlphaGenomeBuildError(
                f"{source.name} is not the AVI artifact. This lane reads "
                f"{AVI_FILENAME} and nothing else: the merged-splicing and feature-importance "
                "artifacts are a different licence class (non-commercial only, and the latter "
                "stacks AlphaMissense/Cactus/phastCons terms), so a snapshot mixing them would "
                "carry a SourceRow that does not describe its own bytes."
            )

list_contigs

list_contigs(source: Path) -> tuple[str, ...]

The contigs the tabix index knows, in the file's own order.

Read from the index rather than assembled from a chromosome list: what the artifact covers is a property of the artifact, and a hand-kept roster is how a build silently skips a contig.

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

def list_contigs(source: Path) -> tuple[str, ...]:
    """The contigs the tabix index knows, in the file's own order.

    Read from the index rather than assembled from a chromosome list: what the artifact covers is a
    property of the artifact, and a hand-kept roster is how a build silently skips a contig.
    """
    try:
        out = subprocess.run(["tabix", "-l", str(source)], capture_output=True, text=True, check=True).stdout
    except FileNotFoundError as exc:
        raise AlphaGenomeBuildError(
            "tabix is not on PATH. The AVI artifact is a bgzipped TSV with a .tbi index, and the "
            "build reads it one contig at a time through `tabix`; install htslib/samtools."
        ) from exc
    except subprocess.CalledProcessError as exc:
        raise AlphaGenomeBuildError(
            f"`tabix -l {source}` failed ({exc.returncode}): {(exc.stderr or '').strip()}. "
            f"Is {source.name}.tbi beside it?"
        ) from exc
    contigs = tuple(line for line in out.splitlines() if line.strip())
    if not contigs:
        raise AlphaGenomeBuildError(f"{source} has an index but no contigs in it")
    return contigs

to_long

to_long(wide)

One row per (chrom, pos, ref, alt) from the artifact's wide rows (RM197).

The inverse of the layout, and it needs nothing but ref. The snapshot stores one row per position with three score columns, and which base each column means is {A,C,G,T} − ref ascending — so long form is recoverable with no stored alt and no lookup table travelling beside the data. That is what made the wide layout adoptable rather than a schema break: a consumer who wants (pos, ref, alt, score) rows calls this and gets them.

It lives here rather than in the reader because it is a fact about the artifact, and three callers now need it to agree — alphagenome_check's join, the tests, and anyone converting a pulled snapshot back to long after download.

Apply it to a filtered frame. Over the whole corpus it is 2.9 billion loci becoming 8.8 billion rows, which is the shape the artifact exists to avoid storing.

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

def to_long(wide):
    """One row per `(chrom, pos, ref, alt)` from the artifact's wide rows (RM197).

    **The inverse of the layout, and it needs nothing but `ref`.** The snapshot stores one row per
    position with three score columns, and which base each column means is `{A,C,G,T} − ref`
    ascending — so long form is recoverable with no stored `alt` and no lookup table travelling
    beside the data. That is what made the wide layout adoptable rather than a schema break: a
    consumer who wants `(pos, ref, alt, score)` rows calls this and gets them.

    It lives here rather than in the reader because it is a fact about the **artifact**, and three
    callers now need it to agree — `alphagenome_check`'s join, the tests, and anyone converting a
    pulled snapshot back to long after download.

    Apply it to a **filtered** frame. Over the whole corpus it is 2.9 billion loci becoming 8.8
    billion rows, which is the shape the artifact exists to avoid storing.
    """
    _require_polars()
    parts = [
        wide.filter(pl.col("ref").cast(pl.String) == base).select(
            "chrom",
            "pos",
            pl.col("ref").cast(pl.String),
            pl.lit(alt).alias("alt"),
            pl.col(f"alt{i}").alias("raw_score_e5"),
        )
        for base in BASES
        for i, alt in enumerate(alts_for_ref(base))
    ]
    present = [p for p in parts if p.height]
    return pl.concat(present).sort(["pos", "alt"]) if present else wide.head(0)

use_restrictions_text

use_restrictions_text(terms_file: Path) -> str

The Output Terms' "Use restrictions" section, verbatim, for LICENSE.txt.

Restriction 3b requires it to travel inside a derivative rather than as a link: anyone who attaches their own terms — which a module's sources.csv is — must carry this section "as an enforceable provision". So the bytes go beside the data, the way ClinPGx's bundled LICENSE.txt already does (SNAPSHOT_LICENSE_FILENAME exists for exactly this).

Sliced from the pinned extraction rather than paraphrased, and bounded at the next heading so the disclaimer sections do not ride along under a name that says "use restrictions".

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

def use_restrictions_text(terms_file: Path) -> str:
    """The Output Terms' **"Use restrictions" section**, verbatim, for `LICENSE.txt`.

    Restriction 3b requires it to travel *inside* a derivative rather than as a link: anyone who
    attaches their own terms — which a module's `sources.csv` is — must carry this section "as an
    enforceable provision". So the bytes go beside the data, the way ClinPGx's bundled `LICENSE.txt`
    already does (`SNAPSHOT_LICENSE_FILENAME` exists for exactly this).

    Sliced from the pinned extraction rather than paraphrased, and bounded at the next heading so
    the disclaimer sections do not ride along under a name that says "use restrictions".
    """
    text = terms_file.read_text()
    start = text.index("Use restrictions")
    end = text.index("Disclaimers and limitations of liability", start)
    section = text[start:end].strip()
    if "commercial organization" not in section:  # pragma: no cover - upstream restructured
        raise AlphaGenomeBuildError(
            f"{terms_file.name} no longer carries the Use restrictions clauses where this builder "
            "slices them; re-pin the document and re-read § 2.7 before shipping a LICENSE.txt."
        )
    return section + "\n"

build_snapshot

build_snapshot(
    source: Path,
    out_dir: Path,
    *,
    contigs: Sequence[str] | None = None,
    workers: int = 12,
    chunk_bytes: int = CHUNK_BYTES,
    hash_source: bool = True,
    terms_file: Path | None = None,
) -> AviBuildResult

Re-encode the AVI artifact into out_dir, and rebuild the knot table from the same pass.

contigs restricts the build — a test builds one, an operator builds all of them. workers is the tabix fan-out: the measured passes ran 24 contigs in 41–46 minutes at twelve, and a single-threaded pass over 8.8 billion rows takes roughly four times as long.

Threads rather than processes because every worker spends its time in a tabix subprocess and in polars, both of which are outside the GIL; a process pool would buy nothing and cost the chunk-sized copies.

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

def build_snapshot(
    source: Path,
    out_dir: Path,
    *,
    contigs: Sequence[str] | None = None,
    workers: int = 12,
    chunk_bytes: int = CHUNK_BYTES,
    hash_source: bool = True,
    terms_file: Path | None = None,
) -> AviBuildResult:
    """Re-encode the AVI artifact into `out_dir`, and rebuild the knot table from the same pass.

    `contigs` restricts the build — a test builds one, an operator builds all of them. `workers` is
    the tabix fan-out: the measured passes ran 24 contigs in 41–46 minutes at twelve, and a
    single-threaded pass over 8.8 billion rows takes roughly four times as long.

    Threads rather than processes because every worker spends its time in a `tabix` subprocess and
    in polars, both of which are outside the GIL; a process pool would buy nothing and cost the
    chunk-sized copies.
    """
    source = Path(source)
    if not source.is_file():
        raise AlphaGenomeBuildError(f"{source} is not a file. This lane never downloads (§ 2.1).")
    check_input_is_the_avi_artifact(source)

    _require_polars()
    data_dir = out_dir / SNAPSHOT_DATA_DIRNAME
    data_dir.mkdir(parents=True, exist_ok=True)

    wanted = tuple(contigs) if contigs is not None else list_contigs(source)
    if contigs is not None:
        known = set(list_contigs(source))
        unknown = [c for c in wanted if c not in known]
        if unknown:
            raise AlphaGenomeBuildError(
                f"{', '.join(unknown)} not in {source.name}'s index. It knows: {', '.join(sorted(known))}"
            )

    results: list[ContigResult] = []
    knot_frames = []
    with ThreadPoolExecutor(max_workers=max(1, workers)) as pool:
        for contig_result, knots in pool.map(
            lambda c: _build_contig(source, c, data_dir, chunk_bytes=chunk_bytes), wanted
        ):
            results.append(contig_result)
            knot_frames.append(knots)

    knot_table = (
        pl.concat(knot_frames)
        .group_by("raw_score_e5")
        .agg(
            pl.col("n").sum().alias("n"),
            pl.col("phred_lo").min().alias("phred_lo"),
            pl.col("phred_hi").max().alias("phred_hi"),
        )
        .sort("raw_score_e5")
        .select(list(KNOT_COLUMNS))
    )
    knot_table.write_parquet(out_dir / KNOT_FILENAME, compression="zstd", compression_level=9)
    # The per-contig aggregates have served their purpose. Removed only after the merged table is on
    # disk, so a crash anywhere before this point leaves the expensive half of the build recoverable.
    shutil.rmtree(out_dir / KNOT_PARTS_DIRNAME, ignore_errors=True)

    rows = sum(c.rows for c in results)
    reconciled = int(knot_table["n"].sum()) if knot_table.height else 0
    if reconciled != rows:
        raise AlphaGenomeBuildError(
            f"the knot table describes {reconciled} rows but the build wrote {rows}. The curve and "
            "the data are two halves of one artifact; a mismatch means they came from different "
            "passes and neither can be trusted for a threshold."
        )

    negative = int(knot_table.filter(pl.col("raw_score_e5") < 0)["n"].sum()) if knot_table.height else 0
    zeros = int(knot_table.filter(pl.col("raw_score_e5") == 0)["n"].sum()) if knot_table.height else 0

    stamp = _artifact_stamp(source)
    result = AviBuildResult(
        contigs=tuple(results),
        rows=rows,
        knots=knot_table.height,
        negative_rows=negative,
        zero_rows=zeros,
        source_path=source,
        source_sha256=_sha256_file(source) if hash_source else None,
        source_mtime=stamp,
        dataset=stamp[:10] if stamp else None,
    )
    _write_release_json(out_dir, result, artifact_stamp=stamp)

    terms = terms_file or _default_terms_file()
    if terms is not None and terms.is_file():
        atomic_write_text(out_dir / SNAPSHOT_LICENSE_FILENAME, use_restrictions_text(terms))
    else:
        logger.warning(
            "no %s written: the Output Terms extraction was not found at %s. Restriction 3b "
            "requires the Use restrictions section to travel inside a derivative, so a snapshot "
            "without it cannot be passed on.",
            SNAPSHOT_LICENSE_FILENAME,
            terms,
        )
    return result