Polars DataFrames¶

vepyr.annotate() without output_vcf returns a polars.LazyFrame. Nothing runs until you collect() or sink_*() it, and every collect starts a fresh annotation stream, so the same LazyFrame can be evaluated more than once.

import polars as pl
import vepyr

lf = vepyr.annotate(
    "input.vcf.gz",
    "/data/vepyr_cache/116_GRCh38_merged",
    reference_fasta="GRCh38.fa",
)
df = lf.collect()

Every column is the --everything value

Annotation always runs Ensembl VEP --everything, so reference_fasta is required and every column holds its --everything value. A select() only decides how much work that takes, see below.

Schema¶

collect_schema() is free: it comes from a probe that opens the cache but annotates nothing. This is the schema for chr22 of HG002 on the release-116 Ensembl cache:

>>> lf.collect_schema()
Schema({
    'chrom': String,
    'start': UInt32,                     # VCF POS, 1-based
    'end': UInt32,                       # POS + len(REF) - 1
    'id': String,                        # '' when the VCF has '.'
    'ref': String,
    'alt': String,
    'qual': Float64,
    'filter': String,
    # the input's own INFO and FORMAT fields, named by their VCF ids
    'DPSum': Int32, 'platforms': Int32, ..., 'GT': String, 'DP': Int32, ...,
    '_vcf_info_keys': String, '_vcf_format_keys': String,   # record layout, see below
    'most_severe_consequence': String,
    'Allele': String,
    'Consequence': List(String),         # one element per CSQ entry, in CSQ order ...
    'IMPACT': List(String),
    'SYMBOL': List(String),
    'Gene': List(String),
    'Feature_type': List(String),
    'Feature': List(String),
    'BIOTYPE': List(String),
    'EXON': List(String),
    'INTRON': List(String),
    'HGVSc': List(String),
    'HGVSp': List(String),
    'cDNA_position': List(String),
    'CDS_position': List(String),
    'Protein_position': List(String),
    'Amino_acids': List(String),
    'Codons': List(String),
    'Existing_variation': List(String),  # distinct values per variant, not per entry
    'DISTANCE': List(Int64),
    'STRAND': List(Int8),
    'FLAGS': List(String),
    'VARIANT_CLASS': String,
    'SYMBOL_SOURCE': List(String),
    'HGNC_ID': List(String),
    'CANONICAL': List(String),
    'MANE': List(String),
    'MANE_SELECT': List(String),
    'MANE_PLUS_CLINICAL': List(String),
    'TSL': List(Int8),
    'APPRIS': List(String),
    'CCDS': List(String),
    'ENSP': List(String),
    'SWISSPROT': List(String),
    'TREMBL': List(String),
    'UNIPARC': List(String),
    'UNIPROT_ISOFORM': List(String),
    'GENE_PHENO': List(String),
    'SIFT': List(String),
    'PolyPhen': List(String),
    'DOMAINS': List(String),
    'miRNA': List(String),
    'HGVS_OFFSET': List(Int64),
    'AF': Float32,                       # frequencies are per variant ...
    'AFR_AF': Float32,
    'AMR_AF': Float32,
    'EAS_AF': Float32,
    'EUR_AF': Float32,
    'SAS_AF': Float32,
    'gnomADe_AF': Float32,
    'gnomADe_AFR_AF': Float32,
    'gnomADe_AMR_AF': Float32,
    'gnomADe_ASJ_AF': Float32,
    'gnomADe_EAS_AF': Float32,
    'gnomADe_FIN_AF': Float32,
    'gnomADe_MID_AF': Float32,
    'gnomADe_NFE_AF': Float32,
    'gnomADe_REMAINING_AF': Float32,
    'gnomADe_SAS_AF': Float32,
    'gnomADg_AF': Float32,
    'gnomADg_AFR_AF': Float32,
    'gnomADg_AMI_AF': Float32,
    'gnomADg_AMR_AF': Float32,
    'gnomADg_ASJ_AF': Float32,
    'gnomADg_EAS_AF': Float32,
    'gnomADg_FIN_AF': Float32,
    'gnomADg_MID_AF': Float32,
    'gnomADg_NFE_AF': Float32,
    'gnomADg_REMAINING_AF': Float32,
    'gnomADg_SAS_AF': Float32,
    'MAX_AF': Float32,
    'MAX_AF_POPS': String,
    'CLIN_SIG': List(String),            # distinct values per variant
    'SOMATIC': String,
    'PHENO': String,
    'PUBMED': List(String),              # distinct values per variant
    'MOTIF_NAME': List(String),          # per entry, MotifFeature entries only
    'MOTIF_POS': List(Int64),
    'HIGH_INF_POS': List(String),
    'MOTIF_SCORE_CHANGE': List(Float32),
    'TRANSCRIPTION_FACTORS': List(String),
    'clin_sig_allele': List(String),     # cache-only columns, no CSQ counterpart
    'clinical_impact': String,
    'minor_allele': String,
    'minor_allele_freq': Float32,
    'clinvar_ids': List(String),
    'cosmic_ids': List(String),
    'dbsnp_ids': List(String),
})

The 80 columns from Allele to TRANSCRIPTION_FACTORS are the CSQ fields in VCF header order. A List column that is aligned with Consequence has one element per CSQ entry; Existing_variation, CLIN_SIG and PUBMED hold the distinct values for the variant instead. Values that VEP repeats on every entry, such as the frequencies, are stored once as scalars. The input's INFO fields and its sample columns are carried too, named by their VCF ids; pass info_fields=[] and format_fields=[] for the 0.7 frame. An id that is also an annotation column arrives as INFO_<id> (fmt_<id> for FORMAT): AF is always VEP's. Add skip_csq=False to get the raw CSQ string as a column.

CSQ is the one carried field this run can replace. Annotating an input that already has INFO/CSQ with skip_csq=False drops the input's, as Ensembl VEP and output_vcf do, so the frame carries this run's consequences and no INFO_CSQ. Without a CSQ of its own (skip_csq=True, the default) there is nothing to replace it with, and the input's is carried and written back untouched.

To narrow the frame, select() the columns you need, see below; the engine then only computes what those columns require.

Plugin columns¶

With a plugin cache configured, every plugin CSQ field is a named column appended after the block above; see Plugins. Its shape and type come from the plugin's manifest, not from the CSQ text it travels in:

A plugin keyed on the variant alone, such as CADD or ClinVar, gives one value per row: CADD_PHRED: String, ClinVar_CLNSIG: String.
A plugin keyed on a transcript feature, such as SpliceAI (by gene symbol) or AlphaMissense (by protein change), gives one value per consequence entry, aligned with Consequence: SpliceAI_pred_DP_AG: List(Int32), am_pathogenicity: List(Float32). Interval plugins such as PhenotypeOrthologous (gene span + Gene) are lists too, empty on entries whose transcript lies outside the matching span.
The element type is the manifest's type: Utf8, Int32 or Float32. CADD's scores are declared Utf8 so the VCF reproduces the source digits exactly, which is why CADD_PHRED is a string; cast it when you need a number.

See String scores and numeric casts in the vepyr-plugins README for per-field numeric types and Polars examples covering scalar columns, consequence-aligned lists and multiple scores inside a dbNSFP element.

The engine only runs the plugin lookup, and only builds the CSQ string the values are carried in, when a query reads one of those columns:

lf = vepyr.annotate(
    "input.vcf.gz", cache, reference_fasta="GRCh38.fa",
    plugin_cache_root="/data/plugin_cache", plugins=["cadd"],
)
lf.select("chrom", "start", "ref", "alt", "CADD_PHRED")   # plugin lookup, no HGVS
lf.select("chrom", "start", "SYMBOL")                     # neither

Selecting a plugin column on a frame without a plugin cache fails when the plan is built, with Polars' ColumnNotFoundError listing the columns the frame has.

What is pushed into the engine¶

The frame is backed by a Polars IO plugin that pulls Arrow batches from the native annotator. Five things reach the engine — the Polars pushdown diagram shows how they travel from the query to the annotator:

