View, query, and interpret VCF/BCF variant files with bcftools and cyvcf2. Use when inspecting variants, extracting fields with query format strings, converting VCF/BCF, or correctly reading a field -- QUAL (site) vs GQ (genotype) vs PL/GL likelihoods, AD vs DP and allele balance, GT phasing/ploidy/PS and missing-vs-hom-ref, INFO/FORMAT Number A/R/G semantics, symbolic alleles (<DEL>, <NON_REF>, spanning *) and END, or telling a raw gVCF apart from a filtered callset.
View, query, and interpret VCF/BCF variant files with bcftools and cyvcf2. Use when inspecting variants, extracting fields with query format strings, converting VCF/BCF, or correctly reading a field -- QUAL (site) vs GQ (genotype) vs PL/GL likelihoods, AD vs DP and allele balance, GT phasing/ploidy/PS and missing-vs-hom-ref, INFO/FORMAT Number A/R/G semantics, symbolic alleles (<DEL>, <NON_REF>, spanning *) and END, or telling a raw gVCF apart from a filtered callset.
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
Format
Description
Use Case
VCF
Text format, human-readable
Debugging, small files
VCF.gz
Compressed VCF (bgzip)
Standard distribution
BCF
Binary VCF
Fast processing, large files
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
Column
Description
CHROM
Chromosome
POS
1-based position of the first base in REF (contrast BED's 0-based half-open)
ID
Variant identifier (e.g., rs number) or . if novel
REF
Reference allele (matches the reference exactly)
ALT
Alternate allele(s), comma-separated. * = allele missing due to an overlapping deletion at this site
QUAL
Phred-scaled quality of the ALT assertion, -10*log10 P(no variant); higher = more confident a variant exists (site-level, NOT per-sample)
FILTER
PASS or semicolon-separated filter names. . means filters were not applied
Colon-separated format keys defining per-sample field order
SAMPLE
Colon-separated values matching FORMAT order
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
Field
Level
Scale
Question answered
QUAL (col 6)
Site
Phred: -10*log10 P(no variant)
"Is there ANY variant at this site?"
GQ (FORMAT)
Genotype
Phred, capped at 99
"Is THIS sample's assigned genotype correct?"
PL (FORMAT)
Genotype
Phred, rebased to min=0
Relative likelihood of every possible genotype
GL (FORMAT)
Genotype
log10, <=0, raw
Same info as PL, unscaled (PL = -10*GL, rebased)
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:
Number
One value per
Examples
On multiallelic split
A
ALT allele
AF, AC
take the k-th element for the k-th ALT
R
allele incl. REF
AD
REF value first, then per-ALT (off-by-one vs A)
G
genotype
PL, GL
re-subset via the k*(k+1)/2+j index formula
.
variable/unknown
--
parser CANNOT auto-subset; carried whole onto every record
0
flag (presence only)
--
--
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
Annotation
Meaning
What It Detects
QD
QUAL / allele depth
Low values suggest variant quality not supported by reads
FS
Fisher strand bias (phred-scaled)
Variant reads predominantly on one strand (artifact)
Very negative = alt reads map much worse than ref (suspicious)
ReadPosRankSum
Read position: ref vs alt reads
Very negative = variant only at read ends (misalignment artifact)
Genotype Encoding
Genotype
Meaning
0/0
Homozygous reference (confidently called ref)
0/1
Heterozygous
1/1
Homozygous alternate
1/2
Heterozygous for two different ALT alleles (compound het at multiallelic site)
./.
Missing genotype (no confident call)
0|1
Phased heterozygous (allele before | is on haplotype 1)
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.
Multiallelic Genotypes
At multiallelic sites (e.g. ALT = G,T), allele indices reference the comma-separated ALT list: 0=REF, 1=first ALT, 2=second ALT. 1/2 means one copy of each ALT. Splitting multiallelics into biallelic records with bcftools norm -m- converts 1/2 into two 0/1 records, losing compound-heterozygosity information -- see variant-calling/variant-normalization for caveats.
Symbolic Alleles, END, and Spanning Deletions
Not every ALT spells out a sequence. Symbolic alleles are angle-bracketed placeholders for events whose sequence is not given inline:
ALT
Meaning
<DEL><DUP><INS><INV><CNV>
Structural-variant classes (sequence not spelled out)
<NON_REF>
gVCF: "any allele not yet observed" (reference-confidence model)
<*>
Same role as <NON_REF> in some callers' gVCF/mpileup output
* (bare)
Spanning deletion: allele MISSING because an upstream deletion on ANOTHER line overlaps this position
Two parsing traps:
INFO/END gives the end coordinate of a symbolic/large event. A tool that infers a record's span from len(REF) is WRONG for symbolic alleles -- it must read END. gVCF reference blocks also use END to mark the last position of the band.
The bare * ALT is interpretable only relative to the overlapping deletion on another record; it is not a real alternate allele here. Splitting/subsetting can strand a * from the deletion it references (see variant-calling/variant-normalization).
gVCF and the <NON_REF> Reference-Confidence Model
A gVCF (GATK HaplotypeCaller -ERC GVCF) is fundamentally different from a filtered callset: it emits a record for EVERY position or block, not just variant sites. Non-variant stretches are compressed into END-delimited blocks (bands) grouped by GQ, so a gVCF is not one line per base.
Every record carries a symbolic <NON_REF> ALT with PL/AD computed against "any unseen allele." This lets joint genotyping evaluate a site in THIS sample even when the variant was only discovered in ANOTHER cohort sample -- the <NON_REF> likelihood supplies the evidence.
Its purpose is to distinguish, at every site, confident homozygous reference from no-data/no-call -- solving the missing-vs-reference problem when squaring off a cohort matrix.
A gVCF is NOT ready for analysis; it is an intermediate. It must be joint-genotyped (GenomicsDBImport/CombineGVCFs -> GenotypeGVCFs) to yield a normal VCF. Do NOT filter, annotate, or count variants on a raw gVCF, and never build a multi-sample callset by bcftools merge-ing single-sample project VCFs when gVCF joint-genotyping is available -- merging fabricates hom-ref genotypes. See variant-calling/joint-calling.
bcftools view
Goal: View, subset, and convert VCF/BCF files from the command line.
Approach: Use bcftools view with flags for header control, region selection, sample extraction, and format conversion.
Goal: Convert between VCF, compressed VCF, and BCF, and index for region queries.
Approach: Use bcftools view output flags (-Ov/-Oz/-Ou/-Ob), then bgzip + index.
bcftools view -Ob -o output.bcf input.vcf.gz # VCF -> BCF
bcftools view -Ov -o output.vcf input.bcf # BCF -> VCF
bgzip input.vcf # -> input.vcf.gz (bgzip, NOT gzip)
bcftools index input.vcf.gz # -> .csi index
bcftools index -t input.vcf.gz # -> .tbi (tabix) index
Output Format Options
Flag
Format
-Ov
Uncompressed VCF
-Oz
Compressed VCF (bgzip)
-Ou
Uncompressed BCF (fast piping)
-Ob
Compressed BCF
BCF is the binary encoding of VCF: faster to parse and smaller for large callsets. Region queries (chr1:1-1000) require a bgzipped+indexed VCF or a BCF -- plain .gz (gzip) is not seekable and fails.
cyvcf2 Python Alternative
Goal: Read, query, and write VCF files programmatically in Python.
Approach: Use cyvcf2's VCF reader to iterate variants with attribute access to fields, and Writer to emit filtered output.
"Parse this VCF in Python" -> Open with cyvcf2 and iterate variant records.
Open, Iterate, and Access Fields
from cyvcf2 import VCF
vcf = VCF('input.vcf.gz')
for variant in vcf:
# ALT is a list; QUAL is site-level and may be Noneprint(variant.CHROM, variant.POS, variant.REF, variant.ALT)
print(variant.ID, variant.QUAL, variant.FILTER, variant.var_type)
dp = variant.INFO.get('DP') # INFO field, None if absent
af = variant.INFO.get('AF')
break
vcf.close()
Access Genotypes and Per-Sample Fields
from cyvcf2 import VCF
vcf = VCF('input.vcf.gz')
samples = vcf.samples
for variant in vcf:
# gt_types: 0=HOM_REF, 1=HET, 2=UNKNOWN(missing), 3=HOM_ALT
gts = variant.gt_types
depths = variant.format('DP') # numpy array, one row per sample
gqs = variant.format('GQ') # per-sample genotype quality
ad = variant.format('AD') # per-allele depth, Number=Rprint(dict(zip(samples, gts)))
break
vcf.close()
Note: cyvcf2 codes missing genotypes as gt_types == 2 (UNKNOWN) -- treat that as no-call, never as HOM_REF.
Fetch Region and Read the Header
from cyvcf2 import VCF
vcf = VCF('input.vcf.gz')
print(vcf.samples, vcf.seqnames) # sample names, contig namesfor info in vcf.header_iter():
if info['HeaderType'] == 'INFO':
print(info['ID'], info['Description'])
for variant in vcf('chr1:1000000-2000000'): # requires an indexprint(variant.CHROM, variant.POS)
Write Filtered VCF
from cyvcf2 import VCF, Writer
vcf = VCF('input.vcf.gz')
writer = Writer('output.vcf', vcf) # inherit the input headerfor variant in vcf:
if variant.QUAL isnotNoneand variant.QUAL > 30: # QUAL is site-level
writer.write_record(variant)
writer.close()
vcf.close()
Quick Reference
Task
bcftools
cyvcf2
View VCF
bcftools view file.vcf.gz
VCF('file.vcf.gz')
View header
bcftools view -h file.vcf.gz
vcf.header_iter()
Get region
bcftools view file.vcf.gz chr1:1-1000
vcf('chr1:1-1000')
Query fields
bcftools query -f '%CHROM\t%POS\n'
Loop with properties
Count variants
bcftools view -H file.vcf.gz | wc -l
sum(1 for _ in vcf)
VCF to BCF
bcftools view -Ob -o out.bcf in.vcf.gz
Use Writer
Common Errors
Error
Cause
Solution
no BGZF EOF marker
Not bgzipped (plain gzip)
Recompress with bgzip, not gzip
index required / region query fails
Missing index
Run bcftools index (-t for tabix)
sample not found
Wrong sample name
Check with bcftools query -l
INFO/FORMAT field missing or mistyped
Header out of sync with body
Fix ##INFO/##FORMAT Number/Type; use bcftools +fill-tags
Every hom-alt or missing site vanishes on filter
Treated ././. as failing or as ref
Missing != hom-ref; make missing pass, never impute 0/0
Wrong allele's AF/AD after split
Number=. field not re-subset
Declare the true Number (A/R/G) so bcftools reapportions
Related Skills
variant-calling/variant-calling - Generate VCF from alignments
alignment-files/pileup-generation - Generate pileup for calling
References
Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156-2158. doi:10.1093/bioinformatics/btr330 (VCF format definition)
Danecek P, Bonfield JK, Liddle J, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10(2):giab008. doi:10.1093/gigascience/giab008 (bcftools view/query/norm reference)
