3
$\begingroup$

I have two VCF files and I want to compare the missing rate in each of them, but I want to only look at sites that are present in both files. How could I go about getting a list of positions that are included in both VCF files and then filtering those VCF files to only include those positions?

$\endgroup$

3 Answers 3

2
$\begingroup$

Something like vcf-isec? Quoting:

Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.

$\endgroup$
1
  • $\begingroup$ Thank you!!! I had suspected that there would be something like that but I was having trouble finding anything. I really appreciate it $\endgroup$
    – Sarah
    May 31, 2019 at 8:36
3
$\begingroup$

The bedtools intersect command is my go-to for operations like this. From the manual:

By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. bedtools intersect allows one to screen for overlaps between two sets of genomic features. Moreover, it allows one to have fine control as to how the intersections are reported. bedtools intersect works with both BED/GFF/VCF and BAM files as input.

$\endgroup$
2
$\begingroup$

You can try something like this in bcftools, which is a bit quicker than VCFtools, as it uses block gzipping and indexing.

# make sure both files are indexed and in bgzip format 
bcftools index vcfA.vcf.gz
bcftools index vcfB.vcf.gz

# find overlapping positions between the two files
bcftools isec -n~11 -c all vcfA.vcf.gz vcfB.vcf.gz > sharedPositions.txt

# filter vcf files for the overlapping positions
bcftools view -T sharedPositions.txt vcfA.vcf.gz -Oz > vcfA.sharedPositions.vcf.gz

bcftools view -T sharedPositions.txt vcfB.vcf.gz -Oz > vcfB.sharedPositions.vcf.gz

The perhaps use something like plink2 to look at the missing rate per variant and individual for each file:

plink2 --vcf vcfA.vcf.gz --missing --out vcfA.missingStats
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.