how to filter a multi-VCF for sex-specific genotypes?

Question

I have a single VCF file of variants for a population that includes males & females and I would like to pull out all variants that are sex-specific. For example, sites where all males are 0/1 and all females are 0/0. Is there a way to do this using bcftools or vcftools? I could enumerate each individual name, or ask for all individuals named *.M and *.F.

I first tried running smoove on the male and female populations separately then comparing the two VCFs, but I'm working with a draft genome with loads of contigs. As only contigs containing variants are output in each VCF, and these are different between the two files, vcftools couldn't handle the comparison.

Someone here https://github.com/samtools/bcftools/issues/118 managed to get it working ike this:

$bcftools view -f PASS -i 'GT~"0/0.*0/1.*0/1.*0/1.*0/1"' combined.vcf.gz

but the same syntax gives me no variants (yet I can see some that meet my criteria in IGV).

Answer 1

Found a great tool with SnpSift!

From the examples: I want to keep samples where the genotype for the first sample is homozygous variant and the genotype for the second sample is reference:

cat variants.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] ) & isVariant( GEN[0] ) & isRef( GEN[1] )" > filtered.vcf 

Answer 2

if you could provide seperate list of male and female samples names you can something like this:

$ bcftools view -e 'GT[@males.txt]!="0/1" && GT[@females.txt]!="0/0"' combined.vcf.gz

or if you don't want to take care of the exact genotype you can write it more generic:

$ bcftools view -e 'GT[@males.txt]!="alt" && GT[@females.txt]!="ref"' combined.vcf.gz

fin swimmer

Answer 3

Given an input file input.vcf, I would take this approach:

1. Create two text files: one with the male sample names (males.txt) and one with the female sample names (females.txt) (one sample per line)

2. Create a BCF (for speed) for all the samples in each list while dropping the INFO field so that INFO/AC is recalculated (I'm not sure dropping the INFO field is strictly necessary). Then only keep sites with a minor allele count greater than zero:

   bcftools view --samples-file males.txt input.vcf -Ou | bcftools annotate -x INFO  -Ou | bcftools view --min-ac 1:minor -Ob -o male_variants.bcf
   bcftools view --samples-file females.txt input.vcf -Ou | bcftools annotate -x INFO -Ou | bcftools view --min-ac 1:minor -Ob -o female_variants.bcf
   

3. Use isec to create intersects and complements of the two variant files and store them in dir/*:

   bcftools isec male_variants.bcf female_variants.bcf -p dir -Oz
   

Double-check dir/README.txt, but the files you seek should be named dir/0000.vcf.gz and dir/0001.vcf.gz

Disclaimer: This code is untested.

Answer 4

using vcffilterjdk http://lindenb.github.io/jvarkit/VcfFilterJdk.html

or ask for all individuals named *.M and *.F.

java -jar dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().allMatch(G->(G.getSampleName().endsWith("M") && G.isHet()) || (G.getSampleName().endsWith("F") && G.isHomRef())); ' input.vcf

variant.getGenotypes().stream() get a stream of genotypes 'G' for the current variant.

allMatch all genotypes must match the statements below:

G.getSampleName().endsWith("M") && G.isHet() sample 's name ends with 'M' and is heterozygous

OR

G.getSampleName().endsWith("F") && G.isHet() sample 's name ends with 'F' and is Hom-ref