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, 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", scale=6),  # quantized: to a stated accuracy
    Field("cadd_raw", alias="cadd", scale=2),  # quantized + renamed
    Field(
        "vep",
        fields=(  # struct: member-by-member
            Field("consequence", categorical=True),
            Field("biotype", categorical=True),
            Field("polyphen", scale=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:

# 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 arrive as Decimals at their own scale, 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()}"
        f"  cadd={t['cadd'][i].as_py()}"
    )
    print(f"  vep={t['vep'][i].as_py()}")

# chr1:100 A>T
#   filters='PASS'  af=0.003100  cadd=3.21
#   vep=[{'consequence': 'missense_variant', 'biotype': 'protein_coding', 'polyphen': Decimal('0.910')},
#        {'consequence': 'intron_variant',   'biotype': 'lncRNA',         'polyphen': Decimal('0.050')}]
# chr1:200 AC>A
#   filters='PASS'  af=0.410000  cadd=-1.44
#   vep=[{'consequence': 'intron_variant', 'biotype': 'protein_coding', 'polyphen': Decimal('0.010')}]

# No division and nothing to look up: the values ARE Decimals, at the scale
# the schema records.

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 2 output is one file per chromosome, several row groups per chunk, so a chunk lookup is a row-group read 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", scale=6)  # float -> DECIMAL, to a stated accuracy
Field("consequence", categorical=True)  # dictionary-encoded (echtvar's string table)
Field("cadd_raw", alias="cadd", scale=2)  # 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", scale=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.
• scale applies per member, and the DECIMAL cast reaches inside — csq.polyphen at scale=3 becomes decimal128(18, 3) while its sibling members are untouched.

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", scale=3))),
        Field("n"),
    ),
)

a.b.c       BYTE_ARRAY  RLE_DICTIONARY
a.b.score   INT64       PLAIN            # DECIMAL(18,3), reads back as 0.125
a.n         INT32       PLAIN

# arrow type: struct<b: struct<c: string, score: decimal128(18, 3)>, n: int32>
# the nested leaf a.b.score is quantized; the scale is in the schema

scale is an accuracy budget, not a performance knob — it is Parquet's DECIMAL scale, a count of decimal places, passed straight through. 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:

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.

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

Stage 1, and it applies no encoding — scale and categorical are both acted on by finalize_parquet, so staging carries values as they were, and a missing value is a Parquet null. There is no pre-flight aggregation over your data; Field rejects an impossible scale at construction, and the DECIMAL cast in stage 2 raises on a value that will not fit. 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          64k-row row groups, nested in chunks
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.

Row groups hold ~64,000 rows and nest inside chunks — a chunk becomes several of them, and none ever spans a chunk boundary. A lookup is two stages, and the order matters:

f = pq.ParquetFile("out/short/chr1.parquet")
md = f.metadata
ci, ki = md.schema.names.index("chunk"), md.schema.names.index("key")
want, key = chunk_of(pos0), var32.encode(pos0, ref, alt)


def holds(i):
    rg = md.row_group(i)
    k = rg.column(ki).statistics
    return rg.column(ci).statistics.min == want and k.min <= key <= k.max


rg = next(i for i in range(md.num_row_groups) if holds(i))  # footer only, no I/O
keys = f.read_row_group(rg, columns=["key"])["key"].to_pylist()
i = bisect.bisect_left(keys, key)

Stage 1 cannot be skipped. var32 keys overlap across chunks by construction (position is truncated to POSITION_BITS), so matching on key statistics alone lands in another chunk's row group. Within a chunk the rows stay key-ordered, so its row groups have disjoint ascending key ranges and exactly one can hold the key. And never index positionally — a chunk with no variants produces no row group at all.

Row-group size is what sets lookup cost for a pyarrow reader, which decodes a whole row group either way — it prunes on the page index for neither read_row_group() nor a filtered scan. That is a pyarrow limitation, not a property of the format: arrow-rs prunes to the page via RowSelection and the offset index, which is what _DATA_PAGE_SIZE exists for.

In the pyarrow regime the cost is linear in row-group size, measured on chr22:

The reader is now DataFusion, and in that regime row-group size does not set lookup cost at all. Pages are bounded by bytes (_DATA_PAGE_SIZE, 128 KiB), not rows, so a 4× wider row group yields 4× more pages of the same size — and arrow-rs reads about one page per projected column regardless. Measured on chr22, all eight families, 326 leaves, SELECT *, in parallel, under three independent samples (even stride and two random seeds):

Flat on latency, 0.26× on footer, and no compression penalty. The "below 16,000 compression collapses" finding is real but applies only going down; it is really a page effect — zstd's window resets per page, and tiny row groups force tiny pages — which is why it reappears at constant row-group size when _DATA_PAGE_SIZE is lowered: 128 KiB costs 2.4%, 32 KiB 11%, 8 KiB 30%.

Scan time alone does rise ~18% at 64,000, because building the RowSelection decodes 4× more pos and pos is DELTA_BINARY_PACKED. It is invisible end to end: once the footer is cached, scan is a minority of a lookup (0.2 ms of 21.6 ms serial across the eight families is metadata; the rest is scan).

What this buys is metadata. The full export measured 2.883 GB of footer over 470.1 GB of data at 16,000 — histograms alone is 1153.5 MB, and chr1 is 234.5 MB across all eight families, which is what one process holds resident to serve that chromosome. At 64,000 those become ~0.74 GB and ~60 MB.

An earlier revision of this section reported the opposite — that with page pruning 64,000 was ~2× slower than 16,000 (consequences 8.8 ms vs 4.5 ms). That measurement was taken on the keyed schema with a different reader, and has not been reproduced on the natural schema under DataFusion 54, where all three sizes are within noise. If the keyed path is ever revived, re-measure it there rather than assuming either result carries over.

Two caveats survive the change. _ROW_GROUP_ROWS sets row-group boundaries, so changing it re-breaks cross-family alignment unless every family is re-finalized together — exactly what the 0.7.0 → 0.8.0 move did. And any consumer still reading with pyarrow takes the 3.4× above.

One more, for anyone reading keys directly: do not read the key column with to_pylist(). It builds one Python int per row under the GIL — 5.6× on a single lookup, and it serialized concurrent lookups across files (10 threads gave 1.15×). Read .to_numpy(), which is zero-copy for uint32, and keep the row in Arrow; that gave 4.8× from concurrency instead.

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.

The files carry no metadata of ours. A quantized column is a Parquet DECIMAL, so its scale is in the schema; a categorical one announces itself through its column-chunk encoding; a missing value is a Parquet null. Nothing has to know what wrote the file to read it correctly:

schema = pq.ParquetFile("/tmp/echtvar/short/chr1.parquet").schema_arrow
schema.field("af").type  # decimal128(18, 6)  -> reads back as 0.010000
schema.field("cadd").type  # list<decimal128(18, 2)>

store_decimal_as_integer keeps the physical column an INT64, so being self-describing costs nothing on disk.

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.