Applied - Reading VCF with pandas

In one line: VCF is a tab file with two levels of nesting hidden inside it — the INFO column and the per-sample FORMAT columns — and untangling them is the exercise.

Uses: Reading and Writing Data · String Accessor · Reshaping with pivot and melt · pandas dtypes · Categorical dtype · Tidy Data

When pandas is the wrong tool

For millions of variants, or when you need region queries, use cyvcf2 or pysam — they parse in C and understand the format properly. Use pandas when the file is small enough to fit in memory and you want to join, group and plot it.

from cyvcf2 import VCF
for v in VCF("cohort.vcf.gz")("chr17:7668000-7688000"):
    v.CHROM, v.POS, v.REF, v.ALT, v.INFO.get("AF"), v.gt_types

The rest of this note assumes pandas is the right call.

Reading the fixed columns

import gzip, pandas as pd

def read_vcf(path):
    opener = gzip.open if str(path).endswith(".gz") else open
    with opener(path, "rt") as fh:
        for line in fh:
            if line.startswith("#CHROM"):
                cols = line.lstrip("#").rstrip("\n").split("\t")
                break
        else:
            raise ValueError("no #CHROM header found")

    return pd.read_csv(
        path, sep="\t", comment="#", names=cols,
        dtype={"CHROM": "category", "POS": "int32", "REF": "category",
               "ALT": "str", "FILTER": "category"},
        na_values=["."],
    )

vcf = read_vcf("cohort.vcf.gz")
vcf.info()

The header scan is necessary because the sample columns are not known in advance. comment="#" then skips all the ## meta lines when the data is read.

Untangling INFO

def parse_info(series):
    """'AF=0.3;DP=44;DB' → a DataFrame of columns."""
    parts = (series.fillna("")
                   .str.split(";")
                   .explode()
                   .str.split("=", n=1, expand=True))
    parts.columns = ["key", "value"]
    parts["value"] = parts["value"].fillna(True)      # flag fields have no '='
    parts["idx"] = parts.index
    return parts.pivot_table(index="idx", columns="key",
                             values="value", aggfunc="first")

info = parse_info(vcf["INFO"])
vcf = vcf.join(info)

# INFO values arrive as strings — cast the ones you need
vcf["AF"] = pd.to_numeric(vcf["AF"], errors="coerce").astype("float32")
vcf["DP"] = pd.to_numeric(vcf["DP"], errors="coerce").astype("Int32")

Flag fields (DB, SOMATIC) have no =, hence the fillna(True).

Faster, when you only want a few keys:

vcf["DP"] = vcf["INFO"].str.extract(r"(?:^|;)DP=(\d+)").astype("Int32")
vcf["AF"] = vcf["INFO"].str.extract(r"(?:^|;)AF=([\d.eE+-]+)").astype("float32")

The (?:^|;) anchor matters — without it, DP= also matches inside MQ_DP= or HaplotypeDP=. See String Accessor.

Untangling the sample columns

def parse_samples(vcf, sample_cols):
    long = vcf.melt(
        id_vars=["CHROM", "POS", "REF", "ALT", "FORMAT"],
        value_vars=sample_cols,
        var_name="sample", value_name="values",
    )
    keys = long["FORMAT"].str.split(":")
    vals = long["values"].str.split(":")
    long["fields"] = [dict(zip(k, v)) for k, v in zip(keys, vals)]
    return pd.concat([long.drop(columns=["FORMAT", "values", "fields"]),
                      pd.json_normalize(long["fields"])], axis=1)

sample_cols = vcf.columns[9:].tolist()
calls = parse_samples(vcf, sample_cols)

calls["DP"] = pd.to_numeric(calls["DP"], errors="coerce").astype("Int32")
calls[["ref_ad", "alt_ad"]] = (calls["AD"].str.split(",", expand=True)
                                          .apply(pd.to_numeric, errors="coerce"))
calls["vaf"] = calls["alt_ad"] / (calls["ref_ad"] + calls["alt_ad"])

Melting the sample columns is what makes the data tidy — one row per (variant, sample), which is the right observational unit for per-sample analysis.

Multi-allelic sites

vcf["ALT"].str.contains(",").sum()        # ← check this immediately

A multi-allelic record packs several ALT alleles into one row, and the AF, AC and AD fields then carry comma-separated values indexed to the ALT list. Most analysis code assumes one ALT per row and quietly gets it wrong.

