Question: how to extract unique variants from GVCF
0
gravatar for harold.smith.tarheel
8 months ago by
United States
harold.smith.tarheel4.3k wrote:

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

ADD COMMENTlink modified 8 months ago by Pierre Lindenbaum119k • written 8 months ago by harold.smith.tarheel4.3k
2
gravatar for Pierre Lindenbaum
8 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum119k wrote:

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
ADD COMMENTlink modified 8 months ago • written 8 months ago by Pierre Lindenbaum119k
1

As usual, @PierreLindenbaum has the right tool for the job. Thanks!

ADD REPLYlink written 8 months ago by harold.smith.tarheel4.3k

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

ADD REPLYlink written 8 months ago by RamRS21k

Will do. I cross-posted here only b/c the response time there can be... slow.

ADD REPLYlink written 8 months ago by harold.smith.tarheel4.3k
1

Sorry, what? I was saying that this post would be linked in Pierre's post page, which he has done now.

ADD REPLYlink written 8 months ago by RamRS21k
1

@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.

ADD REPLYlink written 8 months ago by harold.smith.tarheel4.3k

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.

ADD REPLYlink written 8 months ago by RamRS21k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1108 users visited in the last hour