Bio Vcf Basics

Version Compatibility

Reference examples tested with: bcftools 1.19+, cyvcf2 0.30+, numpy 1.26+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

VCF/BCF Basics

"Show me and extract fields from this VCF" -> Parse the VCF/BCF format, then view, subset, or pull specific columns into a flat table.

• CLI: bcftools view / bcftools query -f
• Python: cyvcf2.VCF (iterate records with attribute access)

The governing principle

A VCF field is only meaningful once its LEVEL and its Number are known. QUAL is a site-level property; GQ, PL, AD, DP, GT are per-sample. QUAL and GQ answer different questions and are NOT interchangeable. A field's header Number (A/R/G/.) dictates how many values it carries and how it must be re-subset after a multiallelic split. And several encodings are load-bearing traps: . (missing) is never 0/0 (hom-ref); a bare * ALT is a spanning-deletion placeholder, not an allele; a gVCF <NON_REF> record is a reference-confidence intermediate, not a filtered call. Read the header, read the Number, read the level -- a structurally valid VCF read at the wrong level silently produces wrong numbers with no error.

Format Overview

VCF Format Structure

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE1
chr1    1000    rs123   A       G       30      PASS    DP=50   GT:DP   0/1:25

Header Lines (##)

• ##fileformat - VCF version
• ##INFO / ##FORMAT - INFO / FORMAT field definitions (ID, Number, Type)
• ##FILTER - Filter definitions
• ##contig - Reference contigs (required for region indexing and contig order)
• ##reference - Reference genome

The Header Contract

Every INFO/FORMAT tag used in the body MUST be declared in a ##INFO/##FORMAT line giving its ID, Number, and Type; parsers (bcftools, cyvcf2, pysam) read these declarations to know how many values a field holds and how to type it. An out-of-sync header -- a tag used but not declared, or declared with the wrong Number/Type -- silently breaks parsing: a Number=1 declaration over data that holds a vector, or a missing ##contig, makes tools drop, mistype, or mis-subset values with NO error thrown. After any hand-edit or annotation that adds a field, update the header to match (bcftools +fill-tags and bcftools annotate manage this automatically).

Data Columns

Critical Field Interpretation

What each field actually measures -- and what it does not -- drives every filtering and interpretation decision. QUAL, GQ, and PL answer three DIFFERENT questions.

QUAL vs GQ vs PL/GL: three different confidences

QUAL is computed once across all samples and SCALES with total depth, so a high-coverage artifact can carry a large QUAL -- hence QD (QUAL normalized by depth) is preferred for filtering. GQ is per-sample and does not scale with cohort size. They are NOT interchangeable: QUAL can be high while an individual genotype is uncertain (low GQ), and a sample can have a confident genotype (high GQ) at a site with only moderate QUAL. Filter site-level junk on QUAL/QD; no-call untrustworthy genotypes on GQ.

PL/GL, and how GQ is derived

PL holds phred-scaled genotype likelihoods, rebased so the CALLED (most likely) genotype is exactly 0 and every other value is its phred penalty relative to that call. For a biallelic diploid site PL is ordered [PL(0/0), PL(0/1), PL(1/1)] -- the index of the 0 IS the genotype the caller assigned. GL is the same information as raw log10 likelihoods (<=0, larger is better). GQ = the difference between the two SMALLEST PL values, i.e. phred confidence in the call versus the next-best genotype; GQ=0 means the top two genotypes are tied (uninformative), GQ is capped at 99 by convention.

For a site with n alleles, diploid genotype j/k (j<=k) sits at PL index k*(k+1)/2 + j (this is the Number=G ordering). Getting this index formula wrong is the classic bug when re-parsing PL after a multiallelic split -- the vector must be re-subset by the formula, never sliced positionally.

AD vs DP, and allele balance

AD (FORMAT, Number=R) is per-allele read depth [ref_depth, alt1_depth, ...], REF first. DP is total depth. sum(AD) is often LESS than DP -- expected, not an error:

• DP counts all reads spanning the position, including uninformative reads (low base quality, ambiguous alignment, filtered reads).
• AD counts only reads that confidently support a specific allele.
• INFO/DP (site-level, summed across samples) differs from FORMAT/DP (per-sample).

Allele balance for a het is DERIVED from AD (GATK does not emit it directly): AB = alt_AD / (ref_AD + alt_AD). A true het sits near 0.5; hets far from 0.5 (e.g. <0.2 or >0.8) suggest a mapping artifact, CNV, or contamination.

INFO/FORMAT Number semantics (A / R / G / .)

Every ##INFO/##FORMAT header declares a Number telling a parser how many values a field holds AND how to re-subset it when a multiallelic record is split:

Load-bearing for correctness: bcftools norm -m- uses these codes to reapportion fields on split. A field mis-declared Number=. when it is really A keeps its full multiallelic vector on every split record, so downstream tools read the WRONG allele's value with no error. See variant-calling/variant-normalization for split/join reapportionment.

Key INFO Annotations for Filtering

Genotype Encoding

Phased vs Unphased

• / separates unphased alleles -- the two chromosomal copies are known, but which came from which parent is not
• | separates phased alleles -- haplotype assignment is known (read-backed phasing, trio analysis, or long-read sequencing)
• Phasing matters for compound heterozygosity: two variants in a gene are pathogenic together only if on different haplotypes (in trans), not the same haplotype (in cis)

Phase Sets (PS)

A | is only meaningful WITHIN a phase set. The FORMAT/PS tag (an integer, usually the POS of the block's first variant) groups variants phased relative to EACH OTHER; 0|1 in two different PS blocks are not guaranteed to lie on the same physical haplotype. Read-backed phasers (WhatsHap) and trio phasing emit PS. A | with no consistent PS across records carries no global phase -- a subtle trap when merging phased VCFs.

Ploidy and Missing vs Hom-Ref

Ploidy is read from the NUMBER of alleles in GT: 0/1 diploid, 0 haploid (chrY, chrM, male chrX outside the PAR), 0/1/1 triploid. Per-region ploidy (PAR, chrX in males, mito) must match the sample karyotype.

. (missing) is NOT reference. ./. = no-call (genotype could not be determined, usually low depth); 0/0 = confidently called homozygous reference. Treating ./. as 0/0 inflates the reference-allele count and biases allele frequencies, missingness, and burden tests. This is load-bearing: never impute ./. as reference. In a gVCF, the ABSENCE of a record also does not mean reference -- see the reference-confidence model below.