How to Analyze a VCF File From Your Own Genome

By the end of this guide you will have a bgzip-compressed, tabix-indexed, normalized, annotated VCF for a single human genome, a set of reproducible filter expressions you trust, and a few flat TSV extracts small enough to open in a spreadsheet or load into pandas. VCF stands for Variant Call Format, the standard text format for recording the positions where a sequenced genome differs from the reference.

To follow along you need a Unix shell with roughly 100 GB of free disk, bcftools 1.19 or newer, and htslib for bgzip and tabix. You also need the reference FASTA your file was called against, either GRCh38 or hg38 with the same contig naming as your VCF. For annotation you need either Ensembl VEP with a local cache or SnpEff. Everything below runs on a laptop with 16 GB of RAM except the VEP cache download, which wants about 30 GB of disk. A single-sample whole-genome VCF is roughly 1 to 2 GB uncompressed and contains on the order of 4 to 5 million variant records, so plan on working with it through indexed queries rather than by scrolling.

1. Confirm which kind of .vcf file you have

Before anything else, make sure the file in front of you is the kind you think it is. The extension .vcf is overloaded, and this is the source of most confusion. A vCard file also uses .vcf, and it is what you get when you export contacts from a phone or mail client. It begins with the line BEGIN:VCARD, and it is what Windows tries to import into your address book when you double-click it. A Variant Call Format file, by contrast, begins with ##fileformat=VCFv4.2 or something similar. One command settles the question:

head -c 200 myfile.vcf
# or, if it is compressed:
zcat myfile.vcf.gz | head -c 200

If you see BEGIN:VCARD, you have a contacts file and nothing in this guide applies to it. You can open it in a text editor to view its contents without importing it, or convert it to CSV with a vCard parser. If you see ##fileformat=VCF, you have genomic variant calls and should keep reading.

One expectation to set early: a genomic VCF will not open usefully in Excel by double-clicking. It is too large for the 1,048,576-row limit, the header lines confuse the delimiter sniffer, and the INFO column is a semicolon-packed key-value blob that no spreadsheet will parse. Step 7 covers the right way to get a subset into a spreadsheet.

2. Compress, index, and validate before you touch the data

The first real step is to put the file into the form every downstream tool expects, then check that its contents are plausible. Tools worth using expect a block-gzipped file with a companion index, which allows random access by genomic coordinate instead of a linear scan through the whole file. Standard gzip will not work here, because bgzip writes a series of independently compressed blocks whose offsets the index records.

# if you received a plain .vcf
bgzip -@ 4 myfile.vcf                  # -> myfile.vcf.gz
tabix -p vcf myfile.vcf.gz             # -> myfile.vcf.gz.tbi

# if you received a .gz that is not bgzipped, recompress
gunzip -c myfile.vcf.gz | bgzip -@ 4 > myfile.bgz.vcf.gz

For genomes with contigs longer than 512 Mb, or if you simply prefer the newer index format, use bcftools index --csi instead. With the index in place, run a structural sanity check and collect the basic counts you will refer back to all day:

bcftools index -n myfile.vcf.gz        # number of records
bcftools stats myfile.vcf.gz > stats.txt
grep -E "^SN|^TSTV" stats.txt

Three numbers in stats.txt tell you whether the file is plausible. A single human genome should show something near 4 to 5 million total records, roughly 4 to 4.5 million SNVs, and a transition/transversion ratio around 2.0 to 2.1 genome-wide. A Ts/Tv much below 1.8 usually means the call set is contaminated with false positives arising from low coverage or an aggressive calling threshold. An exome or targeted panel will look different by design, showing tens of thousands of variants instead, with a higher Ts/Tv because coding sequence is enriched for transitions.

3. Read the header, because it tells you what you are allowed to conclude

The header is both the file’s provenance record and its data dictionary, and reading it in full once will save you a great deal of guesswork later.

bcftools view -h myfile.vcf.gz | less

