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?
3 Answers
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.
-
$\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$– SarahMay 31, 2019 at 8:36
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.
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