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

Question

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?

Answer 1

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

Answer 2

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.

Answer 3

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:

https://gatk.broadinstitute.org/hc/en-us/articles/360041849831-SelectVariants

For this case (subsetting by chromosome) you would simply do:

gatk SelectVariants \
     -R reference.fasta \
     -V input.vcf \
     -L chr2 \
     -O output.chr2.vcf

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:

https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists