echtvar-hail 0.4.0

Pure-Python and Hail-expression implementations of the echtvar var32 and kmer16 variant encodings

echtvar-hail

echtvar's variant encodings in Hail, and an export that writes them to Parquet instead of echtvar's zip container.

Two things live here:

• var32 / kmer16 — Python ports of src/lib/{var32,kmer16}.rs, each with an hl_* twin so the encodings run inside a Hail pipeline rather than by collecting rows to Python.
• parquet — an export in echtvar's per-(chrom, chunk) shape, as standard Parquet. The archive becomes a self-describing columnar file that DuckDB, Polars, pandas or Arrow can read without echtvar.

Install

pip install echtvar-hail

Python 3.11 or 3.12. Hail brings PySpark, which needs a JVM — Hail is built and tested against Java 11.

Quick start

One Hail table, one field config, covering every shape Field supports — copied as-is, quantized, renamed, and a nested struct — plus the routing that splits variants across the two output tables.

import hail as hl
import pyarrow.parquet as pq
from echtvar_hail import var32
from echtvar_hail.parquet import Field, check_precisions, export_parquet, finalize_parquet

hl.init()

ht = hl.Table.parallelize(
    [
        {"locus": hl.Locus("chr1", 100, "GRCh38"), "alleles": ["A", "T"],
         "filters": "PASS", "af": 0.0031, "cadd_raw": 3.21,
         "vep": [hl.Struct(consequence="missense_variant", biotype="protein_coding", polyphen=0.91),
                 hl.Struct(consequence="intron_variant", biotype="lncRNA", polyphen=0.05)]},
        {"locus": hl.Locus("chr1", 200, "GRCh38"), "alleles": ["AC", "A"],
         "filters": "PASS", "af": 0.41, "cadd_raw": -1.44,
         "vep": [hl.Struct(consequence="intron_variant", biotype="protein_coding", polyphen=0.01)]},
        {"locus": hl.Locus("chr1", 300, "GRCh38"), "alleles": ["A", "ACGTACGT"],
         "filters": "low_qual", "af": 0.0001, "cadd_raw": 22.6,
         "vep": [hl.Struct(consequence="stop_gained", biotype="protein_coding", polyphen=1.0)]},
    ],
    hl.tstruct(
        locus=hl.tlocus("GRCh38"), alleles=hl.tarray(hl.tstr),
        filters=hl.tstr, af=hl.tfloat64, cadd_raw=hl.tfloat64,
        vep=hl.tarray(hl.tstruct(consequence=hl.tstr, biotype=hl.tstr, polyphen=hl.tfloat64)),
    ),
)

fields = [
    Field("filters"),                                  # plain: copied as-is
    Field("af", precision=1e-6),                        # quantized: to a stated accuracy
    Field("cadd_raw", alias="cadd", precision=0.01),    # quantized + renamed
    Field("vep", fields=(                               # struct: member-by-member
        Field("consequence", categorical=True),
        Field("biotype", categorical=True),
        Field("polyphen", precision=1e-3),
    )),
]

Every Hail column not named in fields is dropped — locus/alleles become the key, and anything else absent from this list simply doesn't reach the file. What each Field does to its column:

Hail column	Field	on disk
filters	Field("filters")	copied through unchanged
af	Field("af", precision=1e-6)	float64 → int64, 0.0031 → 3100
cadd_raw	Field("cadd_raw", alias="cadd", precision=0.01)	renamed cadd, -1.44 → -144
vep	Field("vep", fields=(...))	stays array<struct<...>>, each member its own leaf

# One aggregation pass. Checks each precision against the range it will be
# applied to, and reports what it will do.
for report in check_precisions(ht, fields):
    print(report)
# {'column': 'af', 'multiplier': 1000000, 'min': 0.0001, 'max': 0.41, 'steps_used': 410000, ...}
# {'column': 'cadd', 'multiplier': 100, 'min': -1.44, 'max': 22.6, 'steps_used': 2260, ...}
# {'column': 'vep.polyphen', 'multiplier': 1000, 'min': 0.01, 'max': 1.0, 'steps_used': 1000, ...}

# Stage 1 (Hail + Spark, distributed): route short/long, encode, partition.
short, long = export_parquet(ht, "/tmp/staging", fields=fields)

# Stage 2 (pyarrow, one node): choose encodings, consolidate.
finalize_parquet(short, "/tmp/echtvar/short", fields=fields)
finalize_parquet(long, "/tmp/echtvar/long", fields=fields)

Reading it back — a key plus its chrom/chunk is a whole variant, quantized values divide back out by their multiplier, and the struct comes back as a list of dicts:

t = pq.read_table("/tmp/echtvar/short/chr1.parquet")
print(t.schema.names)
# ['key', 'filters', 'af', 'cadd', 'vep', 'chrom', 'chunk']

for i in range(t.num_rows):
    d = var32.decode(t["key"][i].as_py(), t["chrom"][i].as_py(), t["chunk"][i].as_py())
    print(f"{d.chrom}:{d.position + 1} {d.reference}>{d.alternate}")
    print(f"  filters={t['filters'][i].as_py()!r}  af={t['af'][i].as_py() / 1_000_000}"
          f"  cadd={t['cadd'][i].as_py() / 100}")
    print(f"  vep={t['vep'][i].as_py()}")