head(n) and limit(n) become a SQL LIMIT, so previewing is fast.
A filter() on chrom, start or end restricts the input before annotation; see Region filters below.
A narrowing select(), together with the columns a pushed-down filter() reads, decides how much of the work runs. The engine skips the HGVS lookup when neither HGVSc nor HGVSp is needed, and the co-located lookup when none of its columns (Existing_variation, CLIN_SIG, SOMATIC, PHENO, PUBMED, the AF family and the cache-only columns) is needed. Everything else --everything adds (SIFT, PolyPhen, MANE, DOMAINS, the motif columns, the per-population gnomAD columns, …) runs only when one of those columns is needed. The selected columns are value-identical to a plain collect().
The input columns a query names are the only ones parsed: select("chrom", "SYMBOL") reads no INFO and no samples.
The plugin columns a query names decide which plugins run. Only the plugins that own a column the query reads are looked up, and only their columns are parsed out of CSQ: select("chrom", "start", "CADD_PHRED") runs the CADD lookup and no other.

Every other filter() is applied to each batch after annotation. It bounds memory, because a batch is dropped as soon as it has been reduced, but it does not reduce the engine's work. The raw CSQ string (skip_csq=False) needs every lookup and every configured plugin, so a query that reads it runs like a plain collect(). A query that reads neither CSQ nor a plugin column runs no plugin lookup at all.

preview = lf.head(20).collect()              # LIMIT 20 in the engine

high = (
    lf.filter(pl.col("IMPACT").list.contains("HIGH"))
      .select("chrom", "start", "ref", "alt", "SYMBOL", "Consequence")
      .collect()
)                                            # runs without HGVS or the co-located lookup

lf.select("chrom", "start", "SYMBOL", "Consequence").collect()   # neither HGVS nor co-located
lf.select("chrom", "start", "HGVSc").collect()                    # HGVS, no co-located lookup
lf.select("chrom", "start", "AF", "CLIN_SIG").collect()           # co-located lookup, no HGVS
lf.select("chrom", "start", "SIFT").collect()                     # neither HGVS nor co-located
lf.collect()                                                      # no projection: all of --everything

On the release-116 Ensembl cache with a FASTA and workers=1, measured on an Apple Silicon M3 Max (16 cores, 64 GiB):

Input	Query	Wall time
chr22, 50,861 variants	`collect()`	2.3 s
	`select(chrom, start, ref, alt, SYMBOL, Consequence, IMPACT)`	1.3 s
	`select(chrom, start, HGVSc, HGVSp)`	1.4 s
	`select(chrom, start, Existing_variation, AF, MAX_AF, CLIN_SIG, PUBMED)`	1.7 s
chr1, 323,430 variants	`collect()`	13.5 s
	`select(chrom, start, ref, alt, SYMBOL, Consequence, IMPACT)`	6.1 s
	`select(chrom, start, HGVSc, HGVSp)`	6.7 s
	`select(chrom, start, Existing_variation, AF, MAX_AF, CLIN_SIG, PUBMED)`	11.0 s

Region filters¶

A filter() on chrom, start or end is pushed into the engine before annotation: contigs outside the filter are never prepared, and an indexed input (bgzip + .tbi/.csi) is read by seek.

df = lf.filter(
    (pl.col("chrom") == "chr22") & pl.col("start").is_between(20_000_000, 25_000_000)
).collect()

The result is always identical to filtering after collect(); only the work changes. Coordinates are the frame's own start/end columns (1-based, closed). Recognised shapes:

chrom conjuncts that name specific contigs: ==, is_in, str.starts_with and boolean combinations of them. A conjunct that would accept any name (!=, is_not_null()) narrows nothing on its own and is only pushed together with a range.
start/end conjuncts: comparisons with an integer literal and is_between. end <= b bounds the range; end >= a does not.
Several regions: an | of (chrom & range) groups, one region per group.

Anything else (a float literal, a range compared to another column, a cast, an | inside a range conjunct) is not pushed down and is applied by Polars after annotation, which is still correct, just not faster.

Without a tabix/CSI index next to the input a RuntimeWarning is emitted: the whole file is parsed once to find which contigs it carries and once more to filter it before annotation, and only the selected rows are annotated. On Merged and RefSeq caches a range costs one extra positional pass over each selected contig, which keeps the result byte-identical to a whole-file run.

On the release-116 caches with a FASTA and workers=1 (HG002 slices, indexed input), measured on an Apple Silicon M3 Max (16 cores, 64 GiB):

Input	Query	Ensembl	Merged	RefSeq
chr22, 50,861 variants	`collect()`	2.3 s	3.1 s	1.7 s
	`filter(chr22:20,000,000-25,000,000)`, 5,457 rows	0.6 s	1.2 s	0.7 s
	`filter(chr22:30,000,000-30,100,000)`, 60 rows	0.2 s	0.6 s	0.4 s
chr1, 323,430 variants	`collect()`	12.8 s	17.7 s	11.8 s
	`filter(chr1:20,000,000-25,000,000)`, 7,994 rows	1.1 s	2.3 s	1.4 s
	`filter(chr1:30,000,000-30,100,000)`, 278 rows	0.7 s	1.4 s	0.9 s

One row per consequence¶

Most downstream work, and every filter_vep expression, is about one CSQ entry: a variant paired with one transcript, regulatory or motif feature. Explode the aligned list columns together to get that long format. The helper pads a null list to the row's entry count first, so a column that is null on some rows cannot break the alignment.

UNALIGNED = {"Existing_variation", "CLIN_SIG", "PUBMED",
             "clin_sig_allele", "clinvar_ids", "cosmic_ids", "dbsnp_ids"}


def consequence_rows(df: pl.DataFrame) -> pl.DataFrame:
    """One row per CSQ entry. Keeps every scalar and per-variant column as is."""
    # The annotation columns start at "Allele"; the input's own INFO and FORMAT
    # fields come before it and may be lists of their own length (AD, ...).
    names = list(df.schema)
    cols = [
        c for c in names[names.index("Allele"):]
        if isinstance(df.schema[c], pl.List) and c not in UNALIGNED
    ]
    padded = [
        pl.when(pl.col(c).is_null())
        .then(pl.col("Consequence").list.eval(pl.lit(None)).cast(df.schema[c]))
        .otherwise(pl.col(c))
        .alias(c)
        for c in cols
    ]
    return (
        df.filter(pl.col("Consequence").is_not_null())
          .with_columns(padded)
          .explode(cols, empty_as_null=True)
    )


long = consequence_rows(df)

Rows whose ALT is * have no consequences and are dropped by the filter. On chr22 of HG002 the 50,861 variants become 943,901 consequence rows.

`filter_vep` expressions in Polars¶

Ensembl's filter_vep evaluates an expression per CSQ entry and keeps the variant when any entry matches; with --only_matched it also drops the entries that did not. Both map onto the frame:

On the long frame (long above), filter() keeps matching entries only. This is filter_vep --only_matched.
On the wide frame (df), wrap the per-entry test in list.eval(...).list.any() to keep the whole variant when any entry matches. This is plain filter_vep.

The table uses the long frame. Fields that filter_vep treats as numbers after stripping text, such as SIFT's score in parentheses or the exon number in 2/10, need the same extraction here.

`filter_vep` expression	Polars expression on `long`
`Feature is ENST00000307301`	`pl.col("Feature") == "ENST00000307301"`
`Feature_type is Transcript`	`pl.col("Feature_type") == "Transcript"`
`Consequence is missense_variant`	`pl.col("Consequence") == "missense_variant"`
`Consequence matches stream`	`pl.col("Consequence").str.contains("stream")`
`Consequence match stop`	`pl.col("Consequence").str.contains("stop")`
`SIFT != tolerated`	`pl.col("SIFT").str.replace(r"$.*$$", "") != "tolerated"`
`SIFT match tolerated`	`pl.col("SIFT").str.contains("tolerated")`
`SIFT < 0.1`	`pl.col("SIFT").str.extract(r"$([\d.]+)$").cast(pl.Float64) < 0.1`
`Protein_position < 10`	`pl.col("Protein_position").str.extract(r"^(\d+)").cast(pl.Int64) < 10`
`Exon > 1`	`pl.col("EXON").str.extract(r"^(\d+)").cast(pl.Int64) > 1`
`MANE` (field exists)	`pl.col("MANE").is_not_null()`
`SYMBOL`	`pl.col("SYMBOL").is_not_null()`
`not Existing_variation`	`pl.col("Existing_variation").is_null()`
`AF < 0.01 or not AF`	`pl.any_horizontal(pl.col("AF") < 0.01, pl.col("AF").is_null())`
`AFR_AF > 0.1 or EUR_AF > 0.1`	`pl.any_horizontal(pl.col("AFR_AF") > 0.1, pl.col("EUR_AF") > 0.1)`
`AFR_AF > #EUR_AF`	`pl.col("AFR_AF") > pl.col("EUR_AF")`
`(AFR_AF > 0.1 or EUR_AF > 0.1) and (EAS_AF < 0.1 and SAS_AF < 0.1)`	`pl.any_horizontal(pl.col("AFR_AF") > 0.1, pl.col("EUR_AF") > 0.1) & pl.all_horizontal(pl.col("EAS_AF") < 0.1, pl.col("SAS_AF") < 0.1)`
`Consequence is missense_variant and CCDS and DOMAINS`	`(pl.col("Consequence") == "missense_variant") & pl.col("CCDS").is_not_null() & pl.col("DOMAINS").is_not_null()`
`SYMBOL in BRCA1,BRCA2`	`pl.col("SYMBOL").is_in(["BRCA1", "BRCA2"])`
`Feature in /data/files/motifs_list.txt`	`pl.col("Feature").is_in(Path("/data/files/motifs_list.txt").read_text().split())`
`Consequence is coding_sequence_variant` with `--ontology`	No ontology expansion; list the child terms: `pl.col("Consequence").str.split("&").list.eval(pl.element().is_in(CHILD_TERMS)).list.any()`

