How to Count Variants in a VCF File

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-10

A VCF file does not have one number in it. It has at least three, and they can differ by a factor of five or more on the same file. By the end of this walkthrough you will be able to produce all of them correctly: the number of records (rows), the number of alleles (the split ALT values), and the number of variant calls carried by each sample. You will also be able to restrict any of those counts to PASS records, to a chromosome, to a BED interval, or to the overlap between two VCF files.

You need a VCF file, bcftools and tabix from the htslib project, and a shell with awk, sort, and uniq. A reference FASTA is useful if you plan to validate REF alleles, and GATK4 or vcftools are optional cross-checks. Nothing here requires a cluster. A 100,000-record VCF counts in under a second, and a 5-million-record VCF counts in a few seconds because the record count is stored in the index.

One warning before the commands. A variant count is a property of a file and a pipeline, not a property of a genome. Two teams can sequence the same sample and report 4.1 million and 4.9 million variants, both correctly, because their callers, filters, and reference builds differ. Report the count together with the exact command, the bcftools version, and the filters applied. That habit prevents most of the arguments that variant counts cause.

What "Count the Variants" Actually Means

A VCF (Variant Call Format) file is a plain-text format defined by the GA4GH VCF specification, currently at versions 4.2 and 4.3. The structure is strict and worth knowing before you count anything:

• Lines beginning with ## are metadata header lines. They declare the reference contigs, the FILTER names, the INFO and FORMAT fields, and the file format version.
• One line beginning with a single # is the column header. The first eight columns are fixed: CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO. Everything after INFO is a FORMAT column followed by one column per sample.
• Every remaining line is a data record. One record is one genomic position with one or more alternate alleles.

That last point is the source of the confusion. The ALT column can hold a comma-separated list. A record with ALT=G,T describes two alternate alleles at one position. So:

• The number of records (also called sites) is the number of data lines.
• The number of alleles is the sum of the ALT entries across all records. It is always greater than or equal to the record count.
• The number of genotypes or sample-level calls is the number of sample columns that carry a non-reference genotype at a record. In a 1,000-sample cohort this number is enormous and almost never what a headline figure means.

There is a fourth number that trips people up: gVCF reference blocks. A gVCF (genomic VCF) stores confidence information for positions where the sample matches the reference, usually as records with ALT=<NON_REF>. Those records are not variants. If you count every line in a gVCF you will report millions of "variants" that are actually non-variant blocks.

A record can also carry symbolic alleles such as <DEL>, <INS>, <DUP>, and <NON_REF>, or breakend notation. These are valid VCF but they are not SNPs and not simple indels, and most counting conventions put them in an "other" bin.

The table below summarizes which command answers which question. Each command is explained in detail in the walkthrough.

Use this table as a decision aid, then read the step that applies to you.

What You Need Before You Start

Install the tools. On a Debian or Ubuntu system the packages are in the standard repositories, and on macOS or inside a conda environment the bioconda channel is the usual route:

# Debian or Ubuntu
sudo apt-get update
sudo apt-get install -y bcftools tabix

# conda, mamba, or pixi
mamba install -c bioconda -c conda-forge bcftools htslib tabix

# confirm
bcftools --version
tabix --version

The bcftools, version output prints two lines: the bcftools version and the htslib version it was built against. Write both down. Classification of variant types and the exact behavior of filters have changed across releases, so a count without a version is not reproducible.

Next, decide which file you are counting. Counting works on three input flavors:

1. Uncompressed VCF (.vcf). Readable with grep and awk, but slow for large files and not indexable.
2. bgzipped VCF (.vcf.gz). Compressed with bgzip, which is gzip-compatible. Random access requires a companion index (.tbi from tabix or .csi from bcftools).
3. BCF (.bcf). The binary encoding of the same records, defined alongside the VCF specification in htslib. Smaller and faster to parse.

To convert an uncompressed VCF into the canonical indexed form:

bcftools view -Oz -o example.vcf.gz example.vcf
bcftools index -t example.vcf.gz
ls -lh example.vcf.gz example.vcf.gz.tbi

-Oz means output compressed VCF. bcftools index -t builds a tabix-style .tbi index. Add -f to overwrite an existing index. For contigs longer than 512 Mb, or for BCF files, use bcftools index -c to build a .csi index instead, because tabix indices use 16-bit binning that cannot address very long chromosomes.

If you only ever stream through a file once, you can count without an index. If you want region queries or the instant count from bcftools index -n, you need the index.

Sanity-check the header before counting:

bcftools view -h example.vcf.gz | head -20
bcftools view -h example.vcf.gz | grep -c '^##contig'

The first line should be ##fileformat=VCFv4.2 or VCFv4.3. The contig count tells you whether the header declares chromosomes. If the contig count is zero you will still get a count, but region queries and downstream tools may misbehave, and htslib will emit a warning during parsing.

A Small Worked Example

All the commands below run against this file. Save it as example.vcf, then compress and index it with the two commands above.

##fileformat=VCFv4.2
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	S1	S2	S3
chr1	1000	.	A	G	50	PASS	DP=30	GT	0/1	0/0	1/1
chr1	2000	.	C	T	40	PASS	DP=28	GT	0/1	0/1	0/0
chr1	3000	.	G	A,T	60	PASS	DP=35	GT	1/2	0/1	0/0
chr1	4000	.	T	TA	12	LowQual	DP=9	GT	0/1	0/0	./.
chr2	500	.	A	<DEL>	.	PASS	END=800	GT	0/1	0/0	0/0
chr2	900	.	C	G	99	PASS	DP=41	GT	1/1	0/1	0/0
chr2	1500	.	N	<NON_REF>	.	.	END=1600	GT	0/0	0/0	0/0

This file has 7 data records, 8 ALT alleles (the G>A,T record contributes two), 4 SNP records, 1 indel record, 2 symbolic records, 5 PASS records, and 10 non-reference genotypes. Every command in the next eight steps reproduces one of those numbers. If a command gives you a different answer, the difference itself is diagnostic.

Count Variants in a VCF File Step by Step

Step 1: Count the data records

The portable one-liner is:

bcftools view -H -Ov example.vcf.gz | wc -l

7

-H suppresses the header lines. -Ov forces uncompressed VCF on stdout so that wc -l counts text lines and not binary bytes. Without -Ov, bcftools picks the output type from the input type, and a BCF input would send binary to wc -l.

Once the file is indexed, the index already stores the record count:

bcftools index -n example.vcf.gz

7

This reads the index and returns immediately, even on multi-gigabyte files. It is the single fastest way to count variants in a VCF file that has an index, and it is the number to quote in a methods section when you have not modified the records.

For an uncompressed VCF with no bcftools available, grep works:

grep -vc '^#' example.vcf

7

-v inverts the match and -c counts. On a bgzipped file, zcat example.vcf.gz | grep -vc '^#' also works because bgzip output is valid gzip, but it decompresses the whole file.

Step 2: Count declared alleles instead of records

A multiallelic record is one line and two or more variants. If your question is "how many alternate alleles are described here", count the ALT entries:

bcftools view -H -Ov example.vcf.gz | awk '
  { n = split($5, a, ",")
    s += n }
  END { print s }
'

8

$5 is the ALT column. split returns the number of comma-separated pieces, and the loop accumulates them. This ignores the ALT value ".", which denotes a non-variant record and returns 0 from split. That is the behavior you want for a gVCF block such as chr2 1500.

The bcftools-native equivalent splits the records and then counts:

bcftools norm -m-any example.vcf.gz | bcftools view -H | wc -l

8

bcftools norm -m-any splits records with multiple ALT alleles into one biallelic record per allele. The output is 8 records for the same reason the awk script returns 8.

Report allele counts and record counts as different numbers. A VCF with many multiallelic sites, which is common in multi-sample cohorts, can have 20 percent more alleles than records. Mixing the two inflates or deflates a headline figure depending on which command produced it.

Step 3: Count by variant class

bcftools classifies each record from its REF and ALT alleles into SNPs, indels, MNPs, and others:

for t in snps indels mnps other
do
  n=$(bcftools view -H -v "$t" example.vcf.gz | wc -l)
  printf '%-8s %s\n' "$t" "$n"
done

snps     4
indels   1
mnps     0
other    2

-v means include these types. -V means exclude them, so bcftools view -V snps returns the non-SNP records. The four counts sum to the record total, which is a useful check: 4 + 1 + 0 + 2 = 7.

The two "other" records in the example are the symbolic deletion <DEL> and the gVCF block <NON_REF>. Neither is a SNP or a simple indel. If you want the real variant count from a gVCF, remove the blocks first:

bcftools view -e 'ALT="<NON_REF>"' -H -Ov example.vcf.gz | wc -l

6

The expression engine accepts -e for exclude and -i for include. Filtering on ALT cuts every non-variant block, which is what turns a gVCF line count into a variant count.

An alternative classification route uses the %TYPE format field:

bcftools query -f '%TYPE\n' example.vcf.gz | sort | uniq -c

      2 OTHER
      1 INDEL
      4 SNP

