# How to Extract Variants from a VCF File by Region 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. | Format | What it holds | Coordinates | Typical use here | |--|--|--|--| | VCF | Variant records plus genotypes and per-variant quality | 1-based, REF-anchored | The input and the usual output | | BCF | Binary, compressed equivalent of VCF | Same as VCF | Faster and smaller intermediate files | | BED | Named intervals, tab-delimited, no header required | 0-based, half-open | Regions file for `bcftools -R` or `bedtools -b` | | GFF3 / GTF | Gene and transcript annotation | 1-based, inclusive | Turning gene models into regions | | FASTA | Sequence only, with `>` headers | None | Reference genome for normalization or annotation | | FASTQ | Sequence plus per-base quality, `@` headers | None | Raw reads before alignment, never a region source | | BAM / CRAM | Aligned reads | 0-based internally, displayed 1-based | The layer beneath a VCF | ### 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. ```bash bcftools view -h cohort.vcf.gz | grep '^##contig' | head -5 ``` Expected output looks like this, with UCSC-style names: ``` ##contig= ##contig= ##contig= ``` 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: ```bash 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: ```bash 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: | Setting | Behavior | |--|--| | `0` | Record is included if its POS falls inside the region | | `1` | Record is included if any part of its REF span overlaps the region (default) | | `2` | Record is included only if the whole REF span fits inside the region | 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 ```bash 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 ``` ```bash 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 ``` ```bash 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: ```bash 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. ```bash 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 `` 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: ```bash bcftools view -r chr1:1000000-2000000 -e 'ALT=""' -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. ```bash 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 ```bash 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 ```bash 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. ```bash 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 ```bash 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 ```bash 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 ```bash 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 ```bash 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]. ```r 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 | Parameter | Tool | Effect | Watch out for | |--|--|--|--| | `-r` | bcftools | Region list on the command line, needs index | Chromosome name must match the header exactly | | `-R` | bcftools | Region list from a file | `.bed` suffix switches to 0-based half-open | | `-t` / `-T` | bcftools | Targets, streams without index, `^` for complement | Slower on large files | | `-regions-overlap` | bcftools | 0, 1, or 2 for POS-only, any-overlap, fully-contained | Default 1 includes spanning records | | `-f` | bcftools | Filter status, for example `PASS` | `PASS` alone drops unfiltered `.` records | | `-i` / `-e` | bcftools | Include or exclude by expression | Quote expressions in single quotes | | `-Oz` / `-Ov` / `-Ob` | bcftools | Output as compressed VCF, plain VCF, BCF | `-Ov` piped to tabix fails | | `-no-version` | bcftools | Suppresses provenance header lines | Removes reproducibility metadata | | `-L` | GATK | Interval in 1-based inclusive form | Requires the matching reference FASTA | | `-interval-padding` | GATK | Expands each interval by N bases | Padding can pull in neighboring genes | | `-from-bp` / `-to-bp` | vcftools | Single interval, 1-based inclusive | One chromosome per run | | `-h` | tabix | Include header with output | No filtering or reformatting | ## 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 `` 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: ```bash 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](/knowledge/bioinformatics/linkage-disequilibrium-and-haplotype-mapping) 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.](https://pubmed.ncbi.nlm.nih.gov/27401132/) 2. [Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR.](https://pubmed.ncbi.nlm.nih.gov/26379229/) 3. [Whole-genome resequencing of the Kazakh coarse-wool fat-tailed sheep identifies functional variants in disease-associated genes.](https://pubmed.ncbi.nlm.nih.gov/42840686/) 4. [sim1000G: a user-friendly genetic variant simulator in R for unrelated individuals and family-based designs.](https://pubmed.ncbi.nlm.nih.gov/30646839/) 5. [Novel genetic variants data for adaptation to hypoxia in native chickens.](https://pubmed.ncbi.nlm.nih.gov/37735456/)