Vcf Variant Filtering

Common Pitfalls

1. Computing statistics on raw VCF files without any filtering. This is the most common and most consequential mistake. Raw Ts/Tv ratios can be 1.5 or lower, suggesting much worse data quality than actually exists, and variant counts will be inflated by thousands of false positives.

   • How to avoid: Always check the FILTER column before computing statistics. If FILTER is . for all records, apply QUAL>=30 filtering first.

2. Assuming FILTER='.' means the variant passed. The dot in the FILTER column means "not assessed," not "passed." Treating unassessed variants as high-quality leads to inclusion of many false positives.

   • How to avoid: Distinguish between . (unassessed), PASS (passed all filters), and named filters (failed). Use bcftools query -f '%FILTER\n' | sort | uniq -c to inspect.

3. Writing custom Python parsers for Ts/Tv calculation. Custom parsers frequently mishandle multi-allelic sites, complex variants, or edge cases in the VCF format. They also tend to be slower than optimized C-based tools.

   • How to avoid: Use bcftools stats for Ts/Tv computation. It is well-tested, fast, and handles all VCF edge cases correctly.

4. Double-filtering a VCF that was already filtered. Applying QUAL>=30 to a VCF that already went through VQSR or hard filtering can remove true variants that were assigned lower QUAL scores by the caller but passed model-based assessment.

   • How to avoid: Inspect the VCF header for filter command records (##bcftools_viewCommand, ##GATKCommandLine) and the FILTER column for non-dot values before adding more filters.

5. Using the wrong Ts/Tv expectation for the assay type. Comparing a WES Ts/Tv of 3.0 against a WGS expectation of 2.1 would incorrectly suggest the data is too good, while comparing a WGS Ts/Tv of 2.0 against WES expectations would incorrectly suggest poor quality.

   • How to avoid: Know the expected Ts/Tv range for your assay type. WGS: 2.0-2.1; WES: 2.8-3.3; targeted panels: variable depending on target regions.

6. Forgetting to report the filtering status in results. Presenting a Ts/Tv ratio without stating whether or how the data was filtered makes the number uninterpretable and non-reproducible.

   • How to avoid: Always include a statement like "Ts/Tv computed after QUAL>=30 filtering (N variants passed)" or "Ts/Tv computed on unfiltered raw calls (N total variants)."

7. Filtering on QUAL alone when richer annotations are available. GATK and other callers provide INFO-level annotations (QD, FS, MQ, SOR, MQRankSum, ReadPosRankSum) that are more informative than QUAL alone. Relying solely on QUAL when these are available leaves quality on the table.

   • How to avoid: Check whether INFO annotations are present. If so, consider GATK hard filtering recommendations or VQSR as a more powerful alternative to simple QUAL filtering.

Workflow

1. Inspect the VCF metadata and FILTER column

   • Check the VCF header for caller information and existing filter records
   • Examine the FILTER column distribution

   # View header for caller and filter information
   bcftools view -h input.vcf | grep -E '##(FILTER|source|GATKCommandLine|bcftools)'
   
   # Check FILTER column values
   bcftools query -f '%FILTER\n' input.vcf | sort | uniq -c | sort -rn | head

2. Assess QUAL score distribution

   • Determine the proportion of low-QUAL variants

   # Check QUAL distribution
   bcftools query -f '%QUAL\n' input.vcf | awk '{if($1<30) low++; else high++} END {print "QUAL<30:", low+0, "QUAL>=30:", high+0}'

   • Decision point: If most variants are QUAL>=30 and FILTER=PASS, skip to Step 4. If FILTER=. and many low-QUAL variants exist, proceed to Step 3.

3. Apply quality filtering

   • Use the appropriate threshold (default QUAL>=30 unless otherwise specified)

   # Apply QUAL>=30 filter
   bcftools view -i 'QUAL>=30' input.vcf -Oz -o filtered.vcf.gz
   bcftools index filtered.vcf.gz

4. Compute summary statistics

   • Extract Ts/Tv ratio, variant counts, and other metrics

   # Ts/Tv ratio from filtered VCF
   bcftools view -i 'QUAL>=30' input.vcf | bcftools stats - | grep ^TSTV | cut -f5
   
   # Full summary numbers
   bcftools view -i 'QUAL>=30' input.vcf | bcftools stats - | grep ^SN
   
   # Per-sample statistics
   bcftools stats -s - input.vcf | grep ^PSC

5. Compare before and after filtering

   • Count variants before and after to quantify the impact of filtering

   # Count variants before filtering
   bcftools view -H input.vcf | wc -l
   
   # Count variants after filtering
   bcftools view -i 'QUAL>=30' input.vcf | bcftools view -H | wc -l
   
   # Compare Ts/Tv before and after
   echo "=== Before filtering ==="
   bcftools stats input.vcf | grep ^TSTV
   echo "=== After QUAL>=30 filtering ==="
   bcftools view -i 'QUAL>=30' input.vcf | bcftools stats - | grep ^TSTV

6. Validate and report

   • Check that Ts/Tv is within the expected range for the assay type (WGS: 2.0-2.1, WES: 2.8-3.3)
   • If Ts/Tv is still below expected range after filtering, investigate sample quality or consider stricter thresholds
   • Report the filtering applied, number of variants before and after, and the resulting statistics
   • Include the bcftools command used so the analysis is fully reproducible

Protocol Guidelines

1. Pre-analysis inspection is non-negotiable. Before any VCF analysis, run bcftools query -f '%FILTER\n' and examine the QUAL distribution. This takes seconds and prevents hours of wasted analysis on unreliable data.

2. Default to QUAL>=30 for unlabeled VCFs. When the filtering status is ambiguous and the user has not specified a threshold, QUAL>=30 is the standard community default. Document this choice explicitly in the output.

3. Never silently filter. If you apply filtering, always report it. If you choose not to filter, state why. Transparency in filtering decisions is fundamental to reproducible genomics.

4. Use piped commands for efficiency. Chain bcftools view and bcftools stats in a pipe rather than writing intermediate files. This is faster, uses less disk, and reduces the chance of analyzing the wrong file.

5. Verify with Ts/Tv as a sanity check. After any filtering operation, compute Ts/Tv as a quality control metric. If the ratio does not improve or remains outside the expected range, reassess the filtering strategy or investigate upstream data quality issues.

Further Reading

• bcftools documentation -- Comprehensive reference for bcftools commands, filtering expressions, and statistics output format
• VCF specification (v4.3+) -- Formal definition of VCF format including FILTER column semantics, QUAL field, and INFO annotations
• GATK Best Practices for Germline Short Variant Discovery -- Recommended filtering strategies including VQSR and hard filtering thresholds for GATK-called variants

Related Skills

• bcftools-variant-manipulation -- Detailed bcftools usage for subsetting, merging, annotating, and manipulating VCF files beyond simple filtering
• gatk-variant-calling -- Upstream variant calling with GATK HaplotypeCaller, including VQSR and hard filtering configuration
• samtools-bam-processing -- BAM file quality control and processing steps that precede variant calling and affect downstream VCF quality