Note that this tallies records, so the multiallelic G>A,T record appears once as SNP.

Step 4: Count only FILTER=PASS records

Raw counts include calls that failed quality control. For a headline figure, count PASS records only:

bcftools view -f PASS -H -Ov example.vcf.gz | wc -l

5

-f selects records whose FILTER column matches one of the listed values. The distinction between PASS and . matters. PASS means the record was evaluated and passed. . means no filter was applied, often because the caller did not evaluate it. They are not interchangeable. To include both:

bcftools view -f .,PASS -H -Ov example.vcf.gz | wc -l

6

List the filter values present so you know what you are excluding:

bcftools query -f '%FILTER\n' example.vcf.gz | sort | uniq -c | sort -rn

      5 PASS
      1 .
      1 LowQual

If your VCF has no FILTER beyond PASS and ., the filter step changes little. If it has dozens of soft filters carried through from multiple callers, it can change the count by a large margin, and you should state which filters you kept.

Step 5: Count per chromosome or per region

Because coordinate-sorted VCFs group records by chromosome, uniq -c on the CHROM column gives a per-chromosome table without any sorting:

bcftools query -f '%CHROM\n' example.vcf.gz | uniq -c

      4 chr1
      3 chr2

If the file is not sorted, or you are unsure, use the sorted variant:

bcftools query -f '%CHROM\n' example.vcf.gz | sort | uniq -c

For a region, use the index:

bcftools view -H -r chr1:1-2500 example.vcf.gz | wc -l

2

-r takes one or more regions in chrom:start-end form or a tab-delimited file of regions, and it requires an index. -R takes a file in VCF-style coordinates. A BED file needs care because BED is 0-based and half-open while VCF is 1-based and closed. A BED interval chr1 1000 2000 covers 1-based positions 1001 through 2000, so a variant at POS 1000 is outside it, not inside. If you feed BED coordinates to -r directly you shift the window by one base.

-R on a BED file includes every record that overlaps an interval, including records whose POS is before the interval start when END extends into it. -T (targets) instead streams the file and tests each position, which is faster on unsorted or very large inputs but does not include overlapping records by default.

A per-interval loop with correct BED conversion:

while IFS=$'\t' read -r chrom start end name
do
  from=$((start + 1))
  n=$(bcftools view -H -r "${chrom}:${from}-${end}" example.vcf.gz | wc -l)
  printf '%s\t%s\n' "$name" "$n"
done < regions.bed

The IFS setting keeps tab-delimited fields intact, and start + 1 converts the 0-based BED start to a 1-based VCF position.

Step 6: Count variants carried by each sample

This is the step where "count variants" and "count carriers" diverge. A record exists once in the file but may be heterozygous in one sample, homozygous in another, and reference in a third. Two correct approaches exist.

The first uses genotype strings directly:

bcftools query -f '[%SAMPLE\t%GT\n]' example.vcf.gz | awk '
  $2 !~ /^(0\/0|0\|0|0|\.\/\.|\.\|\.|\.)$/ { print $1 }
' | sort | uniq -c

      6 S1
      3 S2
      1 S3

The bracketed format [ ... ] repeats the enclosed string once per sample, so this expands every record into one line per sample. The awk filter drops homozygous reference, haploid reference, and fully missing genotypes. Anything else counts as a non-reference call.

The second approach subsets and filters with bcftools:

for s in S1 S2 S3
do
  n=$(bcftools view -H -s "$s" -c1 example.vcf.gz | wc -l)
  printf '%s\t%s\n' "$s" "$n"
done

S1	6
S2	3
S3	1

-s keeps only the named sample, and -c1 keeps records with at least one non-reference allele among the samples still in the file. After -s S1, the allele count is computed over S1 alone, so the result is a per-sample carrier count.

Two traps live in this step. The first is partially missing genotypes such as 0/. or ./1. The regex above counts 0/. as non-reference because it is not an exact match for a reference genotype, which inflates counts on low-depth data. The second is the * spanning-deletion allele, which can appear in the ALT column and produces genotypes that are not real alternate calls. For publication figures, prefer bcftools view -c1 per sample, or run bcftools stats -s - and take the per-sample block, which reports counts consistently for all samples in one pass.

bcftools stats -s - prints a PSC block with one line per sample covering records, SNPs, indels, and genotype class counts. It is the fastest way to get a per-sample table without writing a loop.

Step 7: Count shared and unique variants between two VCFs

Comparing two VCFs needs indexes on both and a defined rule for matching. bcftools uses chromosome, position, and REF/ALT by default:

bcftools isec -n=2 -w1 sampleA.vcf.gz sampleB.vcf.gz | bcftools view -H | wc -l

-n=2 keeps sites present in both files, and -w1 writes the matching sites in the representation used by the first file. The pipeline counts the shared records.

For all four categories at once, write to a directory:

bcftools isec -p isec_out -n=2 sampleA.vcf.gz sampleB.vcf.gz
ls isec_out

The directory receives numbered VCF files plus a README.txt that states which set each number holds. Count each file with bcftools view -H <file> | wc -l and you have shared sites, sites unique to A, and sites unique to B. Matching is positional, so a different reference build, different chromosome naming, or a different representation of the same indel will all reduce the shared count. Normalize both files with the same reference before comparing them:

bcftools norm -f reference.fa -m-any -Oz -o sampleA.norm.vcf.gz sampleA.vcf.gz
bcftools norm -f reference.fa -m-any -Oz -o sampleB.norm.vcf.gz sampleB.vcf.gz
bcftools index -t sampleA.norm.vcf.gz
bcftools index -t sampleB.norm.vcf.gz

-f supplies the reference FASTA, and -m-any splits multiallelics into a consistent biallelic representation.

Step 8: Check the number against a second method

Never report a single number produced by a single command. Cross-check it:

# 1. index count versus streamed count
bcftools index -n example.vcf.gz
bcftools view -H -Ov example.vcf.gz | wc -l

# 2. bcftools stats record line
bcftools stats example.vcf.gz | grep 'number of records'

# 3. duplicate positions
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' example.vcf.gz \
  | sort | uniq -d | wc -l

# 4. header integrity
bcftools view -h example.vcf.gz | grep -c '^##contig'

# 5. exact file identity
md5sum example.vcf.gz

If the index count and the streamed count disagree, the index is stale. Rebuild it with bcftools index -f. If duplicate positions appear, you are probably counting a file produced by merging several callers without deduplication. Remove exact duplicates with bcftools norm -d exact and recount. If bcftools stats reports a record count that differs from your streamed count, check whether you ran stats on the same file and whether an expression filter was applied.

Record the md5 checksum of the exact file you counted. Variant counts get challenged, and the checksum plus the command plus the bcftools version settles the question without rerunning anything.

Command Reference Table

Common Mistakes and Limitations

Counting lines in a compressed file with wc -l. Running wc -l variants.vcf.gz counts newline bytes in the compressed stream. The result is a small, arbitrary number with no relationship to the variant count. Use bcftools index -n, bcftools view -H | wc -l, or zcat file.vcf.gz | grep -vc '^#'.

Forgetting that bcftools view can emit BCF. If the input is BCF and you do not pass -Ov, the piped output is binary and wc -l counts bytes. Always pass -Ov when a text line count follows.

Counting header lines. grep -c '' file.vcf counts metadata lines as variants. On a VCF with 200 ## lines, a small file inflates noticeably. Use grep -vc '^#'.

Indexing an uncompressed VCF. The command fails with an error in the family of [E::bcf_index_build] the file is not BGZF compressed, cannot index. bgzip is required. Convert first with bcftools view -Oz -o out.vcf.gz in.vcf.

Region queries without an index. bcftools and htslib report that they cannot retrieve an index, typically with a line beginning [E::idx_find_and_load] Could not retrieve index file. Build the index or drop the region filter.

Chromosome naming mismatches. -r chr1:1-1000 returns zero records with no error when the file uses 1 instead of chr1. Check the header contigs before concluding a region is empty. Rename consistently with bcftools annotate, rename-chrs mapping.txt rather than editing files by hand.

Confusing sites with alleles. A pipeline that splits multiallelics will report more variants than one that does not, on identical data. State which convention you used.

Counting gVCF blocks as variants. A gVCF can easily contain 10 to 100 times more records than variants. Filter on ALT="<NON_REF>" before counting.

Counting before filtering. Raw record counts include low-quality calls. Report both the pre-filter and post-filter numbers, or state the filter explicitly.

Counting a joint-called VCF and calling it per-sample. A cohort VCF with 500 samples and 3 million records contains up to 1.5 billion sample-level genotypes. Those are different units and must be labeled differently.

Duplicate records after merging. Merging outputs from several callers without deduplication inflates counts. Detect duplicates with uniq -d on the CHROM/POS/REF/ALT key and remove them with bcftools norm -d exact.

Off-by-one BED handling. BED is 0-based half-open, VCF is 1-based. Copying BED coordinates into -r silently shifts every window by one base, which matters most for small target intervals.