There are four things to look for. First, ##reference= and the ##contig lines: if the contigs are named chr1 and your annotation resources use 1, every join will silently return nothing. Second, the ##source= or ##GATKCommandLine lines, which name the caller (GATK HaplotypeCaller, DeepVariant, DRAGEN, bcftools call) and often its exact parameters. Third, the ##FORMAT definitions, which tell you which per-genotype fields exist. You want at minimum GT (genotype), DP (depth), GQ (genotype quality), and AD (allelic depths). Fourth, check whether the file is a gVCF, indicated by <NON_REF> alleles and END= records covering reference blocks. A gVCF is not a variant list and must be genotyped first with gatk GenotypeGVCFs before filtering makes sense.

The identity of the caller matters more than people expect. Large consortia had to define a functional equivalence standard precisely because nominally equivalent pipelines produced call sets that disagreed enough to break joint analysis. The agreed specification pins down aligner, caller version, and post-processing steps rather than leaving them to the analyst¹. If you plan to compare two VCFs of the same person produced by different pipelines, expect real discordance concentrated in indels and repetitive regions. Use a comparison tool that normalizes representation rather than matching on coordinates. hap.py against a Genome in a Bottle truth set is the standard for benchmarking. Newer tools built for this problem align the variant representations before declaring a difference, which avoids the large class of false discordances that come from the same indel written two ways².

4. Normalize: split multiallelics and left-align indels

Normalization puts every variant into a single canonical representation, and skipping it is the most common cause of quietly wrong results. The same deletion can be written at different positions with different padding bases, and multiallelic records (an ALT of A,AT, for instance) break naive filtering because the INFO and FORMAT arrays are stored per allele. Normalization gives you one record per alternate allele, each left-aligned and trimmed against the reference.

bcftools norm -m -any -f GRCh38.fa -c w -Oz -o norm.vcf.gz myfile.vcf.gz
bcftools index norm.vcf.gz

Each flag earns its place. -m -any splits multiallelic sites into biallelic records, -f supplies the reference so indels can be left-aligned and trimmed, and -c w warns when the REF allele in your file disagrees with the reference FASTA instead of silently continuing. Watch the summary line that bcftools prints: a handful of REF mismatches usually indicates a contig naming or build mismatch, while thousands means you have the wrong reference entirely. If you need to collapse duplicate records afterward, add a second pass with bcftools norm -d exact.

5. Filter on quality, and know what each threshold throws away

Filtering happens at two levels, and it helps to keep them distinct. The FILTER column carries the caller’s own verdict, and for a well-calibrated pipeline you should trust it as the first gate. Beyond that, the per-genotype fields let you remove calls that are technically PASS but supported by too little evidence.

bcftools view -f PASS -Oz -o pass.vcf.gz norm.vcf.gz

# a conservative single-sample genotype filter
bcftools view -i 'FMT/DP>=10 & FMT/GQ>=20 & (FMT/AD[0:1])/(FMT/DP)>=0.25' \
  -Oz -o pass.q.vcf.gz pass.vcf.gz
bcftools index pass.q.vcf.gz

The allele-balance term is the most useful of the three. A heterozygous site with 30 reads should show roughly 15 supporting the alternate allele, so a PASS het call with an alternate fraction of 0.08 is far more likely to be a mapping artifact or a somatic event in blood than a germline heterozygote. For a 30x whole genome, meaning an average of about 30 reads covering each position, DP>=10 removes something like a few percent of sites, concentrated in high-GC promoters, segmental duplications, and the usual hard regions.

Be clear about what this filter does not do. It says nothing about positions that were never callable in the first place. If a region had zero coverage, there is no record to filter, and absence of a variant call is not evidence of a reference genotype. Compute your callable footprint separately with mosdepth on the BAM or CRAM and keep it alongside the VCF.

Automated exome and genome pipelines have applied essentially this shape of filter for a decade. They chain alignment, calling, and a quality gate before any interpretation step³. The specific thresholds remain a judgment call about your tolerance for false positives versus false negatives, and they should be stated explicitly whenever you report a result.

6. Annotate, because coordinates alone mean nothing

A filtered VCF is a list of positions. To reason about any of them you need gene and transcript consequence, population allele frequency, and clinical assertions. Ensembl VEP with a local cache is a good default. It runs offline, is versioned, and gives you per-transcript consequences rather than a single collapsed guess.

