How to Filter VCF Files by Quality and Depth

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-10

By the end of this guide you will be able to take a raw Variant Call Format file from GATK, DeepVariant, FreeBayes, or bcftools mpileup and reduce it to a set of variants you can defend. You will filter on site quality, read depth, genotype quality, allele balance, strand bias, and missingness, and you will check that the filtering removed noise rather than signal.

You need a few things on hand before starting. A Linux or macOS machine with bcftools (1.18 or newer), htslib, and tabix installed. Optionally vcftools and GATK4. The VCF itself, ideally bgzipped and indexed. The reference FASTA used for calling, plus its .fai index, because normalization needs it. And roughly twice the VCF size in free disk space, since you will write several intermediate files.

Nothing here assumes a human cohort. The same logic applies to companion animal genomes, livestock, crops, and pathogens, though the threshold values differ and the resource databases thin out fast outside human and mouse.

The bioinformatics file formats you touch before filtering

Filtering a VCF is the last step in a chain of formats. Knowing what each one carries explains why a low base quality in a FASTQ file eventually becomes a low QUAL score in a VCF.

FASTQ vs FASTA: which file carries quality scores

FASTQ is the raw read format. Each record uses four lines: an @ header line, the nucleotide sequence, a + separator line, and a quality string with one ASCII character per base. The quality characters are Phred scores offset by 33. The character I is ASCII 73, so the Phred score is 73 minus 33, which equals 40. A Phred score of 40 means a 1 in 10,000 chance that the base is wrong.

FASTA is the sequence-only format. Each record has a > header line followed by sequence, and there is no quality string at all. The reference genome you align against is FASTA, and so are individual gene or transcript sequences. There is no meaningful "quality" for a reference base. This matters practically because a VCF header records the reference file used, and if you feed a FASTQ file to a tool expecting a reference FASTA, or the reverse, you get a failure that is easy to misread as a filtering bug. The distinction is simple: FASTQ carries per-base confidence, FASTA does not.

BAM, CRAM, VCF, and gVCF

BAM and CRAM are aligned reads. They carry the FASTQ base qualities plus a mapping quality (MAPQ) per read, which is a Phred-scaled confidence that the read is placed correctly.

VCF is a text format for variant sites. It has a header of ## metadata lines, a single #CHROM column header line, and then one line per variant. BCF is the compressed binary equivalent with the same logical content. gVCF adds reference blocks that summarize non-variant positions, and those blocks carry a depth but usually no meaningful QUAL.

A caller takes reads from a BAM or CRAM, compares them to a reference FASTA, and writes a VCF. The quality chain runs from FASTQ base quality, to read MAPQ, to the site-level and genotype-level scores in the VCF. That is why a VCF filter is a proxy for read-level quality, and why filters on mapping quality and depth catch problems that base quality alone cannot.

The VCF columns that matter for filtering

Here is a single biallelic record, split across two lines for readability.

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS    ID     REF  ALT  QUAL  FILTER  INFO                                FORMAT        SAMPLE1        SAMPLE2
chr1    10177  .      A    AC   100   PASS    DP=45;MQ=59;QD=2.2;FS=0.0;SOR=0.7  GT:AD:DP:GQ:PL  0/1:20,25:45:99:120,0,300  0/0:44,0:44:99:0,120,360

Read the record left to right. CHROM and POS give the locus. REF and ALT give the alleles. QUAL is 100, a Phred-scaled score meaning the caller thinks there is about a 1 in 10^10 chance this site is not variant at all. FILTER is PASS. INFO holds site-level summaries. FORMAT names the per-sample fields, and each sample column holds one value per field in that order.

QUAL is a site-level statement. It answers one question: is there a variant here in at least one sample. It does not tell you whether a particular sample's genotype is correct. That job belongs to GQ, the genotype quality.

Site-level, sample-level, and derived fields

The table below is the working vocabulary for the rest of this guide.