Contig and INFO headers that are missing. htslib warns with lines like [W::vcf_parse] Contig 'chr1' is not defined in the header. (Quick workaround: index the file with tabix.) and [W::vcf_parse_info] INFO 'DP' is not defined in the header, assuming Type=String. The count still runs. Downstream tools may refuse the file, so treat these warnings as a data quality signal.

Trusting a single command. Two commands that disagree indicate a stale index, a different file, or a filter you forgot. Resolve the disagreement before reporting either number.

The main limitation is conceptual. A count describes a file, not a genome. It is inflated by sequencing artifacts and false positives, deflated by stringent filters, and completely dependent on the caller, the reference build, and the alignment quality. Counts from different pipelines are not directly comparable, and a difference between two counts is not evidence of a biological difference. A variant count also says nothing about pathogenicity or clinical meaning. Whether a specific variant matters for a patient is a question for a clinical geneticist or the treating clinician, not something a line count can answer.

Where the VCF Sits Among Bioinformatics File Formats

Variant counting makes more sense if you know what sits upstream and downstream of the VCF. The common bioinformatics file formats form a chain, and each one answers a different question.

FASTA stores nucleotide or protein sequence. It uses a > header line followed by sequence lines wrapped at a fixed width, and it carries no quality information. In a variant pipeline, FASTA is the reference genome. The REF column of every VCF record is supposed to match the base at that position in the FASTA. You can verify that with bcftools norm -f reference.fa -c w file.vcf.gz, which warns when a REF allele disagrees with the reference. A wrong or mismatched reference build is one of the quiet causes of a count that looks plausible but is wrong.

FASTQ stores raw sequencing reads. Each read occupies four lines: an @ header, the sequence, a + separator, and one ASCII character per base encoding the Phred quality score. FASTQ is the input to alignment. FASTA and FASTQ are often confused because both start with a header character and both hold sequence, but the split is simple. FASTA is one sequence per entry with no qualities, used as a reference or a product. FASTQ is millions of short reads with per-base qualities, used as raw data. If the file has four lines per record and a quality string, it is FASTQ.

BAM and CRAM store reads aligned to the reference. They are the input to variant calling.

VCF and BCF store the calls. BCF is the binary encoding of VCF with the same logical content, produced with bcftools view -Ob and read transparently by bcftools.

gVCF is a VCF with reference blocks, used for joint genotyping across cohorts. It is a VCF with extra records, not a separate format.

Tabix index (.tbi) or CSI index (.csi) sits beside a bgzipped VCF and enables random access. Without it, -r and -R queries cannot run and bcftools index -n cannot return a count.

The practical rule: count records from a VCF or BCF, never from a FASTQ or a FASTA. Reads and reference bases are not variants.

Frequently Asked Questions

How do I count variants in a VCF file with one command?

Run bcftools view -H -Ov file.vcf.gz | wc -l for the record count, or bcftools index -n file.vcf.gz if the file is indexed. The index route is instant and returns the same number.

Does wc -l work on a compressed VCF?

No. wc -l on a .vcf.gz counts newline bytes in the compressed stream and returns a meaningless number. Decompress first with zcat, or use bcftools view -H to stream records.

What is the difference between counting sites and counting alleles?

A site is one data line. An allele is one ALT entry. A record with ALT=G,T is one site and two alleles, so allele counts are always greater than or equal to site counts. Split multiallelics with bcftools norm -m-any to count alleles.

Why does bcftools count fewer variants than my VCF has lines?

The extra lines are header metadata beginning with #. bcftools view -H suppresses them. If the difference is large and the file is a gVCF, most remaining records may be non-variant reference blocks with ALT=<NON_REF>.

How do I count only PASS variants?

Use bcftools view -f PASS -H file.vcf.gz | wc -l. This excludes records filtered as LowQual or similar, and also excludes records with FILTER set to ., which were never evaluated.

How do I count variants per sample?

Expand genotypes with bcftools query -f '[%SAMPLE\t%GT\n]' file.vcf.gz, then filter out homozygous reference and missing genotypes before tallying by sample name. Alternatively loop over samples with bcftools view -s SAMPLE -c1 -H | wc -l.

Why does counting variants in a region return zero?

The most common cause is a chromosome naming mismatch, such as chr1 in the command against 1 in the file. The second most common cause is a missing or stale index, and the third is a coordinate convention error from feeding 0-based BED coordinates to a 1-based region query.

Is a FASTA or FASTQ file needed to count variants?

Neither is needed for counting. FASTA is useful for checking that REF alleles match the reference genome, and FASTQ is the raw read data that produced the calls. Counts come from the VCF or BCF alone.