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?
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.
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 intersectallows 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 intersectworks 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