How to Extract Variants from a VCF File by Region

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

By the end of this guide you will be able to pull any genomic interval out of a VCF file, keep the header and sample columns intact, prove the output is correct, and sidestep the coordinate and chromosome-naming traps that make region queries return empty files with no warning. You will be able to do it with a single command for one interval, with a regions file for hundreds, and with a filter expression attached so that only the variants you want come out.

You need four things on hand before you start. A sorted VCF that has been compressed with bgzip, an index for that file (.tbi or .csi), the exact contig names as they appear in the VCF header, and one of the tools covered below. Install bcftools and the htslib utilities, which give you bcftools, bgzip, and tabix in one package. Optionally install GATK, bedtools, or R with the vcfR package if you prefer a different route. A reference FASTA is not needed to slice a region out of an existing VCF. It is needed only if you plan to normalize, left-align, or annotate the result.

The Mental Model: What a Region Query Really Does

A VCF (variant call format) is a tab-delimited text file. The top of the file holds metadata lines that begin with ##. A single column header line begins with #CHROM and lists CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, and then one column per sample. Every following line describes one variant: where it is, what the reference allele is, what the alternate allele is, how confident the caller was, and the genotype of each sample.

A region query is not a text search through those lines. It is a jump. The file is compressed in small blocks with bgzip (a blocked version of gzip), and a companion index records which genomic interval lives in which block. With that index, a tool can seek to bytes near your interval and read a few blocks instead of streaming the whole file. That is why a region query on a 40 GB VCF takes seconds.

Two facts follow from this and they govern everything else in this article.

First, without an index, bcftools view -r cannot work. It will report that it could not retrieve the index file and will stop. The workaround is -t (targets), which streams the whole file and is far slower, or tabix -p vcf to build the index once.

Second, coordinates must be interpreted consistently. VCF POS is 1-based and refers to the first base of the REF allele. BED files are 0-based and half-open, so a BED interval of chr1 999999 2000000 covers the same bases as a VCF-style interval of chr1:1000000-2000000. Getting this wrong shifts your window by exactly one base at each edge, which rarely breaks a pipeline but silently changes which boundary variants are included. Many bug reports about "missing variants" are boundary bugs.

Before You Start: Files, Indexes, and Naming

The bioinformatics file formats you will meet

Region work touches several formats at once. Knowing which one carries quality scores, which one is coordinate-based, and which one is just sequence prevents most confusion.

FASTQ versus FASTA, and why neither is a region container

The shortest useful distinction is this. FASTA stores sequence. FASTQ stores sequence plus a quality score for every base.

A FASTA record is a header line beginning with > followed by sequence lines that may wrap at any width. A FASTQ record is exactly four lines: a header beginning with @, the sequence, a separator line beginning with +, and a quality string of the same length as the sequence, encoded as ASCII characters that map to Phred scores offset by 33. A base with quality I in FASTQ has a Phred score of 40, meaning a 1 in 10,000 chance the base call is wrong. FASTA has no such information, which is why a reference genome is distributed as FASTA and raw sequencing output arrives as FASTQ.

You cannot extract variants from FASTQ by region. Reads in a FASTQ file have no genomic position until they are aligned to a reference, and the variant calls that follow alignment are what land in a VCF. Region extraction is an operation on positioned data. If your starting point is FASTQ, the path is FASTQ to BAM or CRAM through an aligner, then variant calling, then this article.

Step 0: confirm contig names and build the index

Never guess the contig names. Read them from the header.

bcftools view -h cohort.vcf.gz | grep '^##contig' | head -5

Expected output looks like this, with UCSC-style names:

##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr3,length=198295559>

If your file instead shows ID=1, ID=2, ID=3, then every region string in this article must be written as 1:1000000-2000000, not chr1:1000000-2000000. A mismatch between chr1 and 1 is the single most common reason a region query returns nothing.

To compress and index a plain VCF:

bgzip -c raw.vcf > cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

The tabix -p vcf command builds cohort.vcf.gz.tbi. The bcftools equivalent is bcftools index cohort.vcf.gz, which builds a .csi index by default, and bcftools index -t cohort.vcf.gz forces a .tbi index. Prefer .csi when any contig is longer than 512 megabases, which does happen with some plant and animal assemblies. Both index types work with bcftools view -r.

Check the index reports sensible counts:

bcftools index --stats cohort.vcf.gz

chr1	248956422	214398
chr2	242193529	198844
chr3	198295559	176102

The columns are contig name, contig length, and number of records. If the record counts look far too low, your file may be a gVCF (which stores reference blocks) or may have been filtered already.

Step-by-Step: Extract Variants from a VCF by Region

The walkthrough below uses bcftools because it is fast, keeps the header, and handles every awkward case. The alternatives section shows the same job in other tools.

Step 1: Decide what "in the region" means at the edges

Before typing a command, decide how a record that only partly overlaps the interval should be treated. bcftools exposes this through -regions-overlap with three settings:

The default of 1 is usually right. A 20-base deletion that starts at position 999,990 and extends to 1,000,009 overlaps a region starting at 1,000,000, so it belongs in the output. If you are counting variants per window and want each record assigned to exactly one window, 0 gives you that, at the cost of dropping spanning deletions from every window they cross except the one containing their start.

Step 2: Run the single-region extraction

bcftools view -r chr1:1000000-2000000 -Oz -o chr1_region.vcf.gz cohort.vcf.gz

The flags do specific work. -r takes a comma-separated list of regions in 1-based inclusive form and requires the index. -Oz writes compressed VCF output (-Ov writes plain text, -Ob writes BCF). -o names the output file. If you omit -o and redirect with >, the command still works, but you lose the automatic format detection that comes with -o.

Expected output on a well-formed run is silence, and a new file:

$ ls -lh chr1_region.vcf.gz
-rw-r--r--  1 user  staff   1.9M Nov 10 09:41 chr1_region.vcf.gz

Every bcftools command writes provenance into the header by default, in lines such as ##bcftools_viewVersion and ##bcftools_viewCommand. That is good practice for reproducibility. Add -no-version if you need byte-identical output across machines for a test fixture.

Step 3: Extract many regions at once

For more than a handful of intervals, put them in a file and use -R.

chr1	1000000	2000000
chr17	41196312	41277500
chrX	3100000	3200000

bcftools view -R regions.txt -Oz -o panel_regions.vcf.gz cohort.vcf.gz

The tab-delimited file above is read with 1-based inclusive coordinates, which matches VCF. If you save the same intervals as a BED file, the numbers change because BED is 0-based half-open:

chr1	999999	2000000
chr17	41196311	41277500
chrX	3099999	3200000

bcftools view -R regions.bed -Oz -o panel_regions.vcf.gz cohort.vcf.gz

bcftools detects the BED convention from the .bed suffix and converts internally. This is a real trap. A file of 0-based coordinates named regions.txt will be read as 1-based and every interval will start one base too late. Name BED files with a .bed extension, or convert them, and never rely on memory alone.

You can also pass regions directly on the command line, comma separated, and you can use the caret to mean the complement:

bcftools view -r chr1:1000000-2000000,chr1:5000000-5100000 -Oz -o two_windows.vcf.gz cohort.vcf.gz
bcftools view -T ^regions.bed -Oz -o everything_else.vcf.gz cohort.vcf.gz

-T is the targets option. It streams the file, does not need an index, and takes the complement when prefixed with ^. It is the right choice for one-off filters on a file you have not indexed, and the wrong choice for repeated queries.

Step 4: Filter while you extract

Region and filter combine in one pass, which saves writing an intermediate file.

bcftools view \
  -r chr1:1000000-2000000 \
  -f PASS \
  -i 'INFO/DP>=10' \
  -Oz -o chr1_region_qc.vcf.gz cohort.vcf.gz

-f PASS keeps only records marked as passing filters. -i takes an expression, and -e excludes on an expression. INFO/DP>=10 requires a minimum total depth of 10 reads at the site. Other useful expressions include INFO/AF>=0.05 for allele frequency and TYPE="snp" to keep only single-nucleotide variants.

If you are working from a gVCF produced by a reference-model caller, the file contains <NON_REF> records that describe blocks of reference-matching sequence rather than variants. Extract them if you want the reference model preserved, exclude them if you want variant sites only:

bcftools view -r chr1:1000000-2000000 -e 'ALT="<NON_REF>"' -Oz -o chr1_variants.vcf.gz cohort.g.vcf.gz

Step 5: Index the output if you plan to query it again

A freshly written .vcf.gz has no index.

bcftools index chr1_region.vcf.gz

Skip this and the next person who runs bcftools view -r on your output will hit the missing-index error. If you plan to use the file with IGV or another browser, the .tbi index is what makes scrolling instant.

Step 6: Verify the record count two ways

bcftools index --stats chr1_region.vcf.gz
bcftools view -H chr1_region.vcf.gz | wc -l

chr1	248956422	18422
18422

The two numbers should agree. bcftools index, stats reads counts from the index built over the file, and bcftools view -H | wc -l counts records directly. A disagreement means the index is stale, which happens when a file is overwritten in place without reindexing. Delete the .csi or .tbi and rebuild.

Step 7: Confirm the boundaries

bcftools query -f '%CHROM\t%POS\n' chr1_region.vcf.gz | sort -k2,2n | head -1
bcftools query -f '%CHROM\t%POS\n' chr1_region.vcf.gz | sort -k2,2n | tail -1

chr1	1000001
chr1	1999998

The first position does not have to equal 1,000,000. It is simply the first variant in the file at or after the boundary. What matters is that no position falls outside the interval and that the position you expected to be present, if you checked one by hand, is there. With the default overlap setting of 1, a record with POS slightly below 1,000,000 can legitimately appear because its REF span reaches into the region.

Step 8: Hand the region to the next step

From here the region file is a normal VCF. You can annotate it, split multi-allelic records, subset samples, or convert formats.

bcftools norm -m -any -f ref.fa -Oz -o chr1_region_split.vcf.gz chr1_region.vcf.gz
bcftools view -s sampleA,sampleB -Oz -o chr1_region_two_samples.vcf.gz chr1_region.vcf.gz

bcftools norm -m -any splits multi-allelic records into one record per ALT allele and -f ref.fa left-aligns and normalizes indels against the reference. Do this deliberately. Splitting changes the record count, so re-run the count check afterward and do not expect the numbers to match the pre-split file.

Alternative Tools and What They Change

bcftools is not the only way, and sometimes the choice is dictated by the rest of your pipeline.

tabix directly

tabix -h cohort.vcf.gz chr1:1000000-2000000 > chr1_region.vcf

tabix prints raw text lines and does not rewrite the header. -h includes the header block. The coordinates in the region argument are 1-based inclusive, matching VCF. This is the leanest option when you want text and no bcftools-specific header additions. It does not filter, split, or reformat anything.

vcftools

vcftools --gzvcf cohort.vcf.gz --chr chr1 --from-bp 1000000 --to-bp 2000000 \
  --recode --recode-INFO-all --stdout > chr1_region.vcf

vcftools takes one chromosome and one interval per run, which makes it awkward for panels but simple for scripts that already use it. -recode-INFO-all preserves the INFO field, which is dropped by default in older versions unless requested. The coordinate arguments are 1-based inclusive. vcftools streams the file and does not require an index, so it is slower on large inputs but tolerant of unindexed files.

GATK SelectVariants

gatk SelectVariants \
  -R ref.fa \
  -V cohort.vcf.gz \
  -L chr1:1000000-2000000 \
  -O chr1_region.vcf.gz

GATK uses -L for an interval and accepts either a command-line interval, an interval list file, or a BED file. -interval-padding 100 extends each interval by 100 bases on both sides, which is useful for capture or exome work where reads near the bait edge are undercalled. -exclude-intervals removes intervals. GATK requires the reference FASTA here, matching it to the VCF, so this route is heavier than bcftools but integrates with the rest of a GATK pipeline. Keep BED files in BED convention and confirm the first record's position in the output rather than assuming.

bedtools intersect

bedtools intersect -header -wa -a cohort.vcf.gz -b regions.bed > chr1_region.vcf

-wa writes the entire A record for each overlap, and -header carries the VCF header through. bedtools works on genomic intervals in general, not just variants, so it is the natural choice when your regions come from another interval file such as peaks or exons. The VCF passed to -a should be bgzipped and tabix-indexed, and if it is not, decompress it to a stream first.

R with vcfR

The vcfR package was built precisely because many R readers discarded the per-variant quality data that sits next to each genotype, and it provides functions to read, write, and extract portions of a VCF along with plots of summary statistics [1].

