This is a test version of Biostars. For the public version, visit https://www.biostars.org.
how to extract unique variants from GVCF

[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?

gatk variant-calling gvcf filtering

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 variants
  • final List<Genotype> L=V.getGenotypes().stream() : get a stream of genotpes for variant V
  • filter(G->G.isHet() || G.isHomVar()). get only genotypes HET or HOM_VAR
  • collect(Collectors.toList()); convert to list
  • if(L.size()!=1) return; list must contain a singleton
  • final Genotype g=L.get(0); the singleton
  • println(V.getContig()+" "+V.getStart()+" "+V.getReference()+" "+g.getSampleName()+" "+g.getAlleles()); print the singleton

Now this post needs to be added to the list there, I guess?

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.