The current GATK version is 3.5-0

#### Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

# Variant Quality Score Recalibration (VQSR)

Posts: 122Dev ✭✭✭
edited March 16

This document describes what Variant Quality Score Recalibration (VQSR) is designed to do, and outlines how it works under the hood. The first section is a high-level overview aimed at non-specialists. Additional technical details are provided below.

For command-line examples and recommendations on what specific resource datasets and arguments to use for VQSR, please see this FAQ article. See the VariantRecalibrator tool doc and the ApplyRecalibration tool doc for a complete description of available command line arguments.

As a complement to this document, we encourage you to watch the workshop videos available in the Presentations section.

## High-level overview

VQSR stands for “variant quality score recalibration”, which is a bad name because it’s not re-calibrating variant quality scores at all; it is calculating a new quality score that is supposedly super well calibrated (unlike the variant QUAL score which is a hot mess) called the VQSLOD (for variant quality score log-odds). I know this probably sounds like gibberish, stay with me. The purpose of this new score is to enable variant filtering in a way that allows analysts to balance sensitivity (trying to discover all the real variants) and specificity (trying to limit the false positives that creep in when filters get too lenient) as finely as possible.

The basic, traditional way of filtering variants is to look at various annotations (context statistics) that describe e.g. what the sequence context is like around the variant site, how many reads covered it, how many reads covered each allele, what proportion of reads were in forward vs reverse orientation; things like that -- then choose threshold values and throw out any variants that have annotation values above or below the set thresholds. The problem with this approach is that it is very limiting because it forces you to look at each annotation dimension individually, and you end up throwing out good variants just because one of their annotations looks bad, or keeping bad variants in order to keep those good variants.

The VQSR method, in a nutshell, uses machine learning algorithms to learn from each dataset what is the annotation profile of good variants vs. bad variants, and does so in a way that integrates information from multiple dimensions (like, 5 to 8, typically). The cool thing is that this allows us to pick out clusters of variants in a way that frees us from the traditional binary choice of “is this variant above or below the threshold for this annotation?”

Let’s do a quick mental visualization exercise (pending an actual figure to illustrate this), in two dimensions because our puny human brains work best at that level. Imagine a topographical map of a mountain range, with North-South and East-West axes standing in for two variant annotation scales. Your job is to define a subset of territory that contains mostly mountain peaks, and as few lowlands as possible. Traditional hard-filtering forces you to set a single longitude cutoff and a single latitude cutoff, resulting in one rectangular quadrant of the map being selected, and all the rest being greyed out. It’s about as subtle as a sledgehammer and forces you to make a lot of compromises. VQSR allows you to select contour lines around the peaks and decide how low or how high you want to go to include or exclude territory within your subset.

How this is achieved is another can of worms. The key point is that we use known, highly validated variant resources (omni, 100 Genomes, hapmap) to select a subset of variants within our callset that we’re really confident are probably true positives (that’s the training set). We look at the annotation profiles of those variants (in our own data!), and we from that we learn some rules about how to recognize good variants. We do something similar for bad variants as well. Then we apply the rules we learned to all of the sites, which (through some magical hand-waving) yields a single score for each variant that describes how likely it is based on all the examined dimensions. In our map analogy this is the equivalent of determining on which contour line the variant sits. Finally, we pick a threshold value indirectly by asking the question “what score do I need to choose so that e.g. 99% of the variants in my callset that are also in hapmap will be selected?”. This is called the target sensitivity. We can twist that dial in either direction depending on what is more important for our project, sensitivity or specificity.

## Technical overview

The purpose of variant recalibration is to assign a well-calibrated probability to each variant call in a call set. This enables you to generate highly accurate call sets by filtering based on this single estimate for the accuracy of each call.

The approach taken by variant quality score recalibration is to develop a continuous, covarying estimate of the relationship between SNP call annotations (QD, SB, HaplotypeScore, HRun, for example) and the the probability that a SNP is a true genetic variant versus a sequencing or data processing artifact. This model is determined adaptively based on "true sites" provided as input (typically HapMap 3 sites and those sites found to be polymorphic on the Omni 2.5M SNP chip array, for humans). This adaptive error model can then be applied to both known and novel variation discovered in the call set of interest to evaluate the probability that each call is real. The score that gets added to the INFO field of each variant is called the VQSLOD. It is the log odds ratio of being a true variant versus being false under the trained Gaussian mixture model.

The variant recalibrator contrastively evaluates variants in a two step process, each performed by a distinct tool:

• VariantRecalibrator
Create a Gaussian mixture model by looking at the annotations values over a high quality subset of the input call set and then evaluate all input variants. This step produces a recalibration file.

• ApplyRecalibration
Apply the model parameters to each variant in input VCF files producing a recalibrated VCF file in which each variant is annotated with its VQSLOD value. In addition, this step will filter the calls based on this new lod score by adding lines to the FILTER column for variants that don't meet the specified lod threshold.

Please see the VQSR tutorial for step-by-step instructions on running these tools.

### How VariantRecalibrator works in a nutshell

The tool takes the overlap of the training/truth resource sets and of your callset. It models the distribution of these variants relative to the annotations you specified, and attempts to group them into clusters. Then it uses the clustering to assign VQSLOD scores to all variants. Variants that are closer to the heart of a cluster will get a higher score than variants that are outliers.

### How ApplyRecalibration works in a nutshell

During the first part of the recalibration process, variants in your callset were given a score called VQSLOD. At the same time, variants in your training sets were also ranked by VQSLOD. When you specify a tranche sensitivity threshold with ApplyRecalibration, expressed as a percentage (e.g. 99.9%), what happens is that the program looks at what is the VQSLOD value above which 99.9% of the variants in the training callset are included. It then takes that value of VQSLOD and uses it as a threshold to filter your variants. Variants that are above the threshold pass the filter, so the FILTER field will contain PASS. Variants that are below the threshold will be filtered out; they will be written to the output file, but in the FILTER field they will have the name of the tranche they belonged to. So VQSRTrancheSNP99.90to100.00 means that the variant was in the range of VQSLODs corresponding to the remaining 0.1% of the training set, which are basically considered false positives.

### Interpretation of the Gaussian mixture model plots

The variant recalibration step fits a Gaussian mixture model to the contextual annotations given to each variant. By fitting this probability model to the training variants (variants considered to be true-positives), a probability can be assigned to the putative novel variants (some of which will be true-positives, some of which will be false-positives). It is useful for users to see how the probability model was fit to their data. Therefore a modeling report is automatically generated each time VariantRecalibrator is run (in the above command line the report will appear as path/to/output.plots.R.pdf). For every pair-wise combination of annotations used in modeling, a 2D projection of the Gaussian mixture model is shown.

The figure shows one page of an example Gaussian mixture model report that is automatically generated by the VQSR from the example HiSeq call set. This page shows the 2D projection of mapping quality rank sum test versus Haplotype score by marginalizing over the other annotation dimensions in the model.

In each page there are four panels which show different ways of looking at the 2D projection of the model. The upper left panel shows the probability density function that was fit to the data. The 2D projection was created by marginalizing over the other annotation dimensions in the model via random sampling. Green areas show locations in the space that are indicative of being high quality while red areas show the lowest probability areas. In general putative SNPs that fall in the red regions will be filtered out of the recalibrated call set.

The remaining three panels give scatter plots in which each SNP is plotted in the two annotation dimensions as points in a point cloud. The scale for each dimension is in normalized units. The data for the three panels is the same but the points are colored in different ways to highlight different aspects of the data. In the upper right panel SNPs are colored black and red to show which SNPs are retained and filtered, respectively, by applying the VQSR procedure. The red SNPs didn't meet the given truth sensitivity threshold and so are filtered out of the call set. The lower left panel colors SNPs green, grey, and purple to give a sense of the distribution of the variants used to train the model. The green SNPs are those which were found in the training sets passed into the VariantRecalibrator step, while the purple SNPs are those which were found to be furthest away from the learned Gaussians and thus given the lowest probability of being true. Finally, the lower right panel colors each SNP by their known/novel status with blue being the known SNPs and red being the novel SNPs. Here the idea is to see if the annotation dimensions provide a clear separation between the known SNPs (most of which are true) and the novel SNPs (most of which are false).