Two of these are derived. QD is QUAL divided by INFO/DP, so a site with QUAL 100 and depth 80 has QD 1.25. That site has plenty of reads but the caller is not confident relative to that depth, which usually means the alt allele is spread thinly or inconsistently across reads. A site with QUAL 100 and depth 45 has QD 2.22. Same QUAL, different story, and QD is what exposes it. VariFAST uses exactly this style of reasoning, combining several false-positive indicators into a weighted score rather than relying on a single threshold [1].

AB is computed from AD. With AD of 20,25 the het allele balance is 25/45, or 0.56. With AD of 45,2 the balance is 0.04, which is what a mismapped read or a contaminating sample looks like, not a clean heterozygous site.

Three checks before you filter anything

Run these first. Each one has saved a full rerun.

Check the build and contig names. The VCF must match the reference you will normalize against. Compare contig names and lengths.

bcftools view -h cohort.vcf.gz | grep '^##contig' | head
cut -f1,2 reference.fa.fai | head

If the VCF says chr1 and the FASTA says 1, normalization will fail or silently skip. Fix the naming with bcftools annotate, rename-chrs using a two-column map file before you do anything else.

Check the samples. Confirm the names are what you expect and there are no duplicates.

bcftools query -l cohort.vcf.gz | sort | uniq -c | sort -rn | head

Check whether the caller already filtered. A callset that arrives with a populated FILTER column is not raw.

bcftools query -f '%FILTER\n' cohort.vcf.gz | sort | uniq -c | sort -rn

If you see thousands of records labeled LowQual alongside PASS, the caller applied its own soft filter. Honor it, question it, or discard it deliberately. Do not pretend the file is unfiltered. Tools such as vcfView exist specifically so that variants which failed caller filters can be re-examined before you throw away sensitivity [2].

How to Filter VCF Files by Quality and Depth, Step by Step

The order below matters. Normalization before thresholding. Site filters before genotype filters. Soft filters before hard removal. Verification after everything.

Step 1. Normalize and index

bcftools norm -f reference.fa -m -any -Oz -o normalized.vcf.gz raw.vcf.gz
bcftools index -t normalized.vcf.gz

-f left-aligns and normalizes indels against the reference. -m -any splits multiallelic records into biallelic ones. Expected output looks like this.

Lines   total/split/realigned/skipped:  4,500,000/180,000/12,400/0

If you see a nonzero skipped count, the reference does not match the VCF and every indel position is suspect.

Two consequences of splitting. First, a multiallelic site becomes several records that share the same INFO/DP, so a depth filter sees the same total depth several times. Second, an allele balance filter computed after splitting is per-allele, which is what you want, but a filter computed before splitting sees summed alt counts. Normalize first and keep it that way.

If the file is unsorted, run bcftools sort before indexing. If duplicate records are possible, add bcftools norm -d exact.

Step 2. Summarize the input

Never set a threshold without looking at the distribution first.

bcftools stats -s - normalized.vcf.gz > raw.stats
grep '^SN' raw.stats | head -12
grep '^PSC' raw.stats | cut -f3,4,5,6,10,14 | head

The SN lines give total records, SNP and indel counts, and multiallelic counts. The PSC lines give per-sample homozygous reference, heterozygous, and missing counts, plus average depth in the tenth column. A sample with a very low average depth or a very high missing count is a candidate for removal, not for rescue.

Per-sample depth distribution from FORMAT/DP, including the fraction of sites with depth zero, needs one line of awk.

bcftools query -f '[%SAMPLE\t%DP\n]' normalized.vcf.gz \
  | awk '{n[$1]++; s[$1]+=$2; if($2==0) z[$1]++} \
         END{for(k in n) printf "%s\t%.1f\t%.1f%%\n", k, s[k]/n[k], 100*z[k]/n[k]}' \
  | sort -k2,2n

Expected output is one line per sample with mean depth and percent zero-depth sites.

SAMPLE1	31.4	2.1%
SAMPLE2	28.9	3.0%
SAMPLE3	6.2	41.7%

SAMPLE3 is the problem. Its mean depth is under 10 and 42% of its sites have no reads. No site-level filter will fix that. Drop the sample or accept that its genotypes are mostly missing.

Step 3. Apply site-level filters

Start with a hard site filter using a bcftools expression.

bcftools view -i 'QUAL>=30 && INFO/DP>=10 && INFO/MQ>=40' \
  -Oz -o sites_pass.vcf.gz normalized.vcf.gz
bcftools index -t sites_pass.vcf.gz

Each condition is a floor. QUAL 30 means a 1 in 1,000 chance the site is not variant. INFO/DP 10 means at least ten reads across the whole cohort. INFO/MQ 40 means the reads mapping here have decent confidence.

Three warnings about this line. Records with a missing QUAL (a dot) are dropped by the comparison, so a callset where the caller never wrote QUAL will lose every record. Using a tag the caller did not write, such as INFO/QD on a FreeBayes VCF, produces an error naming the undefined tag, so check the header first with bcftools view -h. And INFO/DP is total depth across all samples, so in a 500-sample cohort at 30x each it will never be near 10, which means the filter only removes sites where nearly every sample is empty. For multi-sample data, the per-sample filters in step 4 do the real work.

Step 4. Apply genotype-level filters

Here you have a real choice. You can mask failing genotypes to missing and keep the site, or you can drop the site entirely.

Masking is usually correct for population data. A locus is often fine in most samples and poor in a few.

bcftools +setGT normalized.vcf.gz -Oz -o masked.vcf.gz -- -t q -n . -e 'FMT/DP<10 || FMT/GQ<20'
bcftools index -t masked.vcf.gz

The same operation with vcftools, which applies genotype-level filters by converting failing genotypes to missing before applying site-level filters:

vcftools --gzvcf normalized.vcf.gz --minDP 10 --minGQ 20 --max-missing 0.9 \
  --recode --recode-INFO-all --stdout | bgzip > masked.vcf.gz
tabix -p vcf masked.vcf.gz

Be clear about what bcftools view -i 'FMT/DP>=10' does. It decides whether the record survives and leaves the failing genotypes in place. It does not repair them. If your goal is to zero out bad genotypes, use +setGT or vcftools.

Step 5. Check allele balance on heterozygotes

Allele balance catches several failure modes at once: mismapped reads, sample swaps, contamination, and heterozygous calls made from a handful of alt reads. Compute it from AD and list the outliers.

bcftools query -f '[%SAMPLE\t%AD\t%GT\n]' masked.vcf.gz \
  | awk -F'\t' '$3=="0/1" {split($2,a,","); n=a[3]+a[4]; \
      if(n==0) next; ab=a[4]/n; \
      if(ab<0.2 || ab>0.8) {c[$1]++}} \
    END{for(k in c) printf "%s\t%d\n", k, c[k]}' \
  | sort -k2,2nr

Expected output is a count of off-balance heterozygous calls per sample. A sample with several thousand of them needs investigation before its variants are believed.

This filter assumes diploid, disomic inheritance and a germline event. It is wrong for polyploids, where the expected alt fraction depends on allele dosage, and wrong for somatic samples, where tumor purity pushes the variant allele fraction well below 0.2. In tetraploid genotyping-by-sequencing data, the theoretical work behind expected genotype quality shows 17 reads are needed to call a simplex heterozygote at 95% accuracy, against 61 reads to resolve allele dosage [5]. A flat 0.2 to 0.8 window applied to that data will delete real calls.

Step 6. Control missingness

After masking, decide how much missing data a site may have.

bcftools view -i 'F_MISSING<=0.2' masked.vcf.gz -Oz -o final_sites.vcf.gz
bcftools index -t final_sites.vcf.gz

For genotyping arrays and low-coverage population data, 0.1 to 0.2 is a reasonable ceiling. For sparse non-model sampling, relax it and report it. The MaizeGDB reanalysis of 1,498 maize accessions filtered explicitly on mapping quality and completeness alongside linkage disequilibrium, which is the same three-axis logic applied at scale [6].

Step 7. Soft filter instead of deleting

For exploratory work, mark first and delete later.

bcftools filter -e 'QUAL<30 || INFO/DP<10 || MQ<40' \
  -s LowQualDepth -m + normalized.vcf.gz -Oz -o soft.vcf.gz
tabix -p vcf soft.vcf.gz

Records that fail get LowQualDepth in the FILTER column. Records that pass keep PASS. The -m + mode appends your label to any existing filter rather than overwriting the caller's verdict, which preserves provenance. You can inspect the marked set, count it, and only then extract the passing records.

bcftools view -f PASS soft.vcf.gz -Oz -o pass_only.vcf.gz

Note that -f PASS keeps only records explicitly marked PASS. A record whose FILTER is . means no filter was applied, which is not the same as passing, and it will be excluded. If you want both, use -f PASS,. and know that you are including unevaluated records.

The trade-off between hard and soft filtering is the single most consequential choice in this workflow. A hard filter is irreversible unless you kept the input. A soft filter costs a little disk and keeps every option open.

Step 8. Hard filters in GATK, if you are in GATK

GATK users should filter SNPs and indels separately because the artifact profiles differ. Select the types first.

gatk SelectVariants -R reference.fa -V raw.vcf.gz --select-type-to-include SNP -O snps.vcf.gz
gatk SelectVariants -R reference.fa -V raw.vcf.gz --select-type-to-include INDEL -O indels.vcf.gz

Then apply the published hard-filter thresholds for SNPs.

gatk VariantFiltration -R reference.fa -V snps.vcf.gz \
  -filter "QD < 2.0" --filter-name "QD2" \
  -filter "QUAL < 30.0" --filter-name "QUAL30" \
  -filter "SOR > 3.0" --filter-name "SOR3" \
  -filter "FS > 60.0" --filter-name "FS60" \
  -filter "MQ < 40.0" --filter-name "MQ40" \
  -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
  -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
  -O snps.filtered.vcf.gz

For indels the strand-bias and read-position thresholds loosen, because indels generate more strand and position artifacts legitimately. The documented values are QD below 2.0, QUAL below 30.0, FS above 200.0, SOR above 10.0, and ReadPosRankSum below minus 20.0.

Extract the passing set with:

gatk SelectVariants -R reference.fa -V snps.filtered.vcf.gz --exclude-filtered -O snps.pass.vcf.gz

Step 9. Record what you did

Write the exact command line, the tool version, and the input file checksum into the VCF header. FiNGS was built on the argument that ad hoc post-calling filtering is unreproducible unless the filters are explicit, configurable, and shared [7]. The same principle applies whether you use FiNGS or a shell script. bcftools annotate -h notes.hdr appends header lines. A three-line provenance block beats a forgotten threshold six months later.

Threshold options and what each one buys

The table below collects the starting values worth trying. The GATK values are the tool's published hard-filter defaults, not universal truths. The depth and genotype values are common germline starting points.

A depth floor of 15 for germline work is defensible. In a whole-genome cohort where technical repeats were compared at DNA damage response loci, concordance was near-perfect after restricting to nonsilent SNVs with read depth of at least 15, with discordance concentrated in the low-QC fraction [8].

Checking that the filter did what you think

Verification is not optional. Run these five checks on the before and after files.

Variant counts by class. Re-run bcftools stats and compare SNP, indel, and multiallelic counts. A filter that removes 30% of SNPs and 70% of indels has a different problem than one that removes 30% of each.

Transition to transversion ratio. Ti/Tv should be similar before and after. A sharp drop means you are preferentially deleting true transitions, which is what over-filtering looks like. A sharp rise can mean you are keeping only easy, well-covered sites.

Het to hom ratio per sample. Compare the per-sample heterozygous counts from bcftools stats. A sample whose ratio collapses after filtering has lost most of its heterozygous calls to a depth threshold, which biases every downstream allele frequency estimate.

Depth and missingness after filtering. Repeat the per-sample awk summary from step 2 on the final file. If median depth at retained sites is now 40 when the raw median was 32, you have selected high-coverage sites and biased the callset toward easy regions.