# chr1:100 A>T
#   filters='PASS'  af=0.0031  cadd=3.21
#   vep=[{'consequence': 'missense_variant', 'biotype': 'protein_coding', 'polyphen': 910},
#        {'consequence': 'intron_variant',   'biotype': 'lncRNA',         'polyphen': 50}]
# chr1:200 AC>A
#   filters='PASS'  af=0.41  cadd=-1.44
#   vep=[{'consequence': 'intron_variant', 'biotype': 'protein_coding', 'polyphen': 10}]

The third variant is missing from that output because A>ACGTACGT exceeds var32's four-base budget, so it routed to the long table instead:

t = pq.read_table("/tmp/echtvar/long/chr1.parquet")
print(t.schema.field("key").type, t.num_rows)
# list<element: uint32> 1

Why two stages

No single tool can write this file. Hail and Spark are needed for the encoding and the shuffle, at scale. But Spark writes Parquet through parquet-mr, which cannot be told a column's encoding at all — it infers one from a writer version and a dictionary setting. The choices that make the format worth having are unavailable there.

	stage 1 — export_parquet	stage 2 — finalize_parquet
engine	Hail + Spark	pyarrow
runs	distributed	one node, streaming a chunk at a time
does	routes short/long, encodes keys, applies value transforms, partitions, sorts	picks encodings, narrows the key to uint32, writes the footer config, consolidates
output	staging/short/chrom=chr1/chunk=0/…	out/short/chr1.parquet

Stage 2 output is one file per chromosome, one row group per chunk, so a chunk lookup is read_row_group(i) rather than opening a file. Measured over 400 chunks: same bytes to within 1%, and 2.9× faster on random chunk reads locally — the gap widens on object storage, where opening a file is a round trip.

API

Field

How one annotation is carried. Three shapes, and any of them may be array-valued — VEP-style annotations usually are, and an array stays an array rather than being flattened or joined into a delimited string.

Field("af")                              # carried as-is
Field("af", precision=1e-6)              # float -> int, to a stated accuracy
Field("consequence", categorical=True)   # dictionary-encoded (echtvar's string table)
Field("cadd_raw", alias="cadd", precision=0.01)   # renamed on the way out

Nested fields

The vep field above is the common case — nest fields to describe a struct's members, an array of them being one per transcript. Here's why it's shaped that way:

Field("csq", fields=(
    Field("consequence", categorical=True),
    Field("biotype", alias="tx_biotype", categorical=True),
    Field("polyphen", precision=1e-3),
))

Given [hl.struct(consequence="missense_variant", biotype="protein_coding", polyphen=0.912), ...] that writes one Parquet leaf per member:

csq.list.element.consequence   BYTE_ARRAY  RLE_DICTIONARY
csq.list.element.tx_biotype    BYTE_ARRAY  RLE_DICTIONARY
csq.list.element.polyphen      INT64       PLAIN            # 0.912 -> 912

Three things follow from each member being its own column:

• Each takes its own encoding. The two string members get their own dictionaries; the quantized one correctly gets none.

• fields is an allowlist. Members you don't name are dropped, which is how you take three columns out of a VEP annotation with forty.

• precision applies per member, and check_precisions reaches inside:

  check_precisions(ht, fields)
  # [{'column': 'csq.polyphen', 'multiplier': 1000, 'min': 0.05, 'max': 1.0, ...}]
  

The struct survives to disk as a struct — nothing is flattened or exploded, so a row reads back as a list of dicts:

t["csq"][0].as_py()
# [{'consequence': 'missense_variant', 'tx_biotype': 'protein_coding', 'polyphen': 912},
#  {'consequence': 'intron_variant',   'tx_biotype': 'lncRNA',         'polyphen': 50}]

The array is incidental. A plain struct nests the same way, to any depth — fields describes the shape, and whether it sits inside a list only changes the leaf paths:

Field("a", fields=(
    Field("b", fields=(Field("c", categorical=True), Field("score", precision=1e-3))),
    Field("n"),
))

a.b.c       BYTE_ARRAY  RLE_DICTIONARY
a.b.score   INT64       PLAIN            # 0.125 -> 125
a.n         INT32       PLAIN

# arrow type: struct<b: struct<c: string, score: int64>, n: int32>
# check_precisions -> {'column': 'a.b.score', 'multiplier': 1000, 'min': 0.125, 'max': 0.5}

precision is an accuracy budget, not a performance knob. Quantization shrinks a column by making values repeat, so what you keep is what you don't save. On 300k gnomAD-like allele frequencies, against 1.09MB as float32:

precision	size	vs float32	distinct values
1e-8	753KB	1.45×	112,963
1e-7	628KB	1.74×	58,461
1e-6	520KB	2.10×	12,943
1e-5	362KB	3.02×	1,338
1e-4	230KB	4.75×	135

An allele frequency is AC/AN, so it cannot mean anything below 1/AN — about 6.7e-7 at AN≈1.5M. Asking for finer spends bytes on the denominator wobbling between sites.

check_precisions(ht, fields) -> list[dict]

One aggregation pass. Raises if a precision cannot span its field's range without reaching the missing sentinel; otherwise reports the multiplier, observed min/max, and steps used. Advisory — nothing forces you to call it, but a precision that overflows silently turns real values into missing ones.

check_granularity(ht, fields, max_decimals=8) -> list[dict]

The companion to check_precisions, answering the question a range cannot. check_precisions sees only a min and a max, and [0.0, 1.0] looks identical whether the column holds hundredths or arbitrary float64s — so a precision far finer than the data's own step size passes it cleanly, then costs you twice: wider integers, and no new distinct values to show for them.

One aggregation pass. It reports how many decimal places each numeric field actually uses — a value's count is the smallest d where quantizing at 10**d and dequantizing returns it unchanged, and the field's answer is the widest any of its values needed:

for report in check_granularity(ht, fields):
    print(report["column"], report["decimals"], report["suggested"])
# af          2      0.01    # every value is hundredths; 1e-6 was four decades wasted
# cadd        2      0.01
# vep.revel   None   None    # continuous -- no decimal precision to find

suggested is just 10**-decimals, ready to paste into a Field. None means some value needed more than max_decimals — the honest answer for a continuous model score, and the signal to leave the field unquantized. exact_fraction says how near a miss that was: 0.999 is a handful of outliers, 0.0 is continuous throughout.

Runs whether or not a field has a precision yet — the point is to run it before choosing one. Non-numeric leaves are skipped.

The round trip is compared in the value's own type, which is what makes the answer independent of float width. A float32 cannot hold 12.34 — it holds 12.340000152587890625 — so asking how far v * 100 sits from a whole number measures float32's error rather than the data's, and would call every float32 field continuous. Divide back and compare as float32 and the round trip is exact, because the type's error cancels with itself. No tolerance to tune.

One caveat, mostly for small inputs: float32 has only ~7 significant decimal digits, so a handful of continuous float32 values can all round-trip within max_decimals and report a saturated number instead of None. At realistic row counts enough values need a 9th decimal that it reports None like float64 does. Either way, read an answer at or near max_decimals as "continuous" — the useful answers are the small ones.

It measures the granularity of the values as stored, which is not always the granularity they mean: AC/AN is a full-width float64 quotient and will look continuous here even though it cannot carry information below 1/AN. Take the coarser of this answer and what you know about the field.

export_parquet(ht, base_path, *, locus="locus", alleles="alleles", fields=(), mode="overwrite")

Stage 1. Returns (short_path, long_path). Routes each variant on len(ref) + len(alt) > 4: short variants take a 32-bit var32 key, long ones a kmer16 array. Multi-allelic rows contribute only their first alternate — split them (hl.split_multi_hts) first.

finalize_parquet(staging_path, output_path, *, fields=(), filesystem=None)

Stage 2. Pass the same fields. On GCS, hand it a pyarrow.fs.GcsFileSystem — pyarrow does not use Hadoop's GCS connector.

Decoding

var32.decode(key, chrom, chunk) and kmer16.decode(key, chrom) return a PRA(chrom, chunk, position, reference, alternate), 0-based. A var32 key stores position truncated to 20 bits, so it needs its chunk to reconstruct an absolute coordinate — which is why chunk is a column.

Output layout

out/short/chr1.parquet      key: uint32          row group i == chunk i
out/long/chr1.parquet       key: list<uint32>

Columns are key, chrom, chunk, and one per field. Position and alleles are not stored — the key plus chunk reconstructs them.

chunk is the index into the row groups. chrom identifies rows once several files are read together (ds.dataset(root) concatenates them and does not surface filenames). Both are constant over their scope, so RLE reduces them to almost nothing.

Quantized columns are unreadable without their multiplier, so it is written to the Parquet footer under the fields key:

import json
md = pq.ParquetFile("/tmp/echtvar/short/chr1.parquet").metadata
json.loads(md.metadata[b"fields"])
# [{'name': 'filters', 'column': 'filters'},
#  {'name': 'af', 'column': 'af', 'multiplier': 1000000, 'missing': 4294967295},
#  {'name': 'cadd_raw', 'column': 'cadd', 'multiplier': 100, 'missing': 4294967295},
#  {'name': 'vep', 'column': 'vep', 'fields': [
#      {'name': 'consequence', 'column': 'consequence'},
#      {'name': 'biotype', 'column': 'biotype'},
#      {'name': 'polyphen', 'column': 'polyphen', 'multiplier': 1000, 'missing': 4294967295}]}]

A struct field's own entry has no multiplier — it's a container, so the multiplier sits on the member that needs one, and the nesting in the config is mirrored in the footer.

Dictionary encoding is deliberately not recorded there — Parquet already reports each column chunk's encoding, and a footer note would be the less reliable of the two.

Development

uv sync            # dev tools are a PEP 735 dependency group
uv run pytest
uv run ruff check

The suite starts a real Hail/Spark session; it needs a JVM.