library(vcfR)

vcf <- read.vcfR("cohort.vcf.gz", verbose = FALSE)
fix <- getFIX(vcf)
pos <- as.integer(fix[, "POS"])

keep <- fix[, "CHROM"] == "chr1" & pos >= 1000000 & pos <= 2000000
sub  <- vcf[keep, ]

write.vcf(sub, file = "chr1_region.vcf.gz")

This approach reads the whole file into memory before subsetting. That is fine for a few hundred thousand variants and painful for tens of millions. For large files, run bcftools first and read the small region file into R. The limit argument to read.vcfR can cap how many lines are read, which is handy for inspecting the top of an enormous file without loading it all.

Parameters That Change Your Answer

Common Mistakes and Limitations

Chromosome naming mismatch

The VCF header says chr1 and you typed 1. The query usually returns zero records and exits successfully. Nothing in the console tells you the region was never found, which is why this error survives into downstream analysis. Always read the ##contig lines first. If you need to convert, use a two-column mapping file and bcftools annotate, rename-chrs map.txt, then rebuild the index. Renaming contigs without reindexing leaves an index that points at the wrong blocks.

Mixing 1-based and 0-based coordinates

VCF is 1-based. BED is 0-based half-open. GFF and GTF are 1-based inclusive, which matches VCF on the start but not on how features are described. Convert deliberately, name BED files with .bed, and spot-check one known variant at the left edge of a region by searching for its position in the output.

Running a region query without an index

The error looks like this:

[E::idx_find_and_load] Could not retrieve index file for 'cohort.vcf.gz'
Failed to read from cohort.vcf.gz: could not load .tbi/.csi index

The fix is tabix -p vcf cohort.vcf.gz or bcftools index cohort.vcf.gz. If you cannot index because the file is not sorted, sort it first with bcftools sort or the equivalent, then compress and index.

Indexing a file compressed with plain gzip

tabix -p vcf fails on ordinary gzip with a message about the file not being BGZF:

tbx_index_build failed: not a BGZF file

Recompress with bgzip -c raw.vcf > cohort.vcf.gz. The two formats share a .gz extension and are not interchangeable for random access.

Confusing -r with -t under time pressure

-r jumps using the index and is fast. -t scans the whole file and works without one. On a 40 GB VCF the difference is minutes versus seconds. Using -t because you forgot to index is a habit that gets expensive at scale, and the two options also have separate overlap settings, so a switch between them can change which records appear.

Expecting -s to filter variants

bcftools view -s sampleA subsets the sample columns. It removes samples, not variants. A site that is homozygous reference in the remaining samples is still written unless you add an expression such as -i 'GT="alt"' or an allele count filter with -c.

Forgetting that multi-allelic records are single records

A record with ALT A,C is one line in the VCF but two distinct variants. Region counting, allele frequency filters, and most annotation tools behave differently on split and unsplit records. Decide before you start and run bcftools norm -m -any if you need one ALT per line. Split first if you intend to filter on a per-allele value, because INFO/AF on a multi-allelic record reports both values in a comma-separated list.

Losing the header

Piping bcftools view -H into a file, or using grep -v '^#', removes the ##contig and ##INFO lines. The result is a variant table that most tools refuse to read, with errors about a missing header or unknown contig. Region extraction preserves the header. Deliberate post-processing that strips it should be the exception.

gVCF reference blocks treated as variants

A gVCF contains <NON_REF> records covering blocks of reference-matching positions. A region extraction from a gVCF can return thousands of records that are not variant sites at all. Filter them out unless you specifically need the reference model.

Stale indexes after in-place edits

Editing a bgzipped VCF with a tool that writes in place, or replacing it with a new file of the same name, leaves the old .tbi behind. Queries then return wrong or truncated results. Rebuild the index after any change to the file.

Assuming a small output means no variants

A region may genuinely contain no called variants, particularly in a low-coverage sample or a gene-poor interval. Before you accept that answer, confirm the contig name, confirm the build of the coordinates matches the VCF, and confirm the boundaries are inside the contig length listed in the header. Coordinates from a different genome build land in the wrong place and usually land in a region with different content.

How to Check the Result Is Right

Run three cheap checks on every extraction.