Downstream sensitivity. Sequential filtering changes the site frequency spectrum, and that change propagates into demographic inference. In the burrowing owl resource built with snpArcher, the authors show how the sequence of filtering steps shifts the site frequency spectrum and downstream demographic results across 137 genomes, and their QC dashboard aggregates principal component analysis, relatedness, depth-missingness diagnostics, and admixture to catch batch effects before interpretation [4]. If your filtered callset feeds a demographic model or a GWAS, re-run the model and confirm the conclusion survives the filter.

A visual QC report helps here. VariantQC aggregates variants by dataset, chromosome, sample, and filter type in an interactive report designed to flag outliers, and it reads standard VCF so it slots into an existing pipeline [9]. A table of numbers you never look at is not QC.

Common Mistakes and Limitations

Filtering INFO/DP as if it were per-sample depth. In a cohort, INFO/DP sums across samples. A floor of 10 is meaningless at 200 samples. Use FORMAT/DP for per-sample decisions and reserve INFO/DP for detecting empty loci.

Assuming QUAL says something about a genotype. QUAL is about the site. A site with QUAL 500 can still contain a sample with GQ 3. The Germline VCF Annotator pipeline makes this distinction concrete by attaching read-evidence metrics to each variant and assigning a rules-based QC class, rather than trusting a single site score [8].

Treating a dot in FILTER as PASS. A dot means no filter was applied. It is not a verdict. Use -f PASS deliberately and know what you are excluding.

Filtering before normalization. Multiallelic representation changes how depth and allele counts are interpreted, and unnormalized indels sit at positions the reference does not agree with. Normalize first, always.

Filtering before joint genotyping. Depth is recomputed at joint calling. A per-sample depth filter applied to single-sample GVCFs or raw calls will not carry over correctly. Apply depth filters to the final multi-sample callset.

Running VQSR on a small cohort. Variant Quality Score Recalibration needs a large variant set and truth resources to fit its model, and it is known to fail to run in some settings. VEF was developed partly because VQSR and hard filtering are user-dependent and unreliable in exactly those cases, and it demonstrates that a supervised model trained on gold-standard variants beats both on human whole-genome data [10]. If you have no truth set and a few dozen samples, use hard filters and say so.

Applying GATK thresholds to non-GATK callers. FreeBayes, DeepVariant, and bcftools mpileup emit different INFO tags and different QUAL semantics. A threshold tuned for one caller is a guess for another. VEF's authors explicitly restrict training and testing to VCFs generated by the same tools, because feature characteristics must match [10].

Discarding caller-flagged variants without looking. Some flagged variants are real. Re-examining them is the entire purpose of vcfView, which lets variants that failed caller filters be reviewed to recover sensitivity [2].

Ignoring the polyploid case. SWEEP was built for polyploid crops because standard filters fail when homeologous sequence creates apparent polymorphisms, and it filters on depth, alternate allele ratio, and SNP quality on top of its own haplotype contrast [11]. If your organism is polyploid, default filters are wrong by construction.

Opening a VCF in a spreadsheet. Raw VCF is designed for programmatic parsing. Spreadsheet import distorts values through automatic type conversion, and the resulting table is not the file you started with [8]. Use bcftools query to extract columns and keep the VCF as the source of truth.

Deleting intermediates. Once you hard filter, the removed records are gone. Keep the raw VCF and the normalized VCF, both indexed, for at least the life of the project.

Dropping low-coverage samples late. The snpArcher protocol places sample-level QC, including depth-missingness diagnostics, before downstream analysis rather than after, because outlier samples left in place distort principal components and relatedness estimates [4]. Remove or flag them early.

Over-filtering rare variants in sparse cohorts. A missingness ceiling of 0.1 is reasonable in a densely sequenced cohort and destructive in a scattered sampling of wild individuals. Choose the ceiling from the sampling design, not from a habit.

Filtering by read depth when the reference is the problem. A region with collapsed paralogs will show high depth, low MAPQ, and elevated FS together. Any one of those filters catches it. Filtering only on depth lets it through.

Limitations. No hard filter recovers a variant that the caller never proposed. Thresholds calibrated on one species, one library prep, or one sequencing depth do not transfer cleanly, and non-human species often lack the truth sets and population databases that make supervised filtering possible, which is why tools like variantscanR were written specifically for domestic animals and their limited reference resources [12]. Somatic call sets need their own logic, since tumor purity and subclonality break the diploid germline assumptions behind allele balance, and dedicated somatic filters such as FiNGS address that case with configurable, reproducible rules [7]. Filtering is also a sequence, not a single step. OVAS chains quality control and inheritance filters in a defined order for exactly this reason [13]. Finally, if the VCF comes from a clinical case in an animal, the filtered list is a research artifact, and a veterinarian should interpret it alongside the animal's history and phenotype, not as a diagnosis on its own.

Frequently Asked Questions

What read depth should I require when I filter a VCF?

A per-sample depth of at least 10 is a common germline floor, and 15 or more is safer when genotype concordance matters. Use FORMAT/DP for per-sample decisions and keep INFO/DP as a separate check for empty loci.

Should I filter on QUAL or on GQ?

Filter on both, at different levels. QUAL decides whether a site survives, and GQ decides whether an individual genotype survives. A high QUAL site can still hold a low-GQ genotype.

Is INFO/DP the same as FORMAT/DP?

No. INFO/DP is total depth across all samples at a site, while FORMAT/DP is the depth for one sample. Some callers also compute them differently, so pick one deliberately and state which you used.

Should I hard-filter or soft-filter a VCF?

Soft-filter first whenever you might need the removed records. Soft filtering marks records in the FILTER column and leaves the data intact, while hard filtering deletes them permanently unless you kept the input file.

Can I use VQSR instead of hard filters on a small cohort?

Usually not. VQSR needs a large variant set and truth resources to fit its model, and it can fail to run on small or non-model datasets. Hard filters, or a supervised approach trained on a matching gold standard, are the practical choice.

Why did my variant count collapse after filtering?

The most common causes are a depth filter that is too high for the cohort's actual coverage, a QUAL filter applied to a callset where QUAL is missing, or an undefined INFO tag that causes every record to fail the expression. Check the header and the depth distribution before changing thresholds.

Do I need the reference FASTA to filter a VCF?

You need it to normalize, which should happen before thresholding. Pure threshold filters on QUAL, DP, and GQ can run without it, but you lose left-alignment and indel normalization, and split multiallelic records may sit at inconsistent positions.

What is the difference between FASTQ and FASTA?

FASTQ stores raw reads with a per-base Phred quality string, and FASTA stores sequence only with no quality values. References are FASTA, reads are FASTQ, and the base qualities in FASTQ are the origin of the quality scores you later filter on in the VCF.

Sources

1. VariFAST: a variant filter by automated scoring based on tagged-signatures.
2. vcfView: An Extensible Data Visualization and Quality Assurance Platform for Integrated Somatic Variant Analysis.
3. 123VCF: an intuitive and efficient tool for filtering VCF files.
4. Variant calling in nonmodel organisms with snpArcher.
5. Expected Genotype Quality and Diploidized Marker Data from Genotyping-by-Sequencing of Urochloa spp. Tetraploids.
6. A unified VCF dataset from nearly 1,500 diverse maize accessions and resources to explore the genomic landscape of maize.
7. FiNGS: high quality somatic mutations using filters for next generation sequencing.
8. Germline VCF Annotator: a lightweight pipeline for processing germline VCFs with robust variant extraction and read evidence quality control.
9. VariantQC: a visual quality control report for variant evaluation.
10. VEF: a variant filtering tool based on ensemble methods.
11. SWEEP: A Tool for Filtering High-Quality SNPs in Polyploid Crops.
12. VariantscanR: an R-package as a clinical tool for variant filtering of known phenotype-associated variants in domestic animals.
13. OVAS: an open-source variant analysis suite with inheritance modelling.