exome/vcf merge question

[jhersheson]

I research rare mendelian diseases and generally look for shared variants between my related affected samples. Because of variation in the quality and coverage of exomes, I want to be able to look at a merged variant file for all of my affected cases.

This is relatively straight forward with vcf tools however if for example I am looking at the shared variants of two individuals, some variants may not be shared because of three reasons: 1. one individual has variant allele and one has wild-type i.e. not shared. 2) variant not covered on second exome 3) variant covered on second exome but at allele frequency below cutoff for variant call.

Can anyone help with a method to annotate the merged variant list with the depth and allele call for each of the samples including calling the wild-type allele where true.

p.s I had thought of using bedtools to annotate the depth of read at each start position but this wouldnt help me annotate with wildtype calls.

Thanks

Josh

[kgulukota]

I have a similar problem. I get my sequence data in batches (not all samples at once) and would like to have a running list of variants called on samples thus far.

As you said, the biggest issue with straightforward merging of VCFs is that we need to differentiate between

• evidence of absence ("there is sufficient depth at this locus and this sample is reference homozygous") and
• absence of evidence ("this sample does not have enough coverage to infer whether there is a variant at this locus").

I am still searching for solutions and will post if I find one.

[kgulukota]

Re: Create a VCF with your first bam file, say 1.vcf

OK. There is a 3-step procedure that can accomplish what you want (I think).

Step 1. Create VCF's with your first and second bam files separately, say old.vcf and new.vcf.

Step 2. Next create a combined vcf with the two. I used the CombineVariants walker in GATK like so:

java -jar GenomeAnalysisTK.jar -T CombineVariants -R GRCh37.fa --variant old.vcf --variant new.vcf -o joined.vcf -genotypeMergeOptions  UNIQUIFY 

But presumably you can do similar with bedtools.

Step 3. Finally, run the GATK UnifiedGenotyper by using the joined vcf as the target file i.e. with the -L option, like so:

java -jar GenomeAnalysisTK.jar  -T UnifiedGenotyper -R GRCh37.fa -L joined.vcf -I old.bam -I new.bam -o final.vcf 

I have combined 30 old bams with 50 new bams using this method and seems to work well.

However, allow me to hasten to add that the best practice would be to run variant calling on all samples together. The above procedure might be quick and dirty. I think it will be mostly accurate but there will be differences between this procedure and redoing the whole shebang.