Analyzing multi-sample VCF with vcfeval

957 views
Skip to first unread message

Joseph S Reddy

unread,
Jun 29, 2016, 5:09:18 PM6/29/16
to RTG Users
I am a first time user of RTG tools. Using version RTG Tools 3.6.2.

I am trying to evaluate annotations that go into determining VQSLOD for VQSR analysis. I have a multi-sample VCF for ~10,000 WES samples. I am trying to implement a workflow discussed by GATK to plot ROC curves for annotations (QUAL, QD, FS, etc.) in their most recent workshop. They've used a single sample to demonstrate the efficacy of VQSR filtering and how it compares to hard-filtered  or unfiltered callsets. While I was trying to implement the workflow on my multi-sample VCF, I got an error message " Error: No sample name provided but calls is a multi-sample VCF." 

Can I use vcfeval to analyze more than one sample at a time using --sample=<calls_SampleName>? (I have an underscore "_" in my sample names, is this going to be an issue, if RTG uses it as a delimiter?)
If analysis of more than one sample is permitted, do I have to explicitly provide a complete list of samples that I wish to analyze using --sample=<calls_sampleName> 

My interest lies in specifically looking at each of the annotations (QD, QUAL, FS and others) in the "INFO" column of my VCF file and plot ROC curves for SNPs and INDELs. 

Could you please provide recommendations on how best to get vcfeval to work on my multi-sample VCF? 

Thanks,
Joseph. 

Len Trigg

unread,
Jun 29, 2016, 5:57:35 PM6/29/16
to Joseph S Reddy, RTG Users
Hi Joseph,

Vcfeval effectively compares the specific genomic haplotypes that the baseline and called VCFs provide, for a specific sample, and so it is inherently a pairwise comparison. You can easily specify which sample to select from a multisample VCF via --sample as you suggest, although it is normally assumed that the same sample name is used in both the baseline and the calls file. If they differ for some reason, you can use the form --sample=sample1,sample2, where sample1 is the desired sample from the baseline VCF, and sample2 is the desired sample from the calls VCF (so ',' is the delimiter, and the "_" in your sample names won't cause any problem).

However, for ROC analysis, the samples being compared must refer to the same underlying genome, otherwise you can't make appropriate conclusions about correctness of the calls.  If you were to compare (say) GIAB NA12878 as the baseline, and some non-NA12878 sample from the call set, a FP could actually be a true variant in the new sample, similarly a FN could just be that the new genome doesn't contain the variant that NA12878 has. The more disparity there is between the underlying genomes being compared, the more unreliable the ROC (what you tend to see is that the curves get flattened toward the same diagonal you would get with an unpredictive scoring attribute).

Cheers,
Len.






--
You received this message because you are subscribed to the Google Groups "RTG Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rtg-users+...@realtimegenomics.com.
Visit this group at https://groups.google.com/a/realtimegenomics.com/group/rtg-users/.

Joseph S Reddy

unread,
Jun 30, 2016, 7:10:24 PM6/30/16
to RTG Users
Thank you Len for patiently explaining everything. 

After inputs from you and Geraldine at Broard (GATK team), I was able to understand how to effectively use vcfeval to analyze annotations and plot ROC curves. Since we've decided to use gold standard NA12878 sample as our baseline, I'm going to perform joint genotyping on our samples while including NA12878 in the mix. I will then use vcfeval to perform a pairwise comparison of haplotypes in the baseline to the multi-sample VCF for NA12878. 

Thank you once again,

Cheers,
Joseph.  

KIRTI

unread,
Mar 3, 2018, 2:00:41 PM3/3/18
to RTG Users
Can u provide me the link of gold standard data of NA12878 ?

KIRTI

unread,
Mar 3, 2018, 2:05:28 PM3/3/18
to RTG Users


On Thursday, June 30, 2016 at 2:39:18 AM UTC+5:30, Joseph S Reddy wrote:

Len Trigg

unread,
Mar 4, 2018, 3:10:37 PM3/4/18
to KIRTI, RTG Users
On Sun, 4 Mar 2018 at 08:00 KIRTI <kirti...@gmail.com> wrote:
Can u provide me the link of gold standard data of NA12878 ?

For the GIAB truth sets, you probably want:
https://github.com/genome-in-a-bottle/giab_latest_release

And for Illumina Platinum Genomes:
https://github.com/Illumina/PlatinumGenomes

Cheers,
Len.


--
Reply all
Reply to author
Forward
0 new messages