Applied - Reading VCF with pandas

Core Principle: Parsing Variant Call Format (VCF) files into pandas DataFrames requires parsing metadata headers, unpacking multi-allelic records, and restructuring INFO and FORMAT fields into tidy tabular formats.

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

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 132, in run_cell_capture
    exec(compile(code, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 15, in <module>
NameError: name 'vcf' is not defined

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

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 128, in run_cell_capture
    value = eval(compile(tail, "<cell>", "eval"), namespace)
  File "<cell>", line 1, in <module>
NameError: name 'vcf' is not defined

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

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 132, in run_cell_capture
    exec(compile(code, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 1, in <module>
NameError: name 'vcf' is not defined

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

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 126, in run_cell_capture
    exec(compile(head, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 3, in <module>
NameError: name 'vcf' is not defined

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

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 132, in run_cell_capture
    exec(compile(code, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 4, in <module>
NameError: name 'vcf' is not defined

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.

See also

Applied - Genomic File Formats · Applied - Annotating a Variant Table · String Accessor · Reshaping with pivot and melt · Missing Data in pandas · Tidy Data · hgvs Normalizer

Exercise

Untangle both levels of nesting the lesson describes: the INFO column and the per-sample FORMAT columns, with multi-allelic sites split into one row per alternate allele.

Write parse_vcf(text), returning a list of output lines. Skip lines beginning with ##. The #CHROM line is the header; if there is none, print error: no header line and nothing else. Rows with a different field count from the header are malformed and must be skipped. Split each record's ALT column on commas into one row per alternate allele. INFO is semicolon-separated key=value pairs, and a flag without a value counts as present with the value 1. AF is declared Number=A, so it carries one value per alternate allele: take the entry matching the allele's index, and print NA when AF is absent or has too few entries. DP is a single value for the whole record; print NA when it is missing.

Take the genotype from the first sample column, which is the tenth field onward, by finding GT's position in that record's FORMAT string; print NA when the file has no sample columns and ./. when the sample column is missing that subfield. Print records=<n> alleles=<n> pass=<n> where the pass count is the alleles whose FILTER is exactly PASS. Then print one line per passing allele in file order holding chrom:pos:ref:alt, the depth, the allele frequency as round(value, 3) renders it -- so 0.2 stays 0.2 -- and the genotype, tab separated. If no allele passes, print no variants after the count line.