Operator by operator:

`filter_vep`	Polars
`is`, `=`, `eq`	`==` on the whole value, including `&`-joined fields
`!=`, `ne`	`!=`
`match`, `matches`, `re`, `regex`	`str.contains("(?i)" + pattern)`: `filter_vep` matches case-insensitively
`<`, `>`, `<=`, `>=`	the same, after `cast()` from the string form where needed
`exists`, `ex`, `defined`, bare field	`is_not_null()`
`not`	`~` or `is_null()`
`and`, `or`, parentheses	`&` and the pipe operator, each side in its own parentheses; or `pl.all_horizontal()` / `pl.any_horizontal()`
`in a,b,c` / `in file`	`is_in([...])`

Where a value in the frame is null, the filter_vep field is empty, so is_null() and is_not_null() give the not FIELD and FIELD semantics exactly.

is compares the whole value, also for &-joined fields: an entry whose Consequence is missense_variant&splice_region_variant does not satisfy Consequence is missense_variant, so the Polars form is ==, not a test on the split terms. To find a term inside a joined value, use match, which is str.contains. On HG002 chr22 with --only_matched, Consequence is missense_variant and SIFT < 0.05 keeps 789 entries in 83 variants and Consequence match missense_variant and SIFT < 0.05 keeps 983 entries in 108 variants, in both filter_vep 116.0 and Polars.

filter_vep's parser splits the filter string on spaces and parentheses and drops ^ from every word, so a pattern can contain none of them: Consequence regex ^splice is evaluated as Consequence regex splice.

The same expressions on the wide frame keep whole variants:

any_stream = df.filter(
    pl.col("Consequence").list.eval(pl.element().str.contains("stream")).list.any()
)
any_high = df.filter(pl.col("IMPACT").list.contains("HIGH"))

Plugins¶

Point the frame at a plugin cache and the plugin fields are columns like any other; see Plugin columns for how their shape and type are decided. All examples below use the four published caches:

lf = vepyr.annotate(
    "input.vcf.gz", cache, reference_fasta="GRCh38.fa",
    plugin_cache_root="/data/plugin_cache",
    plugins=["cadd", "clinvar", "spliceai", "alphamissense"],
)

CADD is per variant, and its scores are strings so the VCF keeps the source digits. Cast, then filter:

high_cadd = (
    lf.select("chrom", "start", "ref", "alt", "SYMBOL", "CADD_PHRED")
      .with_columns(pl.col("CADD_PHRED").cast(pl.Float32))
      .filter(pl.col("CADD_PHRED") >= 20)
      .collect()
)

ClinVar is per variant too, so its fields are plain string columns. ClinVar holds the variation id, and the assertion fields keep VEP's & joins, so match with a pattern:

clinvar = (
    lf.select("chrom", "start", "ref", "alt", "ClinVar", "ClinVar_CLNSIG", "ClinVar_CLNREVSTAT")
      .filter(pl.col("ClinVar_CLNSIG").str.contains("(?i)pathogenic"))
      .collect()
)