First, count. Compare bcftools index, stats on the region file with bcftools view -H | wc -l. Disagreement points to a stale index.

Second, bracket the positions. The minimum and maximum POS in the output must fall inside the requested interval, allowing one edge record below the start when the default overlap setting includes spanning deletions.

Third, confirm the header survived. bcftools view -h chr1_region.vcf.gz | grep '^##contig' should match the source file. If you renamed contigs, this is where the mistake shows.

For a stricter check, extract the region a second way, for instance with tabix, and compare sorted records:

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' chr1_region.vcf.gz | sort > a.txt
tabix -h cohort.vcf.gz chr1:1000000-2000000 | grep -v '^#' | cut -f1,2,4,5 | sort > b.txt
diff a.txt b.txt

An empty diff means the two methods agree. Any difference is worth a look before you build on the file.

Sensible Next Steps

Region files are usually intermediates. The most common next move is annotation, where each variant is assigned to genes, regulatory elements, and population frequencies. ANNOVAR provides gene-based, region-based, and filter-based annotation directly on a VCF, and its protocol runs in roughly 5 to 30 minutes of computer time for a typical variant file [2]. Annotation depends on build and contig naming matching between the VCF and the annotation database, which is another reason to fix chromosome names before you leave this step.

A second common use is targeted locus work. In a whole-genome resequencing study of Kazakh coarse-wool fat-tailed sheep, variants within the PRNP and TMEM154 loci were extracted from filtered VCF files with bcftools before downstream analysis [3]. That pattern, call variants across the genome, filter, then pull specific loci for focused inspection, is the standard way to move from a genome-wide call set to a tractable candidate region.

A third use is simulation and linkage disequilibrium work. The sim1000G package takes a raw phased VCF as its only required input and pulls haplotypes from genomic regions to compute LD structure for simulated datasets [4]. If you plan to simulate from real haplotypes, the region file must retain phase and the full sample columns.

Finally, filtered and extracted variant sets are routinely shared as VCF. A study of hypoxia adaptation in native chickens identified 3877 novel single-nucleotide variants and deposited them in VCF format at EBI for reuse [5]. Keeping your region files well named, indexed, and documented with the source file and command in the header makes that reuse possible.

Frequently Asked Questions

Can I extract variants from a VCF by region without an index?

Yes, with bcftools view -t or by streaming the file, but it is slow because the tool must read every record. Building an index once with tabix -p vcf or bcftools index makes every later query fast.

Why does my region query return an empty file with no error?

The chromosome name in your region string does not match the name in the VCF header, or the coordinates come from a different genome build. Check the ##contig lines before anything else.

Should I use bcftools -r or -t?

Use -r when the file is indexed and you will query it more than once. Use -t when the file has no index or you want the complement of a region with the ^ prefix.

How do I extract several regions at once?

Put the intervals in a tab-delimited file or a BED file and pass it with bcftools view -R regions.txt. The .bed suffix tells bcftools to read 0-based half-open coordinates, and a plain text file is read as 1-based inclusive.

Can I extract variants from a FASTQ file by region?

No. FASTQ reads have no genomic coordinates until they are aligned. Align to a reference to produce BAM or CRAM, call variants, then extract regions from the resulting VCF.

What is the difference between a VCF and a BCF?

They hold the same information. VCF is tab-delimited text and BCF is its binary, compressed counterpart. BCF is smaller and faster for bcftools to read, and both use the same 1-based coordinates.

Does region extraction keep the header?

bcftools and GATK keep all header lines, including ##contig and ##INFO definitions, and bcftools adds provenance lines recording the command. Tools that print raw lines, such as tabix without -h, drop the header unless you ask for it.

How do I extract a region from a VCF in R?

Read the file with read.vcfR, pull the fixed columns with getFIX, build a logical index on CHROM and POS, and subset the object. For large files, extract with bcftools first and read only the small region file into R.

Sources

1. vcfr: a package to manipulate and visualize variant call format data in R.
2. Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR.
3. Whole-genome resequencing of the Kazakh coarse-wool fat-tailed sheep identifies functional variants in disease-associated genes.
4. sim1000G: a user-friendly genetic variant simulator in R for unrelated individuals and family-based designs.
5. Novel genetic variants data for adaptation to hypoxia in native chickens.