An example of good clustering for SNP calls from the tutorial dataset is shown to the right. The plot shows that the training data forms a distinct cluster at low values for each of the two statistics shown (haplotype score and mapping quality bias). As the SNPs fall off the distribution in either one or both of the dimensions they are assigned a lower probability (that is, move into the red region of the model's PDF) and are filtered out. This makes sense as not only do higher values of HaplotypeScore indicate a lower chance of the data being explained by only two haplotypes but also higher values for mapping quality bias indicate more evidence of bias between the reference bases and the alternative bases. The model has captured our intuition that this area of the distribution is highly enriched for machine artifacts and putative variants here should be filtered out!

### Tranches and the tranche plot

The recalibrated variant quality score provides a continuous estimate of the probability that each variant is true, allowing one to partition the call sets into quality tranches. The main purpose of the tranches is to establish thresholds within your data that correspond to certain levels of sensitivity relative to the truth sets. The idea is that with well calibrated variant quality scores, you can generate call sets in which each variant doesn't have to have a hard answer as to whether it is in or out of the set. If a very high accuracy call set is desired then one can use the highest tranche, but if a larger, more complete call set is a higher priority than one can dip down into lower and lower tranches. These tranches are applied to the output VCF file using the FILTER field. In this way you can choose to use some of the filtered records or only use the PASSing records.

The first tranche (90) which has the lowest value of truth sensitivity but the highest value of novel Ti/Tv, is exceedingly specific but less sensitive. Each subsequent tranche in turn introduces additional true positive calls along with a growing number of false positive calls. Downstream applications can select in a principled way more specific or more sensitive call sets or incorporate directly the recalibrated quality scores to avoid entirely the need to analyze only a fixed subset of calls but rather weight individual variant calls by their probability of being real. An example tranche plot, automatically generated by the VariantRecalibrator walker, is shown below.

This is an example of a tranches plot generated for a HiSeq call set. The x-axis gives the number of novel variants called while the y-axis shows two quality metrics -- novel transition to transversion ratio and the overall truth sensitivity.

Note that the tranches plot is not applicable for indels and will not be generated when the tool is run in INDEL mode.

### Ti/Tv-free recalibration

We use a Ti/Tv-free approach to variant quality score recalibration. This approach requires an additional truth data set, and cuts the VQSLOD at given sensitivities to the truth set. It has several advantages over the Ti/Tv-targeted approach:

• The truth sensitivity (TS) approach gives you back the novel Ti/Tv as a QC metric
• The truth sensitivity (TS) approach is conceptual cleaner than deciding on a novel Ti/Tv target for your dataset
• The TS approach is easier to explain and defend, as saying "I took called variants until I found 99% of my known variable sites" is easier than "I took variants until I dropped my novel Ti/Tv ratio to 2.07"

We have used hapmap 3.3 sites as the truth set (genotypes_r27_nr.b37_fwd.vcf), but other sets of high-quality (~99% truly variable in the population) sets of sites should work just as well. In our experience, with HapMap, 99% is a good threshold, as the remaining 1% of sites often exhibit unusual features like being close to indels or are actually MNPs, and so receive a low VQSLOD score.
Note that the expected Ti/Tv is still an available argument but it is only used for display purposes.

### Finally, a couple of Frequently Asked Questions

#### - Can I use the variant quality score recalibrator with my small sequencing experiment?

This tool is expecting thousands of variant sites in order to achieve decent modeling with the Gaussian mixture model. Whole exome call sets work well, but anything smaller than that scale might run into difficulties.

One piece of advice is to turn down the number of Gaussians used during training. This can be accomplished by adding --maxGaussians 4 to your command line.

maxGaussians is the maximum number of different "clusters" (=Gaussians) of variants the program is "allowed" to try to identify. Lowering this number forces the program to group variants into a smaller number of clusters, which means there will be more variants in each cluster -- hopefully enough to satisfy the statistical requirements. Of course, this decreases the level of discrimination that you can achieve between variant profiles/error modes. It's all about trade-offs; and unfortunately if you don't have a lot of variants you can't afford to be very demanding in terms of resolution.

#### - Why don't all the plots get generated for me?

The most common problem related to this is not having Rscript accessible in your environment path. Rscript is the command line version of R that gets installed right alongside. We also make use of the ggplot2 library so please be sure to install that package as well. See the Common Problems section of the Guide for more details.

Post edited by Geraldine_VdAuwera on
Tagged:

• Posts: 230Member ✭✭✭

I would suggest starting off with removing MQ from the model. that's usually been my achilles heel when dealing with targeted capture. I've had MQRankSum bite me a couple of times (for a different reason).

• StanfordPosts: 5Member

@Kurt Thanks! I will try that.

mikel

• BarcelonaPosts: 40Member ✭✭

I found this sentence VERY confusing, if not wrong (I am still not 100% sure):

The first tranche (from the bottom, with lowest values) is exceedingly specific but less sensitive, and each subsequent tranche in turn introduces additional true positive calls along with a growing number of false positive calls.

Should it not be the first tranche from the TOP? i.e. in tranche-90 everything is considered TP (at least in this figure).

• United KingdomPosts: 400Member ✭✭✭

@SteveL No, starting from "tranche-90" and going up your specificity will decrease and your sensitivity will increase. The sentence is correct.

• BarcelonaPosts: 40Member ✭✭

Thanks @tommycarstensen - however the issue I have is that the image (immediately following paragraph with quoted sentence) on page one of this thread has the "tranche-90" at the top, so it isn't clear, to novices, what the word "bottom" refers to - also the "lowest" values, could be interpreted as tranche or the ti/tv, which unfortunately go in opposing directions. I just think it could be stated more clearly, but maybe it's just me.

#### Issue · Github April 2015 by Geraldine_VdAuwera

Issue Number
20
State
closed
Last Updated
Closed By
vdauwera

@SteveL I see what you mean, it's a fair point. We'll try to clarify.

Geraldine Van der Auwera, PhD

• Posts: 6Member

Is the model better informed by using a VCF containing all the chromosomes or does the model still work efficiently when the recalibration is done per chromosome?

• United KingdomPosts: 400Member ✭✭✭

@shiroann It is better to maximize the number of data points by running simultaneously across all chromosomes. That being said I think VQSR still might give reasonable results for individual chromosomes.

• Posts: 6Member
edited April 2015

Thanks @tommycarstensen . I'm assuming the per chromosome VQSR will be greatly influenced by the number of variants in that particular chromosome.I wonder if that introduces some underlying variation in the VQSLOD score sensitivity across the genome/exome?

Post edited by shiroann on
• United KingdomPosts: 400Member ✭✭✭
edited April 2015

@shiroann Yes, basically don't run VQSR per chromosome, unless you need the results quickly, before the Death Star reaches Alderaan or something like that. And you would have to run ApplyRecalibration separately per chromosome as well, if you were to do it. JUST don't DO IT (not a shoe company TM).

... Unless you are dealing with chromosomes, which you think might have annotations very different from your other chromosomes. I think I remember @Geraldine_VdAuwera telling me to run across all chromosomes irrespective of these differences. Here is the thread I was thinking of:

http://gatkforums.broadinstitute.org/discussion/2895/vqsr-and-sex-chromosomes


I posted a figure in that thread showing the the chromY variants more frequently fail to the DP annotation for obvious reasons:

http://cd8ba0b44a15c10065fd-24461f391e20b7336331d5789078af53.r23.cf1.rackcdn.com/gatk.vanillaforums.com/FileUpload/99/99c914e1b209e5615f9de0c07ca55f.png

Post edited by tommycarstensen on

We really don't recommend running VQSR per chromosome, no. It's going to cause what are essentially batch effects between chromosomes of different sizes, among other issues.

Geraldine Van der Auwera, PhD

• Warwick University, CoventryPosts: 16Member

I posted a question about VQSR and issues with the VariantRecalibrator and VariantAnnotator about two weeks ago; but it still hasn't appeared anywhere so I'm concerned that it has been blocked outright. I'm not re-posting it here right now; because it's possible it was blocked earlier due to my inadvertently making multiple posts,

William

@WVNicholson
Hi William.

I just verified your account. I don't see your question anywhere, so it is best to re-post again. I should be able to see it now.

-Sheila

• Posts: 43Member

Hi,

This might have been asked before (sorry...) but if I have less than 30 exome samples and matching 1000 genomes samples is going to be difficult (because I am running my samples through my pipeline and samples will differ each time according to enrichment method, intervals, ethnicity etc.) is it better to leave out the VQSR step completely or carry on using it on <30 samples?

Thanks

Kath

@Kath I would recommend running some tests and evaluating your results. Generally speaking we think it's worth the additional effort to put together a VQSR-worthy dataset.

Geraldine Van der Auwera, PhD

• Warwick University, CoventryPosts: 16Member

As mentioned earlier, I attempted to post the message below before; so it could actually still be in the system somewhere. (In particular, I sent it to "Ask the GATK Team" which may be distinct from the forums.)

I'm trying to carry out VQSR on a dataset. When I run the VariantRecalibrator it fails to find annotations, even when added either by running the VariantAnnotator or by re-running GenotypeGVCFs with the arguments for the annotations to be generated. (The variants were originally called using the HaplotypeCaller.) My most recent command line for the VariantRecalibrator is:

java -jar /home/u1374090/software/gatk/gatk_only/GenomeAnalysisTK.jar \
-T VariantRecalibrator -R ../ncbi/ncbi_sorghum_combo_genome.fasta \
-input S_all_modern_combo_jointan.vcf \
-resource:ncbi_variations,known=false,training=true,truth=true,prior=10.0 ../ncbi/variations/nstd63_Zheng_et_al_2011.Sorbi1.submitted.variant_call.germline.wvned.vcf \
-an DP -an FS -an MQRankSum \
-mode SNP \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \
-recalFile recalibrate_SNP.recal \
-tranchesFile recalibrate_SNP.tranches \
-nt ${ntPara} The shell variable ntPara was set to "10". I've tried several variations of the above command line, all without success. What happens is that the final .recal and .tranches file are empty and I get the following output including an error message (which I've edited because it is quite long): INFO 17:36:29,136 HelpFormatter - Program Args: -T VariantRecalibrator -R ../ncbi/ncbi_sorghum_combo_genome.fasta -input S_all_modern_combo_jointan.vcf -resource:ncbi_variations,known=false,training=true,truth=true,prior=10.0 ../ncbi/variations/nstd63_Zheng_et_al_2011.Sorbi1.submitted.variant_call.germline.wvned.vcf -an DP -an FS -an MQRankSum -mode SNP -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 -recalFile recalibrate_SNP.recal -tranchesFile recalibrate_SNP.tranches -nt 10 INFO 17:36:29,141 HelpFormatter - Executing as u1374090@hamilton on Linux 3.13.0-46-generic amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_25-b17. . . . INFO 17:36:29,715 GenomeAnalysisEngine - Strictness is SILENT INFO 17:36:29,832 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000 INFO 17:36:29,951 MicroScheduler - Running the GATK in parallel mode with 10 total threads, 1 CPU thread(s) for each of 10 data thread(s), of 64 processors available on this machine INFO 17:36:30,027 GenomeAnalysisEngine - Preparing for traversal INFO 17:36:30,030 GenomeAnalysisEngine - Done preparing for traversal INFO 17:36:30,031 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING] INFO 17:36:30,031 ProgressMeter - | processed | time | per 1M | | total | remaining INFO 17:36:30,031 ProgressMeter - Location | sites | elapsed | sites | completed | runtime | runtime INFO 17:36:30,035 TrainingSet - Found ncbi_variations track: Known = false Training = true Truth = true Prior = Q10.0 INFO 17:36:48,155 VariantDataManager - DP: mean = NaN standard deviation = NaN INFO 17:36:51,231 GATKRunReport - Uploaded run statistics report to AWS S3 ##### ERROR ------------------------------------------------------------------------------------------ ##### ERROR A USER ERROR has occurred (version 3.3-0-g37228af): . . . ##### ERROR MESSAGE: Bad input: Values for DP annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations. See http://gatkforums.broadinstitute.org/discussion/49/using-variant-annotator • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @WVNicholson Hi, It looks like your training file may not have the DP annotation in it. Can you check and make sure it has it? -Sheila • Warwick University, CoventryPosts: 16Member Indeed it looks like my training (VCF) file does not have the DP annotation or apparently any annotations that would be useful for VQSR. I'm still looking; but at this stage it's possibly none of the publicly available sets of known SNPs, at least as available in VCF format actually have the required annotations. It appears I may be able to get data in SRA format (apparently convertible to BAM) for some of the SNPs in the VCF file from Ensembl Plants. So presumably I can make my own annotated VCF from that one and the SRA files (after converting to BAM) by using VariantAnnotator. This looks like it's going to be an involved process; so I thought I should check if that's a sensible course of action? William • Posts: 9,857Administrator, Dev admin @WVNicholson There seems to have been a miscommunication in our team about what is required for VQSR. It is not necessary for the training resource file to contain any annotations; the program only uses the location of the sites from that file. Geraldine Van der Auwera, PhD • Warwick University, CoventryPosts: 16Member Okay, I'm not sure where to go from here then. Inspection of the input VCF file (from my data rather than the other VCF file used for training) by eye appears to suggest it does actually have the DP annotation. The training file does not have the DP annotation; but you just said the training file is not required to contain any annotations. The only thing I can think of is that somehow I managed to call a lot of variants in the new data and none of them with agree with the known variants in the training set. That seems a bit unlikely; but would it give an error message like the one above? William • Warwick University, CoventryPosts: 16Member On further investigation, my training VCF file may not be suitable. It looks like it has no SNPs - just indels and something described as "DUP" in the ALT field (CNVs, I think), William • Warwick University, CoventryPosts: 16Member It turns out I have another VCF file with known SNPs and MNPs that should be suitable. However, when I attempt to use that I get a different error message: ##### ERROR MESSAGE: Input files /home/u1374090/sorghum/gatkwork/../ensembl/variations/sorghum_bicolor.snps_mnps_only_wvned2.vcf and reference have incompatible contigs: Relative ordering of overlapping contigs differs, which is unsafe. I'm confused by the error message; because the contigs in the various files refer to chromosomes and should be treated as non-overlapping. I can however re-order the data for the contigs in the training VCF to match those in my VCF with my variant calls; but the error message may be a clue that there are other problems, William #### Issue · Github May 2015 by Geraldine_VdAuwera Issue Number 987 State closed Last Updated Closed By vdauwera • Posts: 9,857Administrator, Dev admin @WVNicholson Sorry for the confusion. What the error message means by "overlapping contigs" is "contigs that are present in both files". It doesn't mean to imply that the sequences overlap. We'll try to improve the phrasing to clarify this. Geraldine Van der Auwera, PhD • Posts: 261Member ✭✭ INFO=<ID=VQSLOD,Number=1,Type=Float,Description="Log odds ratio of being a true variant versus being false under the trained gaussian mixture model"> This doc does not say how the VQSLOD is computed. Is VQSLOD the log odds ratio or log odds? My understanding is that the term "odds ratio" always associates with at least two factors. • Posts: 9,857Administrator, Dev admin It's the log odds ratio of being a true positive vs. being a false positive. See this presentation (~slide 12): https://drive.google.com/open?id=0BwTg3aXzGxEDanNtNFgyamtaMXM&authuser=0 Geraldine Van der Auwera, PhD • Posts: 9,857Administrator, Dev admin Whoops sorry wrong version, see this one instead, which has the new slides: https://drive.google.com/open?id=0BwTg3aXzGxEDNW5tWUhSM2hZRlk&authuser=0 Geraldine Van der Auwera, PhD • Posts: 261Member ✭✭ edited June 2015 @Geraldine_VdAuwera Thanks so much for the slides! If my understanding is correct, the VQSLOD is the log(p/q) where p is the prob of being TP and q is the prob of being FP (although p and q are from different model, so p!= 1-q). I think it should be called log odds since it is not the ratio of two odds, as compared to the log odds ratio which is (p/(1-p)) / (q/(1-q)) and should be stated as the log odds ratio of being true under the Positive model, instead of the log odds ratio of being true. But just my two cents! Post edited by blueskypy on • Posts: 9,857Administrator, Dev admin Oh, I may have been abusing the terminology on this one. I'll consult with our resident mathketeers for a definitive answer Geraldine Van der Auwera, PhD • Posts: 9,857Administrator, Dev admin Yep, I can now confirm that I was wrong to call it the log odds ratio, it's just log odds, period. For a bit of additional mathy detail: It does not matter that p + q != 1 since you can re-normalize: Let p' = p/(p+q) and q' = q/(p+q). Now p' + q' = 1 and p/q = p'/q' so log(p/q) can rightly be called log odds. Blueskypy is right that this is different from the log odds ratio which would be log((p/(1-p)) / (q/(1-q))) but unless Geraldine used that term in her discussion [G: oops] it should be irrelevant. And they remind me that VQSLOD = "VQS-Log-ODds" Geraldine Van der Auwera, PhD • Posts: 261Member ✭✭ @Geraldine_VdAuwera I'm glad that I helped a bit to the great work you guys has been doing! So probably the header of the VCF file needs to be changed a little: INFO=<ID=VQSLOD,Number=1,Type=Float,Description="Log odds ratio of being a true variant versus being false under the trained gaussian mixture model"> #### Issue · Github June 2015 by Geraldine_VdAuwera Issue Number 1015 State closed Last Updated Assignee Array Closed By vdauwera • Posts: 9,857Administrator, Dev admin @blueskypy I think that's the same thing as what's currently there... do you mean to say we should take "ratio" out? I would agree with that. Geraldine Van der Auwera, PhD • Posts: 261Member ✭✭ @Geraldine_VdAuwera Yes, that's what I mean. • Posts: 9,857Administrator, Dev admin OK we'll fix that in the next release. I put in a ticket so it doesn't fall off the radar. Geraldine Van der Auwera, PhD • Posts: 15Member dumb question, but would it be erroneous to combine my cases and controls together for joint genotyping and even VQSR? • Posts: 9,857Administrator, Dev admin @nchuang Not erroneous at all, that is what you should do. Geraldine Van der Auwera, PhD • BeijingPosts: 6Member edited July 2015 Hi, I have been working with GATK recently and I just came across this situation, VariantRecalibrator produced tranches.pdf as shown in the attachment. I am computing four human exome samples and I used UnifiedGenotyper to call SNPs (four samples together). Here is my code: java -XX:+UseParallelGC -XX:ParallelGCThreads=15 -Xmx15g -jar ../GenomeAnalysisTK.jar -T VariantRecalibrator -R hg19.fa -input$one.snp_indel.raw.vcf -resource:hapmap,VCF,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.b37.sites.vcf.gz -resource:omni,VCF,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.b37.sites.vcf.gz -resource:dnsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_135.b37.vcf.gz -an QD -an HaplotypeScore -an MQRankSum -an ReadPosRankSum -an MQ -recalFile var_plots/$one.recal -tranche 100 -tranche 99.9 -tranche 99.0 -tranche 90 -tranchesFile ./var_plots/$one.tranches -rscriptFile ./var_plots/one.plots.R -mode SNP Thanks for any help! Zeyu Post edited by zzy on • United KingdomPosts: 400Member ✭✭✭ @zzy I think you forgot to ask your question. • BeijingPosts: 6Member edited July 2015 @tommycarstensen said: zzy I think you forgot to ask your question. Thank you and I am sorry. My question is, why my attached tranches.pdf file is full of false positives? Is this due to any mistakes I made during the pipeline or just because my data are less qualified or insufficient sample size? Thanks again, @tommycarstensen Zeyu Post edited by zzy on • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @zzy Hi Zeyu, We recommend using at least 30 exome samples in VQSR. If you don't have 30 samples, you can either add data from 1000 Genomes project or use hard filtering. -Sheila • BeijingPosts: 6Member @Sheila said: zzy Hi Zeyu, We recommend using at least 30 exome samples in VQSR. If you don't have 30 samples, you can either add data from 1000 Genomes project or use hard filtering. -Sheila Hi, @Sheila , Thanks for your reply! Zeyu • ItalyPosts: 18Member Hi there, I have some questions for you. I'm trying to use VQSR on my data and I'm working on Target sequencing extracted from a WES (using a .list file). To be clear, I have 3 datasets, the first one with 3 samples (more or less 20M TOTAL reads per sample extracted from TruSight panel), the second one with 48 samples (1,5M TOTAL reads per sample, target) and the last one with 6 samples (20M reads per sample, extracted from TruSight). Reading on your forum and your previous answers, you suggest to use VQSR on at least 1 sample WGS or 30 WES samples. As concern targeted sequencing, I read this link: https://www.broadinstitute.org/gatk/guide/article?id=39 you suggest to work on --maxGaussians options if I want to use VQSR, so my question is: Can I use VQSR on all my datasets by working on --maxGaussians parameter? If no, can I try to merge all my vcf coming out from my 3 datasets (target and exome) and using them as a sort of db to use VQSR on small target? I hope I was clear and sorry for my dummy question. Thank you, Matteo • Posts: 73Member edited September 2015 What does 'culprit=x' mean for variants that get the PASS vqsr filter? I take it that for variants that fail, it indicates the annotation that caused it to not make the PASS filter? Is the PASS filter assigned based on a particular VQSLOD cutoff? I have done two analyses with using UG and HC. The latter has ~ 400 samples less than the former (the rest are identical). I have some variants that pass the vqsr form the HC results, but didn't in the UG analysis (possibly also the other way around). I am trying to determine if I should trust these variants. For example, one particular variant was given ReadPosRankSum=43.912 in the UG analysis, and failed (culprit=ReadPosRankSum), but ReadPosRankSum=1.16 in the HC analysis, and PASSed. Is ReadPosRankSum calculated by ranking distance to end of read for all reads in all samples? If so, then the difference in score could be caused by the 400 sample difference? In such cases, do you recommend any other sanity checks to have confidence in these variants? Thanks Vicky Post edited by vsvinti on • BarcelonaPosts: 40Member ✭✭ Hi Vicky, Yes, I understand that culprit means the same i.e. that the filter applied following your UG analysis was based on ReadPosRankSum - the value in the example you give shows there was a very strong difference within the reads of the variant calls versus the reference calls in the case of the UG analysis. However it seems strange, assuming that you used the same BAMs as input, that the difference should be so great. Unless the ~400 samples extra that you added were called using a different capture kit? (I assume this was WES). You also don't mention how many samples you have in the UG analysis? Also, were all other steps in your workflow the same, alignment and joint-gentoyping? I would not be surprised if different versions of various tools would lead to different findings due to batch effects of some sort. Steve • Posts: 73Member edited September 2015 @SteveL Thanks for your message! I have ~ 2200 samples in the UG analysis, ~ 1800 in the HC. These are exomes, and yes, a mixture of 2 target captures. The 400 excluded were all with one capture, but I have ~ 600 left with the same capture as the excluded ones (in the HC analysis). The VQSR for UG was done with gatk2.8 - I remember I had done everything with 2.7 but the vqsr performed terribly on my data, and had to update to the 2.8 release when it came out. I have done all the re-analysis including the re-alignments for all the samples, so I shouldn't have any differences in workflow. UG was done by polled calling across all the samples, HC via the new single-sample gvcf, batches of 100 and then calling across these. So I guess I have a different ratio of samples coming from the different capture platforms for UG and HC. Do you think that could cause the issue? Is there a way to deal with this? Vicky Post edited by vsvinti on • BarcelonaPosts: 40Member ✭✭ Hi Vicky, I have no experience with UG myself, only HC, so I guess we'll have to ask those who really no. @Geraldine_VdAuwera @Sheila ? It is not surprising that a few variants fail here and there, given that the batches were different, and there are always some extreme outliers in any NGS analysis comparisons in my experience. Which percentage of your variants that were PASS in the smaller dataset are no longer PASS in the larger dataset? Steve • Posts: 73Member 22% of variants from UG not in HC (~ 170,000) 4.5% of variants from HC not in UG (~ 30,000) (comparing chr pos for PASS variants) Thanks for your input I'll see what the others think .. • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @vsvinti Hi, Are you asking about results from version 2.8? Or, when you say you re-ran the analysis, do you mean you used the latest version? Unfortunately, it is really hard to comment on results from version 2.8. I was not even working here at the time! In general, the differences could certainly be due to the way the annotations were calculated by Haplotype Caller and Unified Genotyper. Also remember, Haplotype Caller does a reassembly step which may shift the reads arounds. http://gatkforums.broadinstitute.org/discussion/4146/hc-step-2-local-re-assembly-and-haplotype-determination -Sheila • Posts: 73Member edited September 2015 Thanks, @Sheila I am asking about whether I should trust the new variants that pass VQSR in the HC analysis. I guess if I hadn't done a previous analysis with UG, I might not have questioned the results, but I am getting some signals that look they might be too good to be true with these new variants. Quoting one snp as eg (mentioned above): For example, one particular variant was given ReadPosRankSum=43.912 in the UG analysis, and failed (culprit=ReadPosRankSum), but ReadPosRankSum=1.16 in the HC analysis, and PASSed. Would the difference in score be due to the 400 sample difference, or local reassembly, or a number of other factors? Would you expect to see such a big discrepancy between the calculated value between two runs? 2.8 was not that long ago, when HC was not do-able for the size of my data... Post edited by vsvinti on • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @vsvinti Yes, the reassembly could explain the large difference in the ReadPosRankSum. You can take a look at the bamout file of the region and compare it to the original bam file. https://www.broadinstitute.org/gatk/guide/article?id=5484 Of course, the best thing to do is run again with Haplotype Caller in GVCF mode. https://www.broadinstitute.org/gatk/guide/article?id=3893 -Sheila • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @Matteodigg Hi Matteo, I am not sure if this is the right assumption, but it sounds like you have 57 whole exome samples. You have subset the exome to the targets of interest, and you want to run VQSR on those targets. If that is the case, it is best to simply run VQSR on the whole exomes themselves, then subset the filtered VCF to your targets of interest. -Sheila • Posts: 73Member edited September 2015 Thanks @Sheila I still need some degree of confidence in my results, and I am becoming a bit unsure what what comes out of VQSR as PASS (it seems to change a lot..). Also, the proportion of FP in my 90s tranch seems to be quite high - what could be causing this? Your example shows only TPs in the 90s, so I'm quite far off that! I ploted this using the target_titv 3.2 (which is recommended for exomes).. Thanks! Vicky ps. yes, I am dealing with exomes and two capture platforms, I pushed everything through the same workflow, and I am only analysing the intersection of captures + 50bp padding Post edited by vsvinti on • Posts: 9,857Administrator, Dev admin Hi @vsvinti, It's not unexpected to see substantial differences in results between UG and HC. The variant modeling is more accurate in HC, and reads are realigned locally, so many annotations values will be different, and this in turn impacts the modeling in VQSR. In addition, having such a significant difference in cohort size, and a mix of capture platforms, makes it difficult to do a direct comparison. You also don't mention whether you have done any evaluation against a common truth set, which would at least give you some objective indication of which callset is closer to the expected truth. Relative numbers of variants are just that, relative. In brief, we can't give you any measure of confidence (good or bad) based on this information, aside from saying that we have greater confidence in our newer tools and in our latest version. If you want greater confidence, you'll need to do a more thorough evaluation. Geraldine Van der Auwera, PhD • Posts: 73Member • ItalyPosts: 18Member Hi @Sheila , thank you for your answer, maybe I have not been so clear. I have 9 whole exome sequencing (ok, I can use all the exome, not only target for variant calling) but also 48 target sequencing and I'd like to use VQSR on all of these data. So, is it a good Idea to use 9 exome + 48 target at the same time in VQSR ? I have also another question for you, if I have 2 vcf files with 10 samples each, can I use CombineGVCF to merge them and then use VQSR? I'm quite sure that I cannot use more than one input file in VQSR, right? Thank you very much for your patience! • Posts: 9,857Administrator, Dev admin @Matteodigg No on both. The error modes are different depending on the experimental design; if you combine them you will confuse the model and obtain inferior results. Geraldine Van der Auwera, PhD • EdinburghPosts: 1Member @SteveL said: I found this sentence VERY confusing, if not wrong (I am still not 100% sure): The first tranche (from the bottom, with lowest values) is exceedingly specific but less sensitive, and each subsequent tranche in turn introduces additional true positive calls along with a growing number of false positive calls. Should it not be the first tranche from the TOP? i.e. in tranche-90 everything is considered TP (at least in this figure). I also found this very confusing, even though I think it is now clear what the first tranche is : "The first tranche (from the bottom, with the highest value of truth sensitivity but usually the lowest values of novel Ti/Tv) is exceedingly specific but less sensitive." Shouldn't the tranche with the highest value of truth sensitivity be exceedingly sensitive and less specific? • GreensboroPosts: 50Member I am a little confused/in dillema with what I should do/do not for VQSR. In my situation I don't have database of highly confident variants (that is used for training the recalibration model) for VQSR process. Based on what we did in the BQSR process I think similar model can applied for VQSR process (but with recalibration step only one time). Is it ok if do the following process: 1. Variant Calling in each biological sample using HC and UG using strict filtering options. 2. Merge the variants (from HC and UG) and call the common (intersecting) variants called by both the callers for this particular biological sample. https://www.broadinstitute.org/gatk/guide/article?id=53 3. Do the above two steps for all other biological samples (I have 6 biological sample from same population). 4. Now, merge the variants from different biological samples and call the common (intersecting) variants. If the biological samples belong to same family and/or population my hunch is that we should now have a very high and confident variants that belong to that family and/or population. So, can I use this highly confident variant (at least as non-true sites) in further VQSR process for training the recalibration model. https://www.broadinstitute.org/gatk/guide/article?id=1259 Thanks in advance ! • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin Sorry for the confusion! You are correct. I just tried to make the sentence a little clearer. Hopefully it makes sense now. Thank you for the catch. -Sheila • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin Are you working with 6 whole exome samples or 6 whole genome samples? We recommend using VQSR with at least 30 whole exome samples or 1 whole genome sample. You might be better off using hard filtering. We have some basic recommendations here: http://gatkforums.broadinstitute.org/discussion/2806/howto-apply-hard-filters-to-a-call-set However, you may need to play around with the thresholds to better suit your data. -Sheila • SeattlePosts: 2Member edited November 2015 I have a VCF file that has several custom annotations(ABHet, blat score, etc) that were calculated with the in-house software and that I've added to VCF INFO field using bcf-annotate. I would like to re-run VQSR recalibration on this VCF file using those custom annotations. So far I had no luck getting it run without the error when using custom annotations (the same command with GATK annotations like DP, QD, FS, etc works just fine). Here is the error message I get when attempting to run VQSR with custom annotations: ##### ERROR MESSAGE: Invalid command line: No tribble type was provided on the command line and the type of the file 'file_1.recal' could not be determined dynamically. Please add an explicit type tag: NAME listing the correct type from among the supported types  So, I'm wondering if it is even possible to build a Gaussian model using custom annotations or is VQSR recalibration allowed only using annotations generated by GATK e.g.: QD\DP\FS\SOR|MQRankSum\etc? Post edited by aachueva on • Broad InstitutePosts: 3,084Member, Broadie, Moderator, Dev admin @aachueva Hi, Can you tell us which version of GATK you are using and the exact commands you tried to run? Thanks, Sheila • SeattlePosts: 2Member edited December 2015 Hi Sheila, thanks for the response. Im using GenomeAnalysisTK-3.4-46. Here are the commands im using (NOTE: -an foo.ABHet and -an foo.uniq20 are my custom annotations that are in the vcf file): PREFIX="file_1" VERSION="v1" java -jarGATK_DIR/GenomeAnalysisTK.jar \
-T VariantRecalibrator \
-R $REF_GENOME \ --input:VCF$PREFIX".vcf" \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 $REF_DIR/hapmap_3.3.b37.vcf.gz \ -resource:omni,known=false,training=true,truth=true,prior=12.0$REF_DIR/1000G_omni2.5.b37.vcf.gz \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 $REF_DIR/1000G_phase1.snps.high_confidence.b37.vcf.gz \ -resource:dbsnp,known=true,training=false,truth=false,prior=2.0$REF_DIR/dbsnp_138.b37.vcf.gz \
-mode SNP \
-an QD \
-an MQRankSum \
-an MQ \
-an SOR \
-an FS \
-an BaseQRankSum \
-an GQ_MEAN \
-an foo.ABHet \
-an foo.uniq20 \
-mG 10 \
-mNG 3 \
-tranche 100.0 \
-tranche 99.9 \
-tranche 99.8 \
-tranche 99.7 \
-tranche 99.6 \
-tranche 99.5 \
-tranche 99.4 \
-tranche 99.3 \
-tranche 99.2 \
-tranche 99.1 \
-tranche 99.0 \
-tranche 98.0 \
-tranche 97.0 \
-tranche 96.0 \
-tranche 95.0 \
-tranche 90.0 \
-recalFile $PREFIX"_"$VERSION".recal" \
-tranchesFile $PREFIX"_"$VERSION".tranches" \
-rscriptFile  $PREFIX"_"$VERSION".R"

java -jar $GATK_DIR/GenomeAnalysisTK.jar \ -T ApplyRecalibration \ -R$REF_GENOME \
--input:VCF $PREFIX".vcf" \ --ts_filter_level 99.5 \ -recalFile$PREFIX"_"$VERSION".recal" \ -tranchesFile$PREFIX"_"$VERSION".tranches" \ -mode SNP \ -o$PREFIX"_"\$VERSION"_output.vcf.gz"


Post edited by aachueva on

@aachueva
Hi,

It should be possible to use custom annotations; you just have to make sure the annotations you specify exactly match the annotation names in the VCF. Can you confirm that the .recal file is not empty?

-Sheila

• GreensboroPosts: 50Member

@Sheila said:
@everestial007
Hi,

Are you working with 6 whole exome samples or 6 whole genome samples? We recommend using VQSR with at least 30 whole exome samples or 1 whole genome sample. You might be better off using hard filtering. We have some basic recommendations here: http://gatkforums.broadinstitute.org/discussion/2806/howto-apply-hard-filters-to-a-call-set However, you may need to play around with the thresholds to better suit your data.

-Sheila

HI Sheila,
Sorry for late reply but just had chance to start working on my data again.
I am working with the 6 whole genome sample/population; 2 populations (so 12 whole genome samples in total).
I have done the filtering part (using stringent filtering parameters) and collected highly confident variants for each population. But I am also interested in rare variants within-between populations to work on EHH (extendend haplotype homozygosity).
I have also called variants in g.vcf mode for each sample.
So, I am wondering if it is possible to take the highly confident variants identified from these samples to run VQSR (I am thinking of it similar to BQSR process) on the g.vcf file, which has the possibility of identifying rare variants. I am not sure if its a good question to ask, but I am not sure if I can take VQSR process in this way. Any suggestions would be helpful?