SpliceAI is matched by gene symbol, so every field is a list aligned with Consequence. The delta scores are strings in the cache; cast the lists and take the largest score across the variant's entries:

DS = ["SpliceAI_pred_DS_AG", "SpliceAI_pred_DS_AL", "SpliceAI_pred_DS_DG", "SpliceAI_pred_DS_DL"]

splice = (
    lf.select("chrom", "start", "ref", "alt", "SpliceAI_pred_SYMBOL", *DS)
      .with_columns(pl.col(DS).cast(pl.List(pl.Float32)))
      .with_columns(pl.max_horizontal(pl.col(DS).list.max()).alias("spliceai_max_ds"))
      .filter(pl.col("spliceai_max_ds") >= 0.5)
      .collect()
)

AlphaMissense is matched by protein change, so am_class and am_pathogenicity (already List(Float32)) belong to the consequence entry that produced the change. Explode them together with the transcript columns to see which transcript each score refers to:

am = (
    lf.select("chrom", "start", "ref", "alt", "Feature", "Protein_position",
              "Amino_acids", "am_class", "am_pathogenicity")
      .filter(pl.col("am_class").list.contains("likely_pathogenic"))
      .explode(["Feature", "Protein_position", "Amino_acids", "am_class", "am_pathogenicity"])
      .filter(pl.col("am_class") == "likely_pathogenic")
      .collect()
)

Plugin columns combine freely with the base annotation, and the projection still decides what runs. This query runs the co-located lookup for gnomADg_AF and the CADD lookup, and nothing else:

candidates = (
    lf.select("chrom", "start", "ref", "alt", "SYMBOL", "Consequence", "gnomADg_AF", "CADD_PHRED")
      .with_columns(pl.col("CADD_PHRED").cast(pl.Float32))
      .filter(
          pl.any_horizontal(pl.col("gnomADg_AF").is_null(), pl.col("gnomADg_AF") < 0.001)
          & (pl.col("CADD_PHRED") >= 25)
      )
      .collect()
)

The consequence_rows helper works unchanged on a frame with plugins: per-feature plugin lists explode alongside the transcript columns, per-variant scalars such as CADD_PHRED repeat on each row, like the frequencies do.

A query that reads no plugin column skips the plugin lookup entirely, and one that reads plugin columns runs only the plugins that own them. On chr22 with the four caches, select(chrom, start, SYMBOL) takes 1.3 s and select(chrom, start, CADD_PHRED), which runs the CADD lookup alone, 4.2 s, so keep plugin columns out of queries that do not need them.

Workers¶

workers=N splits each contig into grid-aligned runs that are annotated concurrently and released in order. It applies to the LazyFrame path as well as output_vcf, and the frame is identical to workers=1, row order included.

lf = vepyr.annotate(
    "input.vcf.gz",
    "/data/vepyr_cache/116_GRCh38_ensembl",
    reference_fasta="GRCh38.fa",
    workers=8,
)
df = lf.collect()

workers > 1 needs an indexed input

The input VCF must be bgzip-compressed with a .tbi or .csi index. An unindexed input raises rather than silently falling back.

See Performance for a sweep across worker counts and cache profiles.

Writing results to disk¶

collect() holds the whole result in memory. On chromosome 1 of a whole-genome sample (323,430 variants) the collected frame is about 1.2 GB, but the process peaks at 8.4 GB because Polars buffers the Arrow batches on top of the frame. Two things bring that down.

Leave the CSQ string off. It is off by default (skip_csq=True); with it on, the chr1 frame more than doubles to 2.7 GB and collect() peaks about 2 GB higher. Turn it on only when you need the exact VEP string.

Stream with sink_parquet and a small row group. Polars holds one row group in memory before writing it, and its default group is far larger than an annotation batch. Setting the group to the engine's buffer size keeps the stream flat:

lf.sink_parquet("annotated.parquet", row_group_size=5000)

Measured on the release-116 Ensembl cache with workers=1 (median of three runs after a warm-up, one process per run):

Input	Path	Wall time	Peak RSS
chr22, 50,861 variants	`output_vcf` (bgzf)	3.5 s	1.1 GB
	`collect()`	2.4 s	1.9 GB
	`sink_parquet()` default	2.4 s	2.5 GB
	`sink_parquet(row_group_size=5000)`	2.3 s	1.5 GB
