Extract Sub-Set Of Regions From Vcf File
6
13
Entering edit mode
9.2 years ago
Rubal7 ▴ 800

Hello all,

Sorry if this has been answered before (though I have not found an appropriate answer already) What is an effective way to extract mutliple chromosomal regions from a vcf.gz file? I see that it can be done for single regions using tabix, but I have a file with several regions in a text file in the format:

chromosome:startposition-endposition eg: 19:47366525-47380539

Perhaps there is a way to refer to a file of regions using tabix, or vcf-tools? If anyone could provide an example command line that would be really helpful.

vcf tabix genome filter • 53k views
0
Entering edit mode
$echo -e "chrN\t123\t456" | bedops -e 1 <(vcf2bed variants.bed) - > answer.bed  ADD REPLY 24 Entering edit mode 9.2 years ago I find tabix to be much easier / much sorter command line arguments. You don't need a fasta for this method. Tabix is screaming fast. bgzip your.vcf tabix -p vcf your.vcf tabix your.vcf.gz chr1:10,000,000-20,000,000  ADD COMMENT 0 Entering edit mode Thanks, I agree it works fast. Shame it doesn't take a file with a whole list of regions. But I can put it into a loop. ADD REPLY 4 Entering edit mode Two ways: 1) if there are many regions: tabix -fB my.vcf.gz reg.bed 2) if you only have a handful of regions: awk '{print$1":"($2+1)"-"$3}' reg.bed | xargs tabix my.vcf.gz.

7
Entering edit mode

For anyone arriving years later, the command appears to have changed. The command you probably want is: tabix -R reg.bed my.vcf.gz.

0
Entering edit mode

None of these work anymore.

6
Entering edit mode
tabix -R region.txt my.vcf.gz


worked for me

1
Entering edit mode

Would anyone be able to let me know what the input is for the region.txt file? I have put 57646425-57803821 (all in one line). Linux is telling me could not parse bed line. Thank you for your help.

0
Entering edit mode

Have you tried to add the chromosome position in the beginning? eg: chr1:57646425-57803821

Another way is separate by tab. eg: chr1 57646425 57803821

0
Entering edit mode

Hello geocarvalho, thank you for the suggestion. I'll give that a try. Really appreciate your help. Julia

0
Entering edit mode

Thanks so much this worked :)

18
Entering edit mode
7.7 years ago
gourneau ▴ 210

For others if you like to use BED files, one way to do this would be like so:

vcftools --vcf input.vcf --bed bed_file_describing_the_range.bed --out output_prefix --recode --keep-INFO-all


Where bed_file_describing_the_range.bed has the positions of interest.

0
Entering edit mode

I am trying to extract specific regions from a VCF file that has SNP information for 16 individuals. My bed file has three columns (scaffold name, start of region, end of region).

When I run vcftools on these files, it runs through the entire vcf file and determines that there are no regions from the bed file contained in the vcf file. A recode vcf file is output, but it only contains the header information from the vcf file.

0
Entering edit mode

Check to ensure your scaffold names match in the BED and VCF files. For example, if your BED file has "chr1" but your VCF file just has "1" for chromosome 1 it will return an empty file.

0
Entering edit mode

i met the same problem with u , so I am wondering how did you solved this problem?

1
Entering edit mode

Answering to help others you are going to need to use sed on the bed file ie sed 's/chr//g' {in_path} > {out_path}

2
Entering edit mode
3.8 years ago
geocarvalho ▴ 180

bedtools intersect worked better to me

intersectBed -a input.vcf -b /path/to/my.interval.bed -header > output.vcf

0
Entering edit mode

For some reason by intersected out VCF have more line compare to input vcf. :( :( No idea what could be the possible reason for this !! Any guess ?

1
Entering edit mode
9.2 years ago
Johan ▴ 870

You could use the GATK SelectVariants Walker, specifying the -L option to set the intervals your are interested in.

Here is an example of how to do that, taken from the link above:

#Select a sample and restrict the output vcf to a set of intervals:
java -Xmx2g -jar GenomeAnalysisTK.jar \
-R ref.fasta \
-T SelectVariants \
--variant input.vcf \
-o output.vcf \
-L /path/to/my.interval_list \
-sn SAMPLE_1_ACTG


Note that you have to have a fasta reference file with the same names and coordinates as you use in our region file and vcf, and that you can remove the -sn option if you want to extract the variations for all samples.

Hope this helps.

0
Entering edit mode

thanks, I can potentially use this, but having some issues with the reference fasta, so a method that does not require a fasta reference would be preferred.

0
Entering edit mode

I'm not sure if it requires uncompressing the vcf files. I'm not that familiar with vcftools, so I couldn't say if there is an easier way to do it using that suite.

1
Entering edit mode
7.7 years ago

I did something like what you describe, but to create multiple VCF files. Uses bedtools: