How to Filter a VCF by QUAL with bcftools

By the end of this guide you will have a filtered and indexed VCF of your own variants. You will also have a soft-filtered version that records why each site failed, plus a small set of summary numbers that tell you whether the filter helped or hurt. You need bcftools 1.19 or newer (bcftools --version), a bgzip-compressed and indexed VCF from a caller such as DeepVariant, GATK, or bcftools call, and the same reference FASTA the variants were called against. A 30x human whole genome produces roughly four to five million variants in a 1 to 2 GB compressed VCF, so every command below should finish in minutes on a laptop. If you are working with your own clinical or personal data, treat the output as measurement rather than interpretation, and take anything medically consequential to a clinical geneticist.

1. Find out what QUAL means in your file

QUAL is the sixth column of a VCF. It is a Phred-scaled estimate of the probability that the alternate allele is not really present at that site, so QUAL 30 corresponds to a one in a thousand chance that the variant call is wrong, and QUAL 20 to one in a hundred ¹. The definition is shared across callers. The numeric scale is not.

The consequence is that a threshold copied from a tutorial may mean something different in your file. In bcftools call output, QUAL grows with read depth, so a deep site with a weak variant can outscore a shallow site with a clean one. In GATK output, QUAL is the variant confidence summed over the callset. That is why GATK’s own hard-filtering recommendations pair QUAL < 30 with QD < 2.0, where QD is QUAL divided by the depth of informative reads. In DeepVariant output, QUAL is derived from the genotype likelihood of the called genotype and is already well calibrated, and the caller writes its own FILTER values of PASS or RefCall.

Start by reading the header rather than guessing.

bcftools view -h calls.vcf.gz | grep -E '^##(source|INFO|FILTER|FORMAT)' | head -40

Note which INFO tags exist. If you see QD, FS, MQ, MQRankSum, and ReadPosRankSum, you have GATK annotations and should use them. If the only quality signals are QUAL and FORMAT/GQ, you have a DeepVariant-style file and your options are narrower, which is fine because such callers are designed to do the filtering internally.

2. Look at the QUAL distribution before choosing a cutoff

A threshold is only defensible against a distribution. Extract QUAL and depth with bcftools query, which uses the %TAG syntax.

bcftools query -f '%CHROM\t%POS\t%QUAL\t%TYPE\t%INFO/DP\n' calls.vcf.gz \
  | awk '$3 != "." {print}' \
  > qual_dp.tsv

# quick quantiles of QUAL
cut -f3 qual_dp.tsv | sort -g | awk '{a[NR]=$1} END {
  printf "n=%d p1=%s p5=%s p50=%s p95=%s max=%s\n",
  NR, a[int(NR*0.01)], a[int(NR*0.05)], a[int(NR*0.5)], a[int(NR*0.95)], a[NR]}'

For a 30x human genome from a modern caller, the QUAL histogram is strongly bimodal. There is a dense mode well above 30 containing most true variants and a low tail of noise. Look at where those modes separate in your own file, and split the distribution by variant type, because indels are consistently harder and their QUAL distribution sits lower than that of SNVs ¹.

Also count how many records carry no QUAL at all.

bcftools view -H -i 'QUAL="."' calls.vcf.gz | wc -l

That number matters for the next step.

3. Write the expression, and know how missing values behave

Filtering happens through two options shared by bcftools view, bcftools filter, and most other subcommands: -i includes records matching an expression and -e excludes them. They are not mirror images, because an expression that references a missing value evaluates to false.

# keeps only records with a QUAL value strictly above 30
bcftools view -i 'QUAL>30' calls.vcf.gz -Oz -o inc.vcf.gz

# drops records with QUAL at or below 30, but KEEPS QUAL="."
bcftools view -e 'QUAL<=30' calls.vcf.gz -Oz -o exc.vcf.gz

If your file has records with QUAL=".", these two commands produce different outputs. We prefer -i for quality thresholds, because it fails closed: anything the caller could not score is removed rather than passed through silently. If you want to keep unscored records deliberately, say so in the expression with -i 'QUAL>30 || QUAL="."'.

Use bare tag names in filter expressions. %QUAL was accepted by older versions and is no longer valid in bcftools filter or bcftools view from 1.21 onward, which produces an error like Error occurred while processing --include even though %QUAL remains correct in bcftools query -f. The % prefix belongs to formatting. It does not belong to filtering.

The other syntax trap is the difference between single and double logical operators. && and || are evaluated at the site level and may be satisfied by different samples. & and | require the same sample to satisfy both sides.

# any sample with DP>10, any sample (possibly another) with GQ>20
bcftools view -i 'FMT/DP>10 && FMT/GQ>20' cohort.vcf.gz

# at least one sample with BOTH DP>10 and GQ>20
bcftools view -i 'FMT/DP>10 & FMT/GQ>20' cohort.vcf.gz

For a single-sample genome the two behave identically. Get it right anyway, because the habit transfers to any multi-sample file you later touch.

4. Add the filters that outperform QUAL alone

QUAL is a summary. The annotations behind it separate true variants from artifacts far better, which is why hard-filter recipes are written as conjunctions of several INFO tags rather than one cutoff ¹. For a GATK callset where you do not have enough samples for variant quality score recalibration, GATK’s published hard filters are a reasonable starting point. Applied as a soft filter with bcftools filter, SNVs look like this.

bcftools filter -m + -s SNP_HARD \
  -e 'TYPE="snp" && (QD<2.0 || QUAL<30.0 || SOR>3.0 || FS>60.0 || MQ<40.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0)' \
  calls.vcf.gz -Oz -o soft1.vcf.gz

Indels need their own thresholds, because strand bias and read position statistics behave differently around them.

bcftools filter -m + -s INDEL_HARD \
  -e 'TYPE="indel" && (QD<2.0 || QUAL<30.0 || FS>200.0 || ReadPosRankSum<-20.0)' \
  soft1.vcf.gz -Oz -o soft2.vcf.gz

Here -s NAME writes NAME into the FILTER column instead of deleting the record, and -m + appends to any FILTER value already there rather than replacing it. Chaining two bcftools filter calls this way gives you a per-site audit trail, which is worth far more than a smaller file.

Depth is the other axis worth filtering, in both directions. Very low depth gives unreliable genotypes and very high depth marks collapsed repeats and other regions where the reference is wrong, so both tails should go ². Compute your mean depth first rather than importing a number.

bcftools query -f '%INFO/DP\n' calls.vcf.gz \
  | awk '$1!="." {s+=$1; n++} END {print "mean DP:", s/n}'

bcftools filter -m + -s DEPTH \
  -e 'INFO/DP<10 || INFO/DP>90' soft2.vcf.gz -Oz -o soft3.vcf.gz

A common convention is to cap at roughly two to three times the mean, which for a 30x genome lands near 90. Choose the cap from your own histogram.

For a DeepVariant or similar deep-learning callset, none of the GATK annotations exist and inventing substitutes makes results worse. Keep FILTER="PASS" and add a modest genotype quality floor. Benchmarking of whole-genome pipelines shows that caller choice and the annotations available dominate downstream concordance, so the filter should be matched to the caller rather than to a habit ³.

5. Mask bad genotypes instead of deleting sites

Dropping an entire record because one sample is poorly genotyped throws away information. bcftools filter -S . sets failing genotypes to missing and leaves the site in place.

bcftools filter -S . -e 'FMT/DP<10 | FMT/GQ<20' soft3.vcf.gz -Oz -o masked.vcf.gz

Note the single |, which keeps the expression per-sample. A GQ floor near 20 is the usual compromise, and published quality-control protocols for population-scale sequencing use genotype-level depth and GQ thresholds of this kind before any site-level statistic is computed ⁴. Those protocols also apply per-variant missingness and Hardy-Weinberg tests, both of which need many samples and have no meaning in a single genome. For one person the meaningful filters are QUAL or its underlying annotations, depth, GQ, allele balance, and region masks.

Allele balance is easy to add if your caller reports allelic depths in FORMAT/AD. A true heterozygous site should sit near 0.5.

bcftools filter -m + -s AB \
  -i 'GT="het"' -e 'FMT/AD[0:1] / (FMT/AD[0:0] + FMT/AD[0:1]) < 0.2' \
  masked.vcf.gz -Oz -o ab.vcf.gz

The [0:1] index means first sample, second allele. Multiallelic records break this arithmetic, so split them first with bcftools norm -m -any -f reference.fa.

6. Apply proximity and region filters

Clusters of variants and variants adjacent to indels are a classic alignment artifact rather than biology. bcftools filter handles both directly.

bcftools filter -g 5 -G 10 ab.vcf.gz -Oz -o gaps.vcf.gz

-g 5 soft-filters SNVs within five bases of an indel and -G 10 filters clusters of indels within ten bases of each other. Both write SnpGap and IndelGap into FILTER.

Region masks are the highest-value filter per line of code. Large parts of the genome are unmappable with short reads, and variants called there are mostly wrong regardless of QUAL.

bcftools view -T ^encode_blacklist.bed -Oz -o clean.vcf.gz gaps.vcf.gz
bcftools index --tbi clean.vcf.gz

-T streams the file and respects the ^ negation, while -R uses the index and is faster for a small number of target regions. For variants called from RNA sequencing the equivalent step is stricter still, since coverage follows expression and splice junctions generate systematic false positives near exon boundaries. Studies calling variants from RNA-seq apply hard filters plus removal of clustered and junction-proximal calls before drawing any conclusion ⁵, and the profile of detectable variants depends heavily on library protocol and depth ⁶. Do not reuse WGS thresholds on RNA-seq calls.

7. Check what the filter did

A filter you have not measured is a guess. Three checks are enough for a single genome.

# how many sites carry each FILTER value
bcftools query -f '%FILTER\n' clean.vcf.gz | sort | uniq -c | sort -rn

# transition/transversion ratio and counts by type
bcftools stats clean.vcf.gz > stats.txt
grep -E '^(SN|TSTV)' stats.txt

For human whole-genome data the transition to transversion ratio sits around 2.0 to 2.1. A value near 1.0 means your remaining set is dominated by random sequencing error, and a ratio that climbs steeply as you raise the QUAL threshold tells you the threshold is removing noise. Watching a quality metric like this as a function of the cutoff is the standard way to pick one, and the same logic underlies SNP quality control for association studies ⁷.

The stronger check is a truth set. If you have sequenced a Genome in a Bottle sample, or you have a second callset from a different pipeline on the same DNA, compare them with hap.py or rtg vcfeval and read precision and recall before and after filtering. Independent assessments of variant identification tools consistently find that filtering choices move precision and recall in opposite directions, so the right cutoff depends on which error you can tolerate ⁸. Replicate sequencing of the same individual, such as monozygotic twin comparisons, shows how much apparent variant discordance is pipeline artifact rather than biology ³.

8. Put it together

The full chain for a GATK-style single-sample genome, written as one pipeline with a soft-filtered intermediate kept on disk.

REF=GRCh38.fa

bcftools norm -m -any -f $REF -Ou calls.vcf.gz \
| bcftools filter -m + -s SNP_HARD \
    -e 'TYPE="snp" && (QD<2.0 || QUAL<30.0 || SOR>3.0 || FS>60.0 || MQ<40.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0)' -Ou \
| bcftools filter -m + -s INDEL_HARD \
    -e 'TYPE="indel" && (QD<2.0 || QUAL<30.0 || FS>200.0 || ReadPosRankSum<-20.0)' -Ou \
| bcftools filter -m + -s DEPTH -e 'INFO/DP<10 || INFO/DP>90' -Ou \
| bcftools filter -g 5 -G 10 -Ou \
| bcftools filter -S . -e 'FMT/DP<10 | FMT/GQ<20' \
    -Oz -o soft.vcf.gz --write-index

bcftools view -f PASS -T ^encode_blacklist.bed \
  soft.vcf.gz -Oz -o pass.vcf.gz --write-index

Use -Ou for uncompressed BCF between piped stages, since compressing and decompressing at every step is the single biggest waste of time in bcftools pipelines. Use -Oz only at the end. --write-index builds the CSI or TBI index in the same pass and saves a separate bcftools index call. Keep soft.vcf.gz permanently, because it is the file that lets you revisit a threshold later without recalling variants.

Common problems

bcftools view -i 'QUAL>30' returns far fewer records than expected. Check for missing QUAL with bcftools view -H -i 'QUAL="."' | wc -l. GVCF-derived and some imputed files leave QUAL unset, and -i drops every one of them.

An expression using %QUAL errors out. Drop the percent sign. Formatting strings in bcftools query -f use %QUAL, while filter expressions use QUAL.

bcftools view -f PASS returns nothing. Some callers write . rather than PASS in the FILTER column for records they did not evaluate. Confirm with bcftools query -f '%FILTER\n' | sort | uniq -c, then use -i 'FILTER="PASS" || FILTER="."' if you intend to keep the unannotated ones.

A FMT/AD expression gives wrong ratios or errors. The record is multiallelic and the allele indices no longer mean what you assumed. Run bcftools norm -m -any -f reference.fa first, then filter, and note that normalization requires the exact reference build used for calling.

The filter passes almost everything. You likely applied GATK thresholds to a callset that lacks GATK annotations, where missing tags make each clause false and the -e expression never triggers. Confirm each tag exists in the header before relying on it.

Multi-sample results look strange after genotype masking. Check whether you used | or ||. The double form evaluates at the site level and will filter records rather than individual genotypes, which is rarely what -S . is meant to do.

Footnotes

1. Peter N. Robinson, Rosario M. Piro, Marten Jäger. Variant Calling and Quality-Based Filtering. Computational Exome and Genome Analysis, 2017. https://doi.org/10.1201/9781315154770-12 ↩ ↩² ↩³

2. Joel Martin, Wendy Schackwitz, Anna Lipzen. Genomic Sequence Variation Analysis by Resequencing. Methods in Molecular Biology, 2018. https://doi.org/10.1007/978-1-4939-7804-5_18 ↩

3. Chengyue Tan, Weifen Sun, Jiaxin Yu, et al. Comprehensive assessment of whole-genome sequencing-based variant calling pipelines in monozygotic twins. Briefings in Bioinformatics, 2025. https://doi.org/10.1093/bib/bbaf652 ↩ ↩²

4. Julia M. Sealock, Franjo Ivankovic, Calwing Liao, et al. Tutorial: guidelines for quality filtering of whole-exome and whole-genome sequencing data for population-scale association analyses. Nature Protocols, 2025. https://doi.org/10.1038/s41596-025-01169-1 ↩

5. Anastasia Chumakova, Ivan Vlasov, Elena Filatova, et al. Application of RNA-seq for single nucleotide variation identification in a cohort of patients with hypertrophic cardiomyopathy. Scientific Reports, 2025. https://doi.org/10.1038/s41598-025-03226-x ↩

6. Xia Zhang, Jiawei Liu, Yabing Zhu, et al. A Global Assessment of the Transcription-Dependent Single Nucleotide Variants Relies on the Characteristics of RNA-Sequencing Technologies. Biomolecules, 2026. https://doi.org/10.3390/biom16020211 ↩

7. Philipp E. Bayer, Mitchell Gill, Monica F. Danilevicz, et al. Producing High-Quality Single Nucleotide Polymorphism Data for Genome-Wide Association Studies. Methods in Molecular Biology, 2022. https://doi.org/10.1007/978-1-0716-2237-7_9 ↩

8. Jeffry M. Gaston, Eric J. Alm, An-Ni Zhang. Fast and accurate variant identification tool for sequencing-based studies. BMC Biology, 2024. https://doi.org/10.1186/s12915-024-01891-4 ↩