A Step-by-Step Guide to Implementing Machine Learning-Based Variant Filtering with VQSR
By Dr. Zubair Khalid, DVM, MS, PhD ·

Key Takeaways
- Variant Quality Score Recalibration (VQSR) employs a machine learning approach to distinguish true genetic variants from sequencing and alignment artifacts by modeling the joint distribution of variant annotations using known variant resources as training data.
- Successful VQSR implementation necessitates precise alignment of input call sets and training resources to a specific reference genome version (e.g., GRCh38) and requires normalization of variant representation (left-aligned, trimmed) across all files to ensure model accuracy.
- The selection of appropriate training resources (e.g., HapMap, 1000 Genomes, dbSNP) with assigned confidence weights is critical, with a preference for resources that match the study population and reference version to optimize model performance.
- Tranche thresholds, defined by sensitivity to a high-confidence truth set (e.g., HapMap for SNVs), are crucial for balancing the trade-off between false positives and false negatives, with decisions guided by downstream analytical goals (e.g., 99.9% sensitivity for discovery vs. 99.0% for stringent filtering).
- Robust quality control involves documenting the exact VQSR configuration (tool versions, resource paths, annotation sets, tranche thresholds) and tracking key filtering metrics such as the transition transversion ratio and truth set sensitivity to ensure reproducibility and identify potential model failures.
- Common VQSR failure patterns include insufficient training variants (especially for INDELs), reference version mismatches, annotation name discrepancies, and artifacts in low-complexity regions, often requiring troubleshooting by verifying input integrity and configuration before parameter adjustment.
Variant Quality Score Recalibration (VQSR) is a machine learning approach that uses known variant resources to separate true biological variants from sequencing and alignment artifacts. This guide walks through the complete implementation process, from preparing training resources to interpreting tranche output, with concrete decisions for each stage of the workflow. The intended reader is a researcher or laboratory professional who has generated variant calls and needs to filter them with defensible, reproducible criteria.
Understanding the Role of VQSR in Variant Calling Workflows
Variant calling pipelines produce raw variant calls that contain both genuine genetic variation and technical artifacts. The artifact sources include sequencing errors, mapping errors, PCR duplicates that escape removal, and alignment uncertainties in repetitive or homologous regions. Raw call sets typically contain tens of thousands to millions of candidate variants depending on the sequencing depth, library complexity, and genome size. Filtering these calls is a mandatory step before downstream analysis because unfiltered call sets inflate false positive rates and distort allele frequency estimates.
VQSR differs from hard filtering approaches. Hard filtering applies fixed thresholds to annotations such as quality by depth (QD), mapping quality (MQ), and strand bias (FS). These thresholds are static and do not adapt to the specific characteristics of a dataset. VQSR instead builds a Gaussian mixture model from the joint distribution of multiple variant annotations, using a set of known true variants as training labels. The model learns which annotation combinations characterize true variants in your specific dataset and then assigns each variant a probability of being true.
The machine learning approach has practical advantages for whole exome sequencing (WES) and whole genome sequencing (WGS) projects. WES studies are a major approach for uncovering gene-disease associations, and the filtering step directly affects which variants enter association testing [<a href="#ref-1">1</a>]. A poorly filtered call set can produce spurious associations, while an overly aggressive filter removes genuine rare variants that may be disease-relevant. VQSR provides a data-driven middle ground when sufficient training resources exist.
The 1000 Genomes Project phase three data set provides a useful reference for understanding what a well-filtered call set looks like. That project produced biallelic SNVs and INDELs from 2,548 samples spanning 26 populations, called de novo on the GRCh38 assembly [<a href="#ref-2">2</a>]. The project invested substantial effort in filtering and benchmarking, and the resulting call set serves as a training resource for other projects. Your own filtering should aim for a similar standard of transparency and reproducibility.
At a Glance: VQSR Implementation Decision Table
| Decision Point | Standard Practice | Alternative Approach | Key Consideration |
|---|---|---|---|
| Input call set | Joint multi-sample VCF with unfiltered calls | Per-sample VCF for small cohorts | VQSR requires sufficient samples to learn annotation distributions |
| Training resources | HapMap, Omni, 1000 Genomes, dbSNP | Population-specific panels such as 3.5KJPNv2 | Match resources to study population and reference version |
| Annotation set | QD, MQ, FS, SOR, MQRankSum, ReadPosRankSum | Reduced set for non-GATK callers | Choose before running and do not change after inspecting results |
| Tranche threshold | 99.9% sensitivity for discovery | 99.0% for stringent filtering | Balance false positives against false negatives for downstream analysis |
| Variant types | Separate SNV and INDEL recalibration | Combined approach for small datasets | INDEL resources are more limited and require separate handling |
| Reference version | GRCh38 with matching resources | GRCh37 with lifted resources | Mismatched versions cause model failure or degraded results |
Prerequisites and Input Data Requirements
Raw Variant Call Format (VCF) Requirements
VQSR operates on a VCF file that contains variant calls with their associated annotations. The input VCF must include the annotation fields that the model will use as features. Common annotations include QD, MQ, FS, SOR, MQRankSum, ReadPosRankSum, and InbreedingCoeff. The specific annotations available depend on the variant caller used. GATK HaplotypeCaller produces these annotations when the appropriate options are enabled during calling.
The input VCF should contain only unfiltered variant calls. If you have already applied hard filters, those filtered variants are excluded from the VCF and cannot be recovered by VQSR. The standard workflow is to produce a raw VCF, apply VQSR, and then apply any additional hard filters only after VQSR has assigned its scores. This ordering preserves the full annotation distribution for the model to learn from.
Training Resource Requirements
VQSR requires at least one training resource, and the standard practice uses multiple resources with different levels of confidence. The training resources are VCF files containing variants that are known to be true with high confidence. The model uses these variants as positive labels during training.
The most commonly used training resources include:
- HapMap: A set of high-confidence variants genotyped in multiple populations. These variants have been validated through extensive experimental and computational work.
- Omni: A set of variants from the Omni array genotyping platform. These provide a broad set of common variants.
- 1000 Genomes: The phase three call set provides a large set of variants with high confidence [<a href="#ref-2">2</a>].
- dbSNP: A database of known variants maintained by NCBI [<a href="#ref-3">3</a>]. dbSNP contains both validated and unvalidated variants, so it is typically used with a lower confidence weight than HapMap or Omni.
The choice of training resources depends on your study population and variant type. For example, the 3.5KJPNv2 panel from the Tohoku Medical Megabank Project provides allele frequencies for 3,552 Japanese individuals including X chromosome and mitochondrial variants [<a href="#ref-4">4</a>]. If your study focuses on a Japanese population, this panel could serve as a population-specific training resource. The panel was constructed using standard pipelines including 1KGP and gnomAD algorithms to reduce technical biases and allow comparisons to other populations [<a href="#ref-4">4</a>].
Reference Genome and Known Sites Files
VQSR requires the reference genome FASTA file that was used during variant calling. The model uses the reference to compute certain annotations and to ensure consistency between the training resources and your call set. The reference genome version must match between the calling step and the VQSR step. Using a mismatched reference produces errors or silently degraded results.
Known sites files are VCF files that mark sites known to be polymorphic in the population. These are distinct from training resources. Known sites are used to stratify the model and to avoid penalizing true variants that happen to be at known polymorphic sites. The standard practice is to provide a known sites file such as dbSNP for the model to use during the recalibration process.
Preparing Training Resources
Obtaining and Normalizing Resource Files
Training resources must be in the same reference genome version as your call set. If your calls are on GRCh38, the training resources must also be on GRCh38. The 1000 Genomes Project phase three data set was called de novo on GRCh38, which makes it directly usable for GRCh38-based pipelines [<a href="#ref-2">2</a>]. For other reference versions, you may need to liftover the resources, which introduces some positional uncertainty.
Each training resource VCF must be normalized to match your variant representation. This means that variants should be left-aligned and trimmed to the minimal representation. Inconsistent representation between your call set and the training resources causes the model to miss true variants that are represented differently. The GATK tool LeftAlignAndTrimVariants performs this normalization, and the same normalization should be applied to your call set and all training resources.
Assigning Confidence Weights
Each training resource receives a confidence weight that tells the model how much to trust that resource. The weights are specified as the -resource argument with a known=false or training=true flag and a prior probability. The prior probability reflects the expected false positive rate of the resource.
HapMap variants receive the highest confidence because they have been extensively validated. Omni variants receive a slightly lower confidence. 1000 Genomes variants receive a lower confidence still because the call set contains some errors. dbSNP receives the lowest confidence because it contains many unvalidated submissions.
The exact weights depend on your variant type. SNV and INDEL recalibration are performed separately with separate training resources and parameters. INDEL training resources are more limited because fewer high-confidence INDEL resources exist. The 1000 Genomes phase three call set includes INDELs and provides a useful training resource for INDEL recalibration [<a href="#ref-2">2</a>].
Population-Specific Considerations
Training resources built from one population may not fully represent variation in another population. The 3.5KJPNv2 panel was constructed specifically to provide a Japanese population allele frequency resource [<a href="#ref-4">4</a>]. If your study population is underrepresented in the standard training resources, you may need to supplement with population-specific resources.
The tradeoff is between resource completeness and resource confidence. A population-specific panel may contain variants that are absent from HapMap or Omni, but the panel may have a higher error rate than the extensively validated resources. The model can accommodate multiple resources with different weights, so you can include both the standard resources and a population-specific panel with an appropriate confidence weight.
Configuring the VQSR Model
Selecting Annotations for the Model
The choice of annotations directly affects model performance. The standard annotation set for SNV recalibration includes QD, MQ, FS, SOR, MQRankSum, and ReadPosRankSum. These annotations capture different artifact signatures:
- QD: Quality by depth, which penalizes low-quality variants at high depth. Artifacts often have low QD because they are supported by few reads.
- MQ: Mapping quality, which reflects how uniquely reads map to the reference. Artifacts in repetitive regions have low MQ.
- FS: Fisher strand bias, which detects whether variant-supporting reads are biased toward one strand. Strand bias is a common artifact signature.
- SOR: Strand odds ratio, another strand bias measure that is more robust at low depth.
- MQRankSum: The rank sum test for mapping quality, which compares mapping quality of reference and variant-supporting reads.
- ReadPosRankSum: The rank sum test for read position, which detects whether variant-supporting reads are preferentially located at read ends.
The annotation set should be chosen before running the model and should not be changed after inspecting results. Changing annotations after seeing the output constitutes circular analysis and invalidates the statistical interpretation of the tranches.
Setting the Training Set and Truth Set Parameters
The VQSR model distinguishes between training sets and truth sets. Training sets provide positive labels for the model. Truth sets are used to evaluate the model's sensitivity at different tranche thresholds. The truth set should be a very high confidence set of variants that you expect your call set to contain.
The standard practice uses HapMap as a truth set for SNVs because it has the highest confidence. The truth set is used to calculate the transition transversion ratio and to determine the sensitivity at each tranche. The model reports the number of truth set variants retained at each tranche, which lets you assess whether the model is performing reasonably.
Understanding Tranches
Tranches are the output of VQSR. The model assigns each variant a VQSLOD score, which is the log odds of being a true variant versus an artifact. The tranches divide the call set into tiers based on this score. The first tranche contains the highest confidence variants, and subsequent tranches contain progressively lower confidence variants.
The tranche thresholds are specified as sensitivity levels to the truth set. For example, a tranche at 99.9% sensitivity means that the model retains 99.9% of the truth set variants. The corresponding VQSLOD threshold is applied to the full call set. Variants below the threshold are filtered out.
The choice of tranche threshold depends on your downstream analysis. A study looking for rare disease variants may accept a lower sensitivity tranche to reduce false positives. A study doing population frequency estimation may prefer a higher sensitivity tranche to avoid missing rare variants. The tranche file produced by VQSR contains the cumulative statistics for each tranche, which you can use to make this decision.
Running the VQSR Workflow
Step-by-Step Implementation
The VQSR workflow proceeds through several distinct stages. Each stage produces an output file that feeds into the next stage.
Step 1: Prepare the input VCF
Ensure your raw VCF contains all required annotations. If you used GATK HaplotypeCaller with default settings, the standard annotations are present. If you used a different caller, verify that the annotations are present and correctly named. The VCF should be indexed and compressed with bgzip.
Step 2: Normalize variants
Run LeftAlignAndTrimVariants on your call set and all training resources. This ensures consistent variant representation. The normalization step is often overlooked but is critical for model performance.
Step 3: Build the SNP recalibration model
Run VariantRecalibrator with the SNP mode. Specify the training resources with their confidence weights, the annotation list, and the truth set. The output is a recalibration file that contains the model parameters.
Step 4: Apply the SNP recalibration
Run ApplyVQSR with the SNP mode. Specify the tranche threshold you have chosen. The output is a VCF with VQSLOD scores and tranche assignments for each variant.
Step 5: Build and apply the INDEL recalibration model
Repeat steps 3 and 4 for INDELs. INDEL recalibration uses a different annotation set and different training resources. The INDEL-filtered VCF is the final output.
Step 6: Verify the results
Inspect the tranche file and the filtered VCF. Check the transition transversion ratio, the number of variants retained, and the truth set sensitivity. These metrics indicate whether the model performed reasonably.
Using Workflow Management Systems
The VQSR workflow involves multiple steps with specific input-output relationships. Workflow management systems help ensure reproducibility and reduce manual errors. The nf-core project provides community-developed pipelines with standardized configuration and usage documentation [<a href="#ref-5">5</a>]. These pipelines implement best practices for variant calling and filtering, including VQSR where appropriate.
Using a workflow system has several advantages. The pipeline steps are versioned and documented, which supports reproducibility. The configuration files make parameter choices explicit and auditable. The workflow system handles file dependencies and parallelization automatically. For a laboratory that runs variant calling regularly, a workflow system reduces the risk of inconsistent filtering between batches.
The Galaxy Training Network provides accessible workflow training that covers variant calling and filtering concepts [<a href="#ref-6">6</a>]. The training materials walk through the analysis steps in a graphical interface, which is useful for learning the workflow before implementing it on a high-performance computing cluster. The European Bioinformatics Institute also provides training pathways for bioinformatics data resources and practical analysis education [<a href="#ref-7">7</a>].
Containerization and Environment Management
VQSR tools have specific version requirements. The model parameters and tranche calculations can change between tool versions, so the tool version should be recorded and fixed for a project. Containerization with Docker or Singularity ensures that the tool version and dependencies remain consistent across runs.
The Bioconductor project provides reproducible genomic analysis workflows and package documentation [<a href="#ref-8">8</a>]. While Bioconductor is primarily an R package repository, its documentation emphasizes reproducible analysis practices that apply to the VQSR workflow. Recording the exact tool versions, reference genome version, and resource file versions is essential for reproducing the filtering results.
Interpreting VQSR Output
Reading the Tranche File
The tranche file is a text file that reports the cumulative statistics for each tranche. The file contains the tranche threshold, the number of variants retained, the number of truth set variants retained, and the sensitivity to the truth set. The tranche file also reports the transition transversion ratio for each tranche.
The transition transversion ratio is a useful quality metric. Transitions (A-G and C-T) are more common than transversions (A-C, A-T, C-G, G-T) in the human genome. A typical ratio for whole genome data is around 2.0 to 2.1. A much lower ratio suggests that many false positive variants are present. The ratio should increase as you move to higher confidence tranches.
The tranche file also reports the number of novel variants at each tranche. Novel variants are those not present in dbSNP. A high proportion of novel variants in a low confidence tranche suggests that the tranche contains many artifacts. The proportion of novel variants should decrease as you move to higher confidence tranches.
Evaluating Model Performance
The model performance is evaluated by the truth set sensitivity. The truth set contains variants that are known to be true with very high confidence. If the model fails to retain these variants at the expected sensitivity, the model may be misconfigured.
Common signs of model failure include:
- The truth set sensitivity is much lower than expected at the chosen tranche.
- The transition transversion ratio is below 1.5 for the whole call set.
- The number of variants retained is dramatically different from expectations based on the sequencing depth and sample count.
- The model reports an error during training, such as insufficient training variants.
When the model fails, the first step is to check the input files. Verify that the training resources are in the correct reference version and that the variant representation is normalized. Verify that the annotation names match between the call set and the model configuration. Verify that the training resources contain enough variants for the model to learn from.
Using the VQSLOD Score
The VQSLOD score is the primary output of VQSR. Each variant receives a VQSLOD score that represents the log odds of being a true variant. The score is used to assign variants to tranches. Variants in higher tranches have higher VQSLOD scores.
The VQSLOD score can be used for additional filtering beyond the tranche assignment. For example, you might choose to retain only variants in the top two tranches for a particular analysis. The VQSLOD score provides a continuous measure that allows this flexibility.
The VQSLOD score is dataset-specific. A VQSLOD score from one dataset cannot be directly compared to a VQSLOD score from another dataset because the model parameters differ. The tranche assignment is the appropriate unit of comparison between datasets.
Quality Control Metrics and Records
Documenting the VQSR Configuration
The VQSR configuration should be recorded for each run. The record should include:
- Tool versions for VariantRecalibrator and ApplyVQSR
- Reference genome version and file path
- Training resource file paths and versions
- Confidence weights for each training resource
- Annotation list
- Truth set file path
- Tranche threshold
- Date and operator
This record supports reproducibility and troubleshooting. If a downstream analysis produces unexpected results, the VQSR configuration can be reviewed to determine whether the filtering contributed to the problem.
Tracking Filtering Metrics
The filtering metrics should be tracked for each run. The key metrics include:
- Number of variants before filtering
- Number of variants retained at each tranche
- Transition transversion ratio before and after filtering
- Truth set sensitivity at the chosen tranche
- Number of novel variants retained
These metrics provide a baseline for comparing runs. A sudden change in the transition transversion ratio or the number of retained variants may indicate a problem with the input data or the configuration.
Comparing to Published Standards
The 1000 Genomes Project phase three call set provides a benchmark for comparison [<a href="#ref-2">2</a>]. The project reported detailed filtering metrics for their call set, including sensitivity and specificity estimates. Your call set should achieve similar metrics when using comparable data and methods.
The MAGICpipeline protocol provides an example of a complete WES analysis workflow that includes variant filtering [<a href="#ref-1">1</a>]. The protocol describes steps for gene-based rare-variant association analyses and incorporates multiple variant pathogenic annotations [<a href="#ref-1">1</a>]. Reviewing such protocols can help you identify quality control steps that you may have overlooked.
Common Failure Patterns and Troubleshooting
Insufficient Training Variants
VQSR requires a minimum number of training variants to build a stable model. When the training resources contain too few variants, the model may fail to converge or may produce nonsensical tranche assignments. This situation commonly arises for INDEL recalibration because high-confidence INDEL resources are limited.
The solution is to ensure that the training resources contain enough variants. The 1000 Genomes phase three call set includes INDELs and provides a useful resource [<a href="#ref-2">2</a>]. If the standard resources are insufficient, you may need to combine multiple resources or use a population-specific panel such as 3.5KJPNv2 [<a href="#ref-4">4</a>].
Reference Version Mismatch
A common error is using training resources that are in a different reference version than the call set. The model will fail or produce degraded results because the variant positions do not align. The 1000 Genomes phase three data set was called de novo on GRCh38, which makes it directly usable for GRCh38 pipelines [<a href="#ref-2">2</a>]. For other reference versions, you must liftover the resources or obtain version-matched resources.
The reference version should be recorded in the pipeline configuration and verified before running VQSR. A simple check is to compare the contig names and lengths between the call set and the training resources. Mismatched contig names indicate a reference version problem.
Annotation Name Mismatches
VQSR expects specific annotation names in the VCF INFO field. If the variant caller used different annotation names, the model will fail to find the annotations and will error out. This situation commonly arises when using non-GATK variant callers.
The solution is to either configure the variant caller to produce the expected annotation names or to use a tool that converts the annotations. The annotation names should be verified before running VQSR. The VCF header contains the annotation definitions, which can be inspected to confirm the names.
Strand Bias Artifacts in Low Complexity Regions
VQSR reduces but does not eliminate artifacts in low complexity regions. These regions include homopolymers, dinucleotide repeats, and other repetitive sequences. The model may not have enough training variants in these regions to learn their artifact signatures.
The practical consequence is that some artifacts in low complexity regions will pass the VQSR filter. These artifacts can be identified by their annotation values, such as low QD or high FS. If your downstream analysis is sensitive to these artifacts, you may need to apply additional hard filters after VQSR.
Batch Effects in Multi-Sample Calling
VQSR is designed to be applied to a joint call set from multiple samples. When samples are called in separate batches and then combined, the annotation distributions may differ between batches. The model may not perform equally well across batches.
The solution is to call all samples jointly or to apply VQSR separately to each batch and then combine the filtered results. The choice depends on your downstream analysis. Joint calling is preferred for variant discovery because it provides consistent genotype calls across samples. Separate calling is acceptable when the batches are analyzed independently.
Limitations and Interpretation Boundaries
VQSR Requires Sufficient Sample Size
VQSR is designed for multi-sample call sets. The model learns from the joint distribution of annotations across all samples. When applied to a single sample, the model has limited data to learn from and may perform poorly. The standard recommendation is to apply VQSR to call sets with at least 30 samples.
For single-sample or small cohort studies, hard filtering may be more appropriate. The hard filter thresholds can be chosen based on the expected annotation distributions for the sequencing platform and depth. The tradeoff is that hard filters are less adaptive to the specific dataset.
VQSR Does Not Replace Visual Inspection
VQSR assigns probabilities, but it does not confirm that a variant is real. The model can be wrong, particularly for variants that are rare in the training resources. Visual inspection of candidate variants in a genome browser remains an important validation step for variants that will be reported or followed up experimentally.
The inspection should focus on the read alignment at the variant site. Look for consistent base quality, mapping quality, and read position. A variant supported by reads with high mapping quality and consistent base quality is more likely to be real. A variant supported by reads with low mapping quality or biased read position is more likely to be an artifact.
Population-Specific Limitations
The standard training resources are built primarily from populations of European and African ancestry. Variants that are common in other populations may be underrepresented in the training resources. The model may assign lower confidence to these variants because it has not seen similar variants during training.
Population-specific panels such as 3.5KJPNv2 address this limitation for specific populations [<a href="#ref-4">4</a>]. The panel provides allele frequencies for the Japanese population including X chromosome and mitochondrial variants [<a href="#ref-4">4</a>]. If your study population has a similar panel available, incorporating it as a training resource can improve model performance for population-specific variants.
Somatic Variant Calling Considerations
VQSR is designed for germline variant calling. Somatic variant calling has different artifact signatures because the variant allele fraction is often low and the background is heterogeneous. The standard VQSR training resources are not appropriate for somatic calling because they represent germline variation.
Somatic variant calling requires specialized filtering approaches. The ToTem tool provides automated pipeline optimization for somatic variant calling from ultra-deep targeted gene sequencing data [<a href="#ref-9">9</a>]. The tool generates, executes, and benchmarks different variant calling pipeline settings, allowing the optimal pipeline to be selected based on the user's priorities [<a href="#ref-9">9</a>]. For somatic studies, consider using such optimization tools instead of applying germline VQSR directly.
Professional Escalation Criteria
When to Seek Expert Assistance
VQSR implementation can fail in ways that are not obvious from the error messages. Seek expert assistance when:
- The model fails to converge after multiple configuration attempts.
- The transition transversion ratio is below 1.5 after filtering.
- The truth set sensitivity is much lower than expected.
- The filtered call set contains an unexpected number of variants.
- You are unsure whether the training resources are appropriate for your study population.
Expert assistance may come from a bioinformatics core facility, a collaborator with variant calling experience, or a community support forum. The nf-core community provides documentation and support for community-developed pipelines [<a href="#ref-5">5</a>]. The Galaxy Training Network provides training materials that can help you understand the workflow [<a href="#ref-6">6</a>].
When to Abandon VQSR for Hard Filtering
VQSR is not always the best choice. Consider hard filtering when:
- Your cohort has fewer than 30 samples.
- Your sequencing platform produces annotation distributions that differ substantially from the training resources.
- You are working with a non-human organism that lacks training resources.
- You are doing somatic variant calling.
- The VQSR model fails repeatedly despite troubleshooting.
Hard filtering is a valid approach when VQSR is not applicable. The hard filter thresholds should be chosen based on the expected annotation distributions for your sequencing platform and depth. The thresholds should be documented and justified in the methods section of any resulting publication.
When to Consult a Statistical Geneticist
The choice of tranche threshold affects downstream analysis. If you are unsure which tranche threshold to use for your analysis, consult a statistical geneticist. The choice depends on the balance between false positives and false negatives that is acceptable for your specific research question.
A statistical geneticist can also help you understand the interaction between VQSR filtering and downstream statistical methods. For example, rare-variant association tests are sensitive to the variant call set. The MAGICpipeline protocol describes gene-based rare-variant association analyses that incorporate multiple variant pathogenic annotations [<a href="#ref-1">1</a>]. The filtering decisions directly affect which variants enter these analyses.
Building a VQSR Decision Log and Run Comparison System
A recurring problem in VQSR implementation is that researchers often treat each run as an isolated event. When a run produces unexpected tranche assignments or a downstream analysis reveals filtering problems, the operator has no structured way to determine whether the configuration, the input data, or the model itself caused the issue. A decision log that records the rationale behind each configuration choice, combined with a run comparison system that tracks metrics across attempts, turns troubleshooting from guesswork into a systematic process.
What to Record Before the First Run
The decision log should capture the reasoning behind each configuration choice before you execute VariantRecalibrator. This record is distinct from the configuration file itself, which only stores the final parameter values. The log answers the question of why a particular value was chosen, which becomes essential when you revisit the filtering strategy months later or when a collaborator asks about the filtering approach.
For each training resource, record the source, the version or accession date, the reference genome version, and the confidence weight you assigned. Note whether you verified that the resource was normalized to match your call set representation. For the annotation set, record which annotations you selected and why. If you omitted an annotation such as InbreedingCoeff because your caller did not produce it, note that decision and its potential impact on model performance.
The tranche threshold decision deserves particular attention in the log. Record whether you chose 99.9% sensitivity for discovery purposes or 99.0% for stringent filtering, and document the downstream analysis that motivated this choice. A rare-variant association study has different tolerance for false positives than a population frequency estimation project, and the log should reflect that reasoning.
The NCBI provides official documentation for its databases and search systems, which can help you verify that the training resource files you downloaded are the correct versions [<a href="#ref-3">3</a>]. The European Bioinformatics Institute offers training pathways for bioinformatics data resources that include guidance on documenting analysis decisions [<a href="#ref-7">7</a>]. These resources support the record-keeping process by helping you identify what information about a resource is worth capturing.
The Run Comparison Table
A run comparison table tracks the key metrics across every VQSR attempt for a given dataset. This table turns troubleshooting from a memory exercise into a data-driven process. The table should include the run identifier, the date, the configuration version, and the following metrics:
- Number of variants before filtering
- Number of variants retained at the chosen tranche
- Transition transversion ratio before and after filtering
- Truth set sensitivity at the chosen tranche
- Number of novel variants retained
- Model convergence status
- Any error messages or warnings
The transition transversion ratio is a particularly useful diagnostic metric. For human whole genome data, a ratio around 2.0 to 2.1 is typical. If a run produces a ratio below 1.5, the model may have failed to distinguish true variants from artifacts. The run comparison table makes these anomalies visible across attempts instead of requiring you to remember the metrics from each run.
The 1000 Genomes Project phase three call set provides a benchmark for what well-filtered metrics look like. That project produced biallelic SNVs and INDELs from 2,548 samples spanning 26 populations, called de novo on GRCh38 [<a href="#ref-2">2</a>]. Comparing your metrics to this published standard helps you determine whether your filtering performance is within a reasonable range.
A Structured Troubleshooting Sequence
When a run produces unexpected results, work through the following sequence instead of changing parameters at random. This sequence isolates the cause of the problem and prevents you from introducing new issues while troubleshooting.
Step 1: Verify input integrity
Check that the input VCF and all training resources are in the correct reference version. Compare contig names and lengths between the call set and the training resources. Verify that all files are indexed and uncompressed properly. This step catches the most common causes of model failure.
Step 2: Confirm annotation availability
Inspect the VCF header to confirm that all annotations specified in the model configuration are present in the INFO field. If you used a non-GATK variant caller, the annotation names may differ from what VariantRecalibrator expects. This check takes minutes and prevents a common source of errors.
Step 3: Review training resource statistics
Count the number of variants in each training resource and verify that the resources contain enough variants for the model to learn from. INDEL recalibration frequently fails because high-confidence INDEL resources are limited. The 1000 Genomes phase three call set includes INDELs and provides a useful resource for this purpose [<a href="#ref-2">2</a>].
Step 4: Examine the tranche file
The tranche file contains the cumulative statistics for each tranche, including the number of variants retained and the truth set sensitivity. If the truth set sensitivity at your chosen tranche is much lower than expected, the model may be misconfigured. If the transition transversion ratio is below 1.5, the model may have failed to separate true variants from artifacts.
Step 5: Compare to previous runs
Use the run comparison table to determine whether the current run differs from previous attempts. If a configuration change produced the unexpected result, revert that change and test again. If the problem appeared without a configuration change, the input data may have changed.
Recording Post-Filtering Validation
The decision log should extend beyond the VQSR run itself to include post-filtering validation results. Record the number of variants that passed the filter, the proportion of novel variants retained, and any downstream quality checks you performed. This information helps you assess whether the filtering strategy achieved its intended goal.
For example, if you are using the filtered call set for gene-based rare-variant association analyses, record the number of rare variants that entered the association testing. The MAGICpipeline protocol describes procedures for gene-based rare-variant association analyses that incorporate multiple variant pathogenic annotations [<a href="#ref-1">1</a>]. The filtering decisions directly affect which variants enter these analyses, and the log should capture this connection.
The Galaxy Training Network provides accessible workflow training that emphasizes reproducibility and documentation practices [<a href="#ref-6">6</a>]. The nf-core documentation describes community standards for pipeline configuration and usage that support systematic record-keeping [<a href="#ref-5">5</a>]. These resources reinforce the importance of treating the decision log as a core component of the VQSR workflow instead of an optional extra.
Common Failure Patterns in the Decision Log
Reviewing decision logs across multiple projects reveals recurring failure patterns. The most common pattern is changing the annotation set after inspecting the tranche output. This practice constitutes circular analysis because the model is tuned to produce a desired result instead of learned from the data. The decision log prevents this by recording the annotation set before the run and flagging any changes as deviations from the original plan.
Another common pattern is using training resources that do not match the study population. The standard resources are built primarily from populations of European and African ancestry. Population-specific panels such as 3.5KJPNv2 address this limitation for specific populations by providing allele frequencies for the Japanese population including X chromosome and mitochondrial variants [<a href="#ref-4">4</a>]. The decision log should record whether you considered population-specific resources and why you chose to include or exclude them.
A third pattern is applying VQSR to datasets that are too small for the model to learn from. VQSR is designed for multi-sample call sets, and the model learns from the joint distribution of annotations across all samples. For small cohorts, hard filtering may be more appropriate. The decision log should record the sample count and the reasoning behind the choice to use VQSR despite the sample size limitation.
Using the Log for Professional Escalation
When you need to escalate a VQSR problem to a bioinformatics core facility or a collaborator, the decision log and run comparison table provide the context they need to help efficiently. Instead of describing the problem from memory, you can share the configuration rationale, the metrics from each run, and the troubleshooting steps you have already tried.
The nf-core community provides documentation and support for community-developed pipelines [<a href="#ref-5">5</a>]. The Bioconductor project offers reproducible genomic-analysis documentation that emphasizes the importance of recording analysis decisions [<a href="#ref-8">8</a>]. The Carpentries lessons provide foundational training in computing and data practices that support systematic record-keeping [<a href="#ref-10">10</a>]. These resources help you prepare the documentation that makes professional escalation productive.
The decision log also supports the methods section of any resulting publication. Reviewers increasingly expect detailed filtering documentation, and the log provides the information needed to write a transparent methods description. The log should be archived with the project data so that the filtering decisions remain accessible long after the analysis is complete.
Frequently Asked Questions
What is the minimum sample size for VQSR to work reliably?
VQSR is designed for multi-sample call sets. The model learns from the joint distribution of annotations across all samples, and a minimum of approximately 30 samples is commonly recommended. With fewer samples, the model has limited data to learn from and may produce unstable tranche assignments. For small cohorts, hard filtering is a more appropriate approach.
Can VQSR be applied to a single sample?
VQSR can technically be run on a single sample, but the model will have very limited data to learn from. The annotation distribution from a single sample may not represent the full range of true variant and artifact signatures. The tranche assignments from a single-sample VQSR run should be interpreted with caution. For single-sample studies, hard filtering is generally preferred.
What is the difference between a training resource and a truth set?
A training resource provides positive labels for the model during training. The model learns which annotation combinations characterize true variants from the training resource. A truth set is used to evaluate the model's sensitivity at different tranche thresholds. The truth set should be a very high confidence set of variants that your call set is expected to contain. HapMap is commonly used as both a training resource and a truth set.
How do I choose the tranche threshold for my analysis?
The tranche threshold determines the balance between false positives and false negatives. A lower sensitivity tranche (for example, 99.0%) retains fewer variants but has a lower false positive rate. A higher sensitivity tranche (for example, 99.9%) retains more variants but includes more false positives. The choice depends on your downstream analysis. Rare-variant studies may prefer a higher sensitivity tranche to avoid missing true variants. Studies focused on common variants may prefer a lower sensitivity tranche to reduce false positives.
Why does VQSR fail for INDELs?
INDEL recalibration is more challenging than SNV recalibration because high-confidence INDEL training resources are limited. The model may not have enough training variants to learn the INDEL artifact signatures. The 1000 Genomes phase three call set includes INDELs and provides a useful training resource [<a href="#ref-2">2</a>]. If the model still fails, you may need to combine multiple resources or use a population-specific panel.
Can I use VQSR for somatic variant calling?
VQSR is designed for germline variant calling. Somatic variant calling has different artifact signatures because the variant allele fraction is often low and the background is heterogeneous. The standard VQSR training resources are not appropriate for somatic calling. Specialized filtering approaches are needed for somatic data. The ToTem tool provides automated pipeline optimization for somatic variant calling [<a href="#ref-9">9</a>].
What annotations should I include in the VQSR model?
The standard annotation set for SNV recalibration includes QD, MQ, FS, SOR, MQRankSum, and ReadPosRankSum. These annotations capture different artifact signatures including low quality at depth, poor mapping quality, strand bias, and read position bias. The annotation set should be chosen before running the model and should not be changed after inspecting results.
How do I know if my VQSR run worked correctly?
Check the tranche file for the truth set sensitivity, the transition transversion ratio, and the number of variants retained at each tranche. The truth set sensitivity should match the expected value at your chosen tranche. The transition transversion ratio should be around 2.0 to 2.1 for human whole genome data. A much lower ratio suggests that many false positive variants are present.
Related Bioinformatics Guides
- Machine Learning Bioinformatics Projects: From Idea to Publication
- Benchmarking Machine Learning Models in Bioinformatics: Best Practices and Pitfalls
- Machine Learning for Variant Effect Prediction on Protein Stability
- Reinforcement Learning in Bioinformatics: Emerging Applications and Challenges
- Genomic Data Analysis Tools: A Comparative Guide for Researchers
Related Clinical & Scientific Guides
- A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
- Computational Immunology: Modeling the Immune System
- How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices
References and Further Reading
[1] [Protocol for detecting rare and common genetic associations in whole-exome sequencing studies using MAGICpipeline.](https://doi.org/10.1016/j.xpro.2023.102806). 2024. [2] [Variant calling on the GRCh38 assembly with the data from phase three of the 1000 Genomes Project.](https://doi.org/10.12688/wellcomeopenres.15126.2). 2019. [3] [NCBI Data Resources](https://www.ncbi.nlm.nih.gov/). National Center for Biotechnology Information. [4] [3.5KJPNv2: an allele frequency panel of 3552 Japanese individuals including the X chromosome.](https://doi.org/10.1038/s41439-019-0059-5). 2019. [5] [nf-core Documentation](https://nf-co.re/docs). nf-core. [6] [Galaxy Training Network](https://training.galaxyproject.org/). Galaxy Project. [7] [EMBL-EBI Training](https://www.ebi.ac.uk/training). European Bioinformatics Institute. [8] [Bioconductor](https://bioconductor.org/). Bioconductor Project. [9] [ToTem: a tool for variant calling pipeline optimization.](https://doi.org/10.1186/s12859-018-2227-x). 2018. [10] [The Carpentries Lessons](https://carpentries.org/lessons). The Carpentries.This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.