vep --offline --cache --dir_cache ~/.vep --assembly GRCh38 \
    --fasta GRCh38.fa --input_file pass.q.vcf.gz --format vcf \
    --vcf --compress_output bgzip --output_file annot.vcf.gz \
    --everything --pick_allele_gene --fork 4 \
    --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax

--everything turns on SIFT, PolyPhen, canonical and biotype flags, and existing variant IDs. --pick_allele_gene reduces the output to one consequence per allele per gene, which keeps the file manageable at the cost of hiding alternate transcripts you may care about later. Add --custom blocks for ClinVar and any other coordinate-keyed resource. If VEP’s cache is too heavy for your setup, bcftools csq with an Ensembl GFF3 gives fast, haplotype-aware consequence calls with far less setup. The annotation set is thinner, however.

Population frequency is the highest-yield annotation of the lot. Most of your 4 to 5 million variants are common, and filtering to a gnomAD allele frequency below 0.1 percent typically leaves a few thousand rare variants genome-wide and a few hundred rare coding ones. That reduction, combined with consequence filtering, is what makes a genome interpretable at all⁴.

A practical note on scope: whole-genome data contains far more non-coding variation than exome data, but the additional clinically interpretable yield over a well-analyzed exome is modest, and reanalysis of existing exome data with updated annotation databases recovers much of the difference⁵. Treat the non-coding fraction as a research resource rather than a source of answers. Any variant you are tempted to treat as medically meaningful needs confirmation in a clinical laboratory and interpretation by a qualified clinician or genetic counselor. Research-grade VCFs are not diagnostic, and ClinVar assertions vary widely in evidence quality.

7. Extract a flat table you can open in a spreadsheet

With annotation in place, the spreadsheet question finally has a good answer. Do not convert the whole VCF. Filter first, then project the columns you want into a TSV with bcftools query, which writes exactly the fields you name and nothing else.

bcftools +split-vep annot.vcf.gz -c Consequence,SYMBOL,gnomADg_AF:Float -Ou \
 | bcftools view -i 'gnomADg_AF<0.001 || gnomADg_AF="."' -Ou \
 | bcftools view -i 'Consequence~"missense_variant" || Consequence~"stop_gained" || Consequence~"frameshift_variant"' -Ou \
 | bcftools query -H -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%FILTER\t%INFO/SYMBOL\t%INFO/Consequence\t%INFO/gnomADg_AF\t[%GT\t%DP\t%GQ]\n' \
 > rare_coding.tsv

The -H flag writes a header row with column numbers, which you should strip or rename before import. Open the result with Excel’s Data > From Text/CSV dialog and set the delimiter to Tab and every column type to Text, because otherwise the importer will turn gene symbols like MARCH1 and SEPT9 into dates and coordinates into floats. In pandas, pd.read_csv("rare_coding.tsv", sep="\t", dtype=str) avoids the same trap. For a whole genome this pipeline typically yields a few hundred to a couple of thousand rows, which is a spreadsheet-sized problem.

If you would rather point and click through filters than write expressions, tools exist for exactly this. VCF-Miner provides a graphical interface for building filters over INFO and annotation fields on multi-sample VCFs⁶, and VCF-Explorer is designed to filter and summarize whole-genome-scale VCFs without loading them into a database first⁷. We still prefer the command line, because a filter expression saved in a shell script is reproducible and a sequence of mouse clicks is not.

8. Move to a columnar format once queries get repetitive

Text VCF is a fine interchange format and a poor analysis format. Once you are running many queries, or joining across samples and time points, it pays to convert to something columnar. For R users, SeqArray stores WGS calls in a GDS container with per-field compression and gives order-of-magnitude speedups on genotype-wide scans compared with parsing text, along with a data model that handles millions of variants without loading everything into memory⁸. For Python and PLINK-based workflows, one command gets you there:

plink2 --vcf pass.q.vcf.gz --make-pgen --out mysample

The PGEN format is compact and is what you want for any allele-frequency or polygenic score arithmetic. Whichever route you take, keep the annotated VCF as the source of truth and treat the columnar copies as derived artifacts you can regenerate at will.

