Variant Calling and the VCF Format

Once you have aligned your sequencing reads to a reference genome, the next step is to find where your sample differs from that reference. Those differences are called variants, and finding them is called variant calling. This lesson shows how a variant caller works and how to read the file it produces, the VCF.

From reads to differences

A BAM file is the compressed, binary record of your reads after alignment, storing where each read sits along the reference. A variant caller walks along the genome, stacks up all the reads covering each position into a pileup, and compares the bases they carry against the reference base. Where enough reads consistently disagree with the reference, the caller reports a variant at that position.

Common variant callers

Two callers are used constantly in practice. GATK HaplotypeCaller, from the Broad Institute, reassembles the reads in active regions and is the standard for human germline studies. bcftools, which pairs the mpileup and call commands, is lightweight, fast, and a great place to start.

# mpileup stacks the reads over each position; -f points it at the reference FASTA
# -a FORMAT/DP records each sample's read depth so DP appears in the per-sample column
# call flags the variant sites: -m is the multiallelic caller, -v keeps only variant sites
# -Oz writes a compressed VCF (.vcf.gz); -o names the output file
bcftools mpileup -a FORMAT/DP -f reference.fa sample.bam \
  | bcftools call -mv -Oz -o raw.vcf.gz

# bcftools call leaves FILTER as a dot, so we clean up and label the sites in two steps
# annotate -x FORMAT/PL drops the genotype-likelihood field to keep the sample column simple
# filter fills in FILTER: -e 'QUAL<10' tags low-quality sites q10, and the rest become PASS
bcftools annotate -x FORMAT/PL raw.vcf.gz \
  | bcftools filter -e 'QUAL<10' -s q10 -Oz -o calls.vcf.gz

# Print the #CHROM header line and the variant rows, hiding the ## meta-info lines
bcftools view calls.vcf.gz | grep -v '^##'

#CHROM   POS      ID  REF  ALT   QUAL  FILTER  INFO        FORMAT  sample1
20       14370    .   G    A     29    PASS    DP=14       GT:DP   0/1:14
20       17330    .   T    A     3     q10     DP=11       GT:DP   0/1:11
20       1110696  .   A    G     67    PASS    DP=10       GT:DP   1/1:10
20       1230237  .   CTT  C     50    PASS    INDEL;DP=9  GT:DP   1/1:9
20       1234567  .   G    GACT  42    PASS    INDEL;DP=8  GT:DP   0/1:8

Reading the columns

Every variant row starts with eight fixed columns. CHROM and POS give the chromosome and its 1-based position; ID is a known identifier such as an rs number, or a dot when none exists; REF is the reference base or bases and ALT is the alternate found in the sample; QUAL is a Phred-scaled confidence score where higher means more certain; FILTER reads PASS, or the name of a check the site failed, or a dot before any filtering has run; and INFO packs site-level annotations, shown here trimmed to DP (total read depth) even though a real file carries many more. A ninth column, FORMAT, then defines the per-sample fields that follow it, so GT:DP means every sample column lists a genotype followed by that sample's read depth.

When REF and ALT are both a single base, like G to A, the variant is a SNP, a single-nucleotide change. When they differ in length, like CTT to C (a deletion of TT) or G to GACT (an insertion of ACT), the variant is an indel. The genotype GT reports the two inherited copies, called alleles, as numbers where 0 is the reference and 1 is the first ALT, so 0/1 is heterozygous (one of each), 1/1 is homozygous alternate (both copies changed), and 0/0 would match the reference on both copies.