chr1, 323,430 variants	`output_vcf` (bgzf)	17.9 s	3.3 GB
	`collect()`	12.6 s	8.4 GB
	`sink_parquet()` default	12.6 s	6.6 GB
	`sink_parquet(row_group_size=5000)`	12.4 s	3.8 GB

With skip_csq=False add roughly 1 GB on chr1 for the small-row-group sink and 2 GB for collect(), and about 1.5 s of wall time to either. The engine itself needs about 3 GB on chr1 whichever output you choose, so the small-row-group sink is within half a gigabyte of the VCF writer. The small row group costs no wall time; the Parquet file grows by about 7 %.

List columns cannot be written to CSV directly. Join them with & first, on the wide frame or on the long frame, where the per-variant value sets are the only lists left:

def join_lists(df: pl.DataFrame) -> pl.DataFrame:
    return df.with_columns(pl.col(pl.List).cast(pl.List(pl.String)).list.join("&"))


join_lists(df).write_csv("annotated.tsv", separator="\t")
join_lists(consequence_rows(df)).write_csv("annotated_long.tsv", separator="\t")

Writing a filtered VCF¶

Filter the frame in Polars and write it back as VCF with polars-bio's sink_vcf:

pip install "vepyr[polars-bio]"    # polars-bio >= 0.36.1, Python 3.11+

import polars as pl
import polars_bio as pb
import vepyr

lf = vepyr.annotate(
    "HG002.vcf.gz",
    cache_dir,
    reference_fasta="GRCh38.fa",
    skip_csq=False,           # the CSQ column is what gets written
)
rare_damaging = lf.filter(
    pl.col("IMPACT").list.contains("HIGH")
    & (pl.col("MAX_AF").is_null() | (pl.col("MAX_AF") < 0.01))
    & (pl.col("DP") >= 20)
)
pb.sink_vcf(rare_damaging, "rare_damaging.vcf")

The file is what output_vcf would have written, minus the rows you filtered out. On HG002 chr22 against the release-116 cache the unfiltered body has the same md5 as output_vcf's, which passes the strict md5 gate against Ensembl VEP, and every line of a filtered file is byte-identical to its output_vcf line. The header is the input's own, line for line, followed by vepyr's provenance and the CSQ declaration; its command-line line records the annotation but has no output or compression, and it does not record what you did to the frame afterwards.

That relies on three things:

skip_csq=False. Without the CSQ column there is nothing to write.
The record layout. _vcf_info_keys and _vcf_format_keys hold each record's own INFO key order and FORMAT key list, which the typed columns cannot: every record carries every declared key, and a . value parses to the same null as an absent key. They are carried by default and have to stay in the frame; a select() that drops them writes keys in header order and leaves out a FORMAT key that is missing in every sample. Pass preserve_record_layout=False to leave them out altogether. The default does without them for BCF input and for a file that declares a field with either name; asking for them with True raises there instead.
Canonical values. Values go through typed columns, so QUAL 50.0 is written 50 and AF=0.50 as 0.5 — as output_vcf does too.

Agreement with the VCF output¶

The DataFrame and VCF paths run the same engine. On chr22 and chr1 of HG002 against the release-116 Ensembl cache, the CSQ column (skip_csq=False) matched the VCF's INFO/CSQ byte for byte on every record, the variant columns matched CHROM, POS, ID, REF, ALT, QUAL and FILTER, and every typed column matched its CSQ field element for element, motif entries included. The only representational difference is the three per-variant value-set columns:

Per-variant value sets

Existing_variation, CLIN_SIG and PUBMED are deduplicated per variant rather than repeated per consequence. No values are lost, but their list length does not match the per-consequence columns. Elements are split on the cache's , separator only, so a combined assertion such as benign&likely_benign stays one element; match with str.contains.

When a downstream step needs the string itself, decode it directly:

CSQ_FIELDS = [...]  # the Format: list from the VCF header, or the schema names from "Allele" on

entries = pl.col("CSQ").str.split(",")
decoded = df.with_columns(
    entries.list.eval(pl.element().str.split("|").list.get(i, null_on_oob=True))
           .alias(name)
    for i, name in enumerate(CSQ_FIELDS)
)