9. Run identity and integrity checks before you believe anything

Two failure modes cost more time than any others: analyzing the wrong person’s file, and analyzing a file whose sample metadata is wrong. Both are cheap to rule out, so do it before you draw conclusions.

bcftools query -l pass.q.vcf.gz                    # sample names in the file
bcftools guess-ploidy -g GRCh38 pass.q.vcf.gz      # inferred sex from X/Y
bcftools stats -s - pass.q.vcf.gz | grep "^PSC"    # per-sample het/hom counts

Inferred sex that contradicts the expected value points to a sample swap or a contaminated library, and it invalidates everything downstream. Heterozygosity far above the expected genome-wide rate, which is roughly 0.001 per base for an outbred individual, also suggests contamination. If you have two VCFs that should belong to the same person, check concordance with bcftools gtcheck, which reports discordance rates between sample pairs. Expect greater than 99 percent concordance on common SNVs for a true match, even across different pipelines.

Common problems

The failures below account for most of the time lost in this workflow, and each has a short diagnosis.

tabix fails with “the index file is older than the data file” or “not a BGZF file.” The file was compressed with gzip rather than bgzip. Recompress as shown in step 2 and rebuild the index.

Annotation returns nothing for every variant. This is almost always a contig naming mismatch (1 versus chr1) or a build mismatch (GRCh37 coordinates run against a GRCh38 cache). Rename contigs with bcftools annotate --rename-chrs and a two-column mapping file, and confirm the build by spot-checking a known variant’s coordinates.

Filter expressions silently match zero records. Check whether the field lives in INFO or FORMAT and prefix it accordingly (INFO/DP versus FMT/DP), and remember that after splitting multiallelics some array-valued INFO fields keep their original arity unless the caller wrote Number=A.

Excel mangles the data. Gene symbols become dates, long identifiers become scientific notation, and rows beyond 1,048,576 vanish without warning. Import as Text through the Data tab, and never edit a genomic table in place in a spreadsheet you intend to reuse.

Two VCFs of your own genome disagree on thousands of indels. This is expected. Most of the difference is representation rather than biology, and a normalizing comparison will confirm as much². Different pipelines also genuinely differ, which is why harmonized specifications exist¹.

The file has no GQ or AD fields. Some pipelines strip per-genotype annotations to save space, which removes your ability to filter on evidence. Ask for the original gVCF or the CRAM, from which both can be regenerated.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Allison A. Regier, Yossi Farjoun, David E. Larson, et al. Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects. Nature Communications, 2018. https://doi.org/10.1038/s41467-018-06159-4 ↩ ↩²

2. James M. Holt, Christopher T. Saunders, Egor Dolzhenko, et al. Aardvark: sifting through differences in a mound of variants. Genome Biology, 2026. https://doi.org/10.1186/s13059-026-04165-0 ↩ ↩²

3. Yunfei Guo, Xiaolei Ding, Yufeng Shen, et al. SeqMule: automated pipeline for analysis of human exome/genome sequencing data. Scientific Reports, 2015. https://doi.org/10.1038/srep14283 ↩

4. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

5. Ahmed Alfares, Taghrid Aloraini, Lamia Al subaie, et al. Whole-genome sequencing offers additional but limited clinical utility compared with reanalysis of whole-exome sequencing. Genetics in Medicine, 2018. https://doi.org/10.1038/gim.2018.41 ↩

6. Steven N. Hart, Patrick Duffy, Daniel J. Quest, et al. VCF-Miner: GUI-based application for mining variants and annotations stored in VCF files. Briefings in Bioinformatics, 2015. https://doi.org/10.1093/bib/bbv051 ↩

7. Mete Akgün, Hüseyin Demirci. VCF-Explorer: filtering and analysing whole genome VCF files. Bioinformatics, 2017. https://doi.org/10.1093/bioinformatics/btx422 ↩

8. Xiuwen Zheng, Stephanie M Gogarten, Michael Lawrence, et al. SeqArray—a storage-efficient high-performance data format for WGS variant calls. Bioinformatics, 2017. https://doi.org/10.1093/bioinformatics/btx145 ↩