# How to subset a VCF by chromosome and keep the header?

I would like to subset a VCF which only has chromosome 2.

The problem with using various grep commands, e.g.

grep -w '^#\|^2' my.vcf > my_new.vcf


or if there's a 'chr' prefix

grep -w '^#\|^chr2' my.vcf > my_new.vcf


is that this will remove the header.

I've been trying to use vcftools with the following command:

vcftools --vcf input.vcf  --chr 2 --out subset


which does not output a VCF subset.vcf as expected, but rather subset.log:

VCFtools - 0.1.15
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
--vcf input.vcf
--chr 2
--out subset

After filtering, kept 2 out of 2 Individuals
After filtering, kept 80 out of a possible 720 Sites
Run Time = 0.00 seconds


Given that the run time is Run Time = 0.00 seconds, I guess there's an error here.

How can I subset a VCF and keep the header? Is there an option with bcftools instead maybe?

This seems to work for me with bcftools filter and the -r or -R argument.

   -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]


Comma-separated list of regions, see also -R, --regions-file. Note that -r cannot be used in combination with -R.

   -R, --regions-file FILE


Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the tab-delimited file are: CHROM, POS, and, optionally, POS_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while bgzip-compressed and tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note that chromosome ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE. However, within chromosomes, the VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that overlapping regions in FILE can result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be used in combination with -r.

Here's an example:

$bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000 ##fileformat=VCFv4.2 ##FILTER=<ID=PASS,Description="All filters passed"> ##fileDate=20171222 ##source=PLINKv1.90 ##contig=<ID=1,length=249167691> ##contig=<ID=2,length=242695901> ##contig=<ID=3,length=197800245> ##contig=<ID=4,length=190915651> ##contig=<ID=5,length=180666277> ##contig=<ID=6,length=170877445> ##contig=<ID=7,length=159086805> ##contig=<ID=8,length=146293415> ##contig=<ID=9,length=141018424> ##contig=<ID=10,length=135434552> ##contig=<ID=11,length=134938471> ##contig=<ID=12,length=133763353> ##contig=<ID=13,length=115045730> ##contig=<ID=14,length=107285438> ##contig=<ID=15,length=102369712> ##contig=<ID=16,length=90141356> ##contig=<ID=17,length=81006630> ##contig=<ID=18,length=78014583> ##contig=<ID=19,length=59071322> ##contig=<ID=20,length=62906515> ##contig=<ID=21,length=48077813> ##contig=<ID=22,length=51156934> ##contig=<ID=23,length=154847490> ##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746 ##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 1_125 3_5555 4_7565 5_9 4 71566 rs13125929 T C . PASS PR GT 0/1 0/1 0/0 4 87113 rs2006748 T C . PASS PR GT 0/1 0/0 0/0 4 110646 rs11727494 C T . PASS PR GT 0/1 0/0 1/1 4 142550 rs11735742 T C . PASS PR GT 0/1 0/0 0/0 4 200631 rs6826124 A G . PASS PR GT 0/1 1/1 0/1 4 221623 rs7695945 G T . PASS PR GT 0/0 0/0 0/0  • What is the exact command you're using with bcftools? – EB2127 Jan 25 '18 at 22:27 • as in the command line: bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4, which produces results from chromosome 4. – gringer Jan 26 '18 at 4:43 You can use grep for this. The problem is the -w flag which means "match entire words", but VCF header lines start with ##, not #. Here are some alternative grep commands: grep -E '^(#|2[[:space:]])' my.vcf  Or, if you have GNU grep (default on Linux): grep -P '^(#|2\t)' my.vcf  To match either 2 or chr2: grep -E '^(#|(chr2|2)[[:space:]])' my.vcf  or grep -P '^(#|(chr){0,1}2\t)' my.vcf  Alternatively, you can use awk: awk '/^#/ ||$1=="2" || \$2=="chr2"' my.vcf


Finally, if you want to use vcftools, you need to also use the --recode option:

--recode

--recode-bcf

These options are used to generate a new file in either VCF or BCF from the input VCF or BCF file after applying the filtering options specified by the user. The output file has the suffix ".recode.vcf" or ".recode.bcf". By default, the INFO fields are removed from the output file, as the INFO values may be invalidated by the recoding

(e.g. the total depth may need to be recalculated if individuals are removed). This behavior may be overriden by the following options. By default, BCF files are written out as BGZF compressed files.

But that will remove the INFO fields as explained above, so just use normal text-parsing tools like the ones I suggested.

• grep (at least the version on Linux) can use multiple -e commands as expressions to search for; lines matching any expression will be displayed. – gringer Jan 25 '18 at 20:25
• @gringer yes, that's POSIX. Personally, I prefer ERE or PCRE syntax over multiple -e, but both work. – terdon Jan 25 '18 at 20:45
• These are both great answers. I will probably mark @gringer as the accepted as there are less opportunities for readers to make mistakes, but personally, I'm using these grep commands. Thank you! – EB2127 Jan 26 '18 at 17:12
• @EB2127 you're very welcome, and please feel very free to mark whichever you prefer. I personally prefer rolling my own parser scripts, but with formats as complex as vcf, it is indeed usually better to use dedicated tools. I would go for gringer's answer as well in your place :) – terdon Jan 26 '18 at 17:33

You can do this more easily with the GATK tool SelectVariants, which allows you to subset variants based on various criteria. Here's the doc page:

gatk SelectVariants \

If you wanted to do this for several chromosomes, or regions within chromosomes, you would give the tool a file of intervals with -L instead of the chromosome name. Suitable interval file formats are described here: