Now this post needs to be added to the list there, I guess?
[note: cross-posted on GATK forum - still awaiting a response]
I have a GVCF (generated using GATK's HaplotypeCaller w/ -ERC GVCF parameter) of 36 related samples and would like to determine the (potentially de novo) variants that are unique to each sample. Short of creating 36 N-1 GVCFs for discordance testing, or extracting individual sample VCFs for subtraction, is there a straightforward method to obtain the desired information?
1 answer
merge your g.vcf files with CombineGVCFs https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_CombineGVCFs.php
and then I would use bioalcidaejdk to find the sample containing a singleton: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html
java -jar dist/bioalcidaejdk.jar -e 'stream().forEach(V->{final List<Genotype> L=V.getGenotypes().stream().filter(G->G.isHet() || G.isHomVar()).collect(Collectors.toList());if(L.size()!=1) return;final Genotype g=L.get(0);println(V.getContig()+" "+V.getStart()+" "+V.getReference()+" "+g.getSampleName()+" "+g.getAlleles());});' input.vcf
stream():get a stream of variantsfinal List<Genotype> L=V.getGenotypes().stream() :get a stream of genotpes for variant Vfilter(G->G.isHet() || G.isHomVar()).get only genotypes HET or HOM_VARcollect(Collectors.toList());convert to listif(L.size()!=1) return;list must contain a singletonfinal Genotype g=L.get(0);the singletonprintln(V.getContig()+" "+V.getStart()+" "+V.getReference()+" "+g.getSampleName()+" "+g.getAlleles());print the singleton
Will do. I cross-posted here only b/c the response time there can be... slow.
Sorry, what? I was saying that this post would be linked in Pierre's post page, which he has done now.
@Ram: Sorry for the confusion. I thought your comment was directed at me, since I mentioned cross-posting on GATK. I still intend to share Pierre's solution there, so others will be aware of it.
Ah, OK. No, I was referring to the cross-linking that Pierre does on his tool-page. Sorry, I should have been clearer in my comment.
As usual, @PierreLindenbaum has the right tool for the job. Thanks!
I added -F VCF to run the code while used for gvcf file. Is there any user manual to understand the syntax?
Log in to answer this question.