Thank you,

I think Geraldine's answer here sums up the issue.

Good luck!

-Sheila

• GreensboroPosts: 50Member

Thanks @Sheila
The discussion on the link you posted is obviously going to be helpful.

Thank you !

• GreensboroPosts: 50Member

The drop box link mentioned in the beginning of the article isn't working.

Yes, that does not work anymore. I removed it from the document. The best thing to do is check out the slides from the talks (from the Events Webpage). The latest BroadE workshop slides and presentations will be uploaded soon. But, you can have a look at the BroadE workshop slides and talks from March here.

-Sheila

• GreensboroPosts: 50Member

@Sheila said:
@everestial007
Hi,

Yes, that does not work anymore. I removed it from the document. The best thing to do is check out the slides from the talks (from the Events Webpage). The latest BroadE workshop slides and presentations will be uploaded soon. But, you can have a look at the BroadE workshop slides and talks from March here.

-Sheila

Thank you !

• OstravaPosts: 4Member

Hi,
I work with exome-capture data of several non-model fish species from one genus, and their hybrids. I called variants with HC in GVCF mode and then joint genotyped them using GenotypeGVCFs tool. This first analysis was done on 18 priority samples of 2 species and their hybrids (and one outgroup). Another about 60 samples are currently running through the HC pipeline and I will repeat the whole procedure with them included - I'm just learning this tool on smaller dataset I have at hand.

I tried this VQSR on our 18 samples. For training I used our own database of high-confidence SNPs (~30k) that are diagnostic between the two species (parental species are homozygous either for ref or alt and hybrids are heterozygous). I've set tranches as described in the best practice document, i.e. -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 and also --maxGaussians 4. I also excluded DP and InbreedingCoeff annotations from the call (we don't have just one population so we don't expect HWE).

However the final recalibrated_snps.vcf contains only PASS or VQSRTrancheSNP99.90to100.00 filters, no other tranches seem to be present in my data. How is it possible? Could it be that my training dataset is too small?

I was also wondering if the VQSR is OK with data from more than one species (I've excluded InbreedingCoeff for this reason) and for that matter - what about the rest of GATK pipeline? As far as I understand the tools this shouldn't be a problem for most of them, possibly only VQSR, am I right?

Thanks in advance for any insights.

Jan

• OstravaPosts: 4Member

Hi again,
just realised my previous question was quite dumb because there is only one trenche in my data due to ApplyRecalibration procedure. My bet

The question about multiple species (same genus) in VQSR still stands though

Thanks

Jan

Hi @jena, no worries.

All the GATK tools do require that you run your analysis relative to a single reference genome. If you can meet that condition then it should be fine. But to be honest we don't have any direct experience with such a use case so we won't be able to advise you much.

Geraldine Van der Auwera, PhD

Hi, we have been trying to use VQSR to annotate/classify our callset. I think we understand the nature of VQSLoD score as well as the source of the "true positives" during training, but just wonder what would be the "false positives" that are fed to the GMM?

Thanks,

Chao