Normalise upstream with bcftools rather than in pandas:

bcftools norm -m -any -f reference.fa input.vcf.gz -Oz -o normalized.vcf.gz

-m -any splits multi-allelics; -f reference.fa also left-aligns indels. Doing this before you read the file removes an entire class of bug. See hgvs Normalizer.

Filtering

passing = vcf.query("FILTER == 'PASS'")

rare = vcf.loc[
    (vcf["AF"].fillna(0) < 0.001)          # ← fillna! NA fails every comparison
    & (vcf["FILTER"] == "PASS")
    & (vcf["DP"] >= 20)
]

def report(df, mask, label):
    print(f"{label}: {mask.sum():,}/{len(df):,} ({mask.mean():.1%})")
    return df.loc[mask]

The fillna(0) is critical. A variant absent from gnomAD has NA frequency; NA fails < 0.001, so without the fill you silently exclude the rarest variants — the exact opposite of the intent. See Missing Data in pandas.

Summarising

vcf.groupby("CHROM", observed=True).size()
vcf["FILTER"].value_counts(dropna=False)

vcf["variant_type"] = np.select(
    [vcf["REF"].str.len().eq(1) & vcf["ALT"].str.len().eq(1),
     vcf["REF"].str.len() > vcf["ALT"].str.len(),
     vcf["REF"].str.len() < vcf["ALT"].str.len()],
    ["SNV", "deletion", "insertion"], default="complex",
)
vcf["variant_type"].value_counts()

# transition/transversion ratio — a standard QC metric
TRANSITIONS = {("A","G"), ("G","A"), ("C","T"), ("T","C")}
snv = vcf.query("variant_type == 'SNV'")
ti = snv.apply(lambda r: (r["REF"], r["ALT"]) in TRANSITIONS, axis=1)
print(f"Ti/Tv = {ti.sum() / (~ti).sum():.2f}")

# per-sample call rate
calls.groupby("sample")["GT"].apply(lambda s: (s != "./.").mean())

Ti/Tv is the single most informative VCF QC number. Whole-genome should be ~2.0–2.1, whole-exome ~2.8–3.0. Substantially lower means false positives; the metric is sensitive because sequencing errors are roughly random across substitution types while real variation is not.

Plotting

import seaborn as sns
sns.set_theme(style="ticks", context="paper")

CHROM_ORDER = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY", "chrM"]
vcf["CHROM"] = pd.Categorical(vcf["CHROM"], categories=CHROM_ORDER, ordered=True)

fig, axes = plt.subplots(2, 2, figsize=(10, 7), layout="constrained")

sns.countplot(data=vcf, x="CHROM", ax=axes[0,0])
axes[0,0].tick_params(axis="x", rotation=90, labelsize=6)

sns.histplot(data=vcf, x="AF", bins=50, log_scale=(True, False), ax=axes[0,1])
sns.histplot(data=vcf, x="DP", bins=50, ax=axes[1,0])
sns.countplot(data=vcf, x="variant_type", ax=axes[1,1])

The ordered categorical is what stops the x-axis reading chr1, chr10, chr11, chr2. See Categorical dtype.

Handing off to hgvs

from bioutils.assemblies import make_ac_name_map
AC = {v: k for k, v in make_ac_name_map("GRCh38").items()}

vcf["hgvs_g"] = vcf.apply(
    lambda r: f"{AC[r['CHROM'].removeprefix('chr')]}:g.{r['POS']}{r['REF']}>{r['ALT']}",
    axis=1,
)

Only correct for SNVs — indels need proper construction, not string formatting, because VCF and HGVS represent them differently. See Applied - Annotating a Variant Table and hgvs Pitfalls.

Common mistakes

Ignoring multi-allelic sites.
NA allele frequencies failing filters, dropping the rarest variants.
Unanchored INFO regexes matching the wrong key.
Forgetting INFO/FORMAT values are strings.
chr1 vs 1 breaking joins.
Assuming pandas can handle a whole-genome VCF. Use cyvcf2.
Not checking Ti/Tv.
Alphabetical chromosome order.
Building HGVS strings for indels by concatenation.

Applied - Reading VCF with pandas

Applied - Reading VCF with pandas

When pandas is the wrong tool

Reading the fixed columns

Untangling INFO

Untangling the sample columns

Multi-allelic sites

Filtering

Summarising

Plotting

Handing off to hgvs

Common mistakes

See also