Filtering a VCF Down to PASS Variants (and What to Do Next)

By the end of this you will have three artifacts from one raw VCF: a PASS-only, left-aligned, biallelic VCF with a .tbi index, a TSV of the fields you care about (CHROM, POS, REF, ALT, genotype, depth, genotype quality, allele balance), and a set of counts that tell you whether the filtering did something sane. You need bcftools 1.18 or newer (it bundles tabix and bgzip via htslib), the exact reference FASTA your VCF was called against (usually GRCh38 with or without alt contigs), and about 3x the size of your input VCF in free disk. A single-sample WGS VCF from a 30x genome runs 1–2 GB gzipped and holds roughly 4.5–5 million variants before filtering. Everything below is command line. If your goal is “open this in Excel,” skip to step 7, but read step 2 first, because the thing most people call filtering is a grep that quietly throws away variants they wanted.

1. Read the header before you filter anything

The FILTER column is not standardized across callers. Every value in it is defined in the header of your own file, and that is the only authoritative description of what it means.

bcftools view -h sample.vcf.gz | grep '^##FILTER'

For a GATK HaplotypeCaller VCF run through VQSR you will see something like:

##FILTER=<ID=PASS,Description="All filters passed">
##FILTER=<ID=VQSRTrancheSNP99.90to100.00,Description="Truth sensitivity tranche level ...">
##FILTER=<ID=LowQual,Description="Low quality">

For DeepVariant:

##FILTER=<ID=RefCall,Description="Genotyping model thinks this site is reference.">
##FILTER=<ID=LowQual,Description="Confidence in this variant being real is below calling threshold.">

For DRAGEN you get DRAGENSnpHardQUAL, DRAGENIndelHardQUAL, lod_fstar, and sample-level filters in FORMAT/FT. Illumina’s Strelka adds LowGQX, HighDPFRatio, SiteConflict.

Then get the actual distribution, not the header’s list of possibilities:

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

A typical single-sample DeepVariant WGS run gives you something near 4.4M PASS and 1.5M RefCall. If PASS is a small minority of your rows, you have either a joint-called cohort VCF where most sites are not variant in your sample, or a caller whose output includes non-variant reference blocks. Check which before you proceed.

2. Know what PASS guarantees

PASS means “this record cleared the filters the caller applied,” and nothing more. It is a statement about the site, not about your sample’s genotype at that site, and not about whether the variant means anything.

Two consequences matter in practice. First, a record can be PASS and your sample’s genotype at it can still be ./. or supported by four reads. In a multi-sample VCF the record-level FILTER is shared across all samples while FORMAT/FT is per sample. Second, callers differ in how aggressively they filter, and the FILTER column carries no information about the consequence, frequency, or reproducibility of a variant. Orthogonal evidence helps: k-mer based approaches score each variant against the raw read k-mer spectrum and catch false positives that pass the caller’s own internal model, which is useful when you care more about precision than recall ¹.

Structural variants are a separate problem. SV callers emit records with SVTYPE, END, and imprecise breakpoints, and their FILTER semantics have nothing in common with SNV filters. Large family studies that estimate de novo SV rates apply heavy manual and automated curation on top of caller PASS, because raw SV call sets carry high false positive rates ². If your VCF has an ALT column full of <DEL> and <DUP>, handle it in a separate pipeline.

3. Filter to PASS with bcftools, not grep

The common mistake is grep PASS file.vcf. That matches the string PASS anywhere on the line, including inside INFO fields like CSQ=...PASS... or a sample name, and it also matches header lines. It will also miss records where FILTER is . (missing), which some tools emit for unfiltered-but-valid sites.

Use the FILTER-aware flag:

bcftools view -f PASS -Oz -o sample.pass.vcf.gz sample.vcf.gz
bcftools index -t sample.pass.vcf.gz

-f, --apply-filters understands the FILTER column structure. If you want to keep unfiltered records too:

bcftools view -f PASS,. -Oz -o sample.pass.vcf.gz sample.vcf.gz

If you need the inverse, or want everything except one specific filter:

# everything that is not PASS
bcftools view -i 'FILTER!="PASS"' -Oz -o sample.fail.vcf.gz sample.vcf.gz

# keep PASS plus records whose only failure was LowGQX
bcftools view -i 'FILTER="PASS" || FILTER="LowGQX"' -Oz -o sample.lenient.vcf.gz sample.vcf.gz

Note the difference between -f and -i 'FILTER=...': with -i, FILTER="X" means the FILTER column contains X among possibly several values, while FILTER=="X" (double equals) requires X to be the only value. This distinction is easy to miss and changes counts on multi-filter records.

Always write -Oz (bgzipped VCF) and index. Uncompressed VCF costs you random access, and every downstream tool will want a .tbi or .csi.

4. Normalize before you do anything else

Filtering is not the first operation you should run. Normalization is, because unnormalized indels break every position-based join you will later do against gnomAD, ClinVar, or your own annotation tables.

bcftools norm -m -any -f GRCh38.fa -c w -Oz -o sample.norm.vcf.gz sample.pass.vcf.gz
bcftools index -t sample.norm.vcf.gz

• -m -any splits multiallelic records into one row per ALT allele. A record with ALT=A,AT becomes two rows. This is what you want for any per-variant analysis, and it is what annotation databases assume.
• -f GRCh38.fa left-aligns and trims indels against the reference.
• -c w warns when the REF allele in the VCF does not match the reference FASTA instead of silently exiting. If you get many of these, you have the wrong reference build or the wrong contig naming (chr1 vs 1).

Splitting multiallelics changes FORMAT/AD and INFO/AC semantics. bcftools norm handles the standard fields, but custom INFO fields added by an annotator may not be split correctly. Normalize first, annotate second.

For canonical variant descriptions (HGVS NM_...:c.76A>T style) you need transcript-aware mapping, and the edge cases around intronic offsets and reference discrepancies are worse than they look. Use a validator rather than string-building HGVS yourself ³.

5. Add genotype-level filters on top of PASS

PASS says the site is real. Now say something about your genotype at it. For a single sample WGS at 30x:

bcftools filter \
  -S . \
  -e 'FMT/DP<10 | FMT/GQ<20' \
  -Oz -o sample.gtfilt.vcf.gz sample.norm.vcf.gz

-S . sets failing genotypes to missing rather than dropping the record. Then drop records where your sample now has no genotype:

bcftools view -e 'GT="mis"' -Oz -o sample.clean.vcf.gz sample.gtfilt.vcf.gz
bcftools index -t sample.clean.vcf.gz

Thresholds we use as a default on 30x WGS: DP >= 10, GQ >= 20, and for heterozygous calls an allele balance between 0.25 and 0.75. Allele balance needs the AD field:

bcftools view -i 'GT="het" & (FMT/AD[0:1])/(FMT/DP) > 0.25 & (FMT/AD[0:1])/(FMT/DP) < 0.75' \
  -Oz -o sample.het.vcf.gz sample.clean.vcf.gz

The AD[0:1] indexing is sample 0, allele index 1 (first ALT). After splitting multiallelics this is the right index. Before splitting it is not, which is another reason to normalize first.

Raise DP if you have deeper coverage. Homozygous-alt calls tolerate lower thresholds than hets, because a het at DP 10 with a 2/8 split is much more likely to be an artifact. Also keep in mind that these thresholds interact with your caller: pipelines that combine calling and filtering in a single fast pass make different precision/recall tradeoffs than staged ones, and their PASS sets are not interchangeable ⁴.

6. Subset to regions you care about

Genome-wide filtering leaves you 4M+ rows. Most questions are about a gene list or a panel.

# BED of regions, 0-based half-open
bcftools view -R genes.bed -Oz -o sample.panel.vcf.gz sample.clean.vcf.gz

Use -R (requires an index, does random access, fast) rather than -T (streams the whole file) when the region set is small relative to the genome. -T is faster when you are subsetting to a large fraction of sites. On very large cohort VCFs the cost of repeated full passes dominates, and execution engines that skip untouched blocks rather than streaming everything exist for exactly this reason ⁵.

7. Flatten to TSV, and think twice about Excel

bcftools query is the tool. It prints exactly the fields you name, one line per record:

bcftools query \
  -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT\t%DP\t%GQ\t%AD]\n' \
  -H sample.clean.vcf.gz > sample.tsv

The square brackets loop over samples. -H prints a header line, though it prefixes column names with # [1] style indices that you will want to clean up. To convert to CSV, do it in the same step rather than with sed on the TSV, since ALT fields and annotation strings contain commas:

bcftools query -f '%CHROM,%POS,%REF,%ALT,[%GT,%DP]\n' sample.clean.vcf.gz > sample.csv

If you added VEP annotations, the CSQ INFO field is a pipe-delimited, comma-separated multi-transcript blob. Use bcftools +split-vep to pull named subfields into columns:

bcftools +split-vep sample.vep.vcf.gz \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%gnomADg_AF\t[%GT]\n' \
  -d -A tab > sample.annot.tsv

-d expands one line per transcript consequence. Drop it if you only want the picked consequence.

On Excel: it caps at 1,048,576 rows, so a WGS PASS set does not fit. It also autoconverts gene symbols like MARCH1 and SEPT9 into dates unless you explicitly import every column as text via the Get Data dialog. Excel-based genomic tools exist and are genuinely useful for small, curated sets, mitochondrial variant panels being one example where the whole call set fits comfortably in a worksheet ⁶, and haplotype tables for a locus being another ⁷. For a filtered gene panel of a few hundred rows Excel is fine. For anything genome-scale, load the TSV into DuckDB:

duckdb -c "CREATE TABLE v AS SELECT * FROM read_csv_auto('sample.tsv', delim='\t'); \
           SELECT count(*), filter FROM v GROUP BY filter;"

DuckDB reads a 5M-row TSV in a few seconds and lets you join against gnomAD Parquet files without writing a pipeline. That is the workflow we would use. If you want a GUI over annotated variants without writing SQL, desktop applications built for exactly this filter-summarize-visualize loop are a reasonable middle ground ⁸.

8. Prioritize, then sanity-check

Once you have a clean PASS set, ranking matters more than filtering. Missense predictors give a per-variant score for whether an amino acid substitution is tolerated, precomputed for whole genomes so you do not have to run alignments yourself ⁹. Phenotype-driven tools go further and combine variant-level pathogenicity with the similarity between a patient’s observed phenotype and known gene-disease associations, which is the approach used in diagnostic exome and genome analysis ¹⁰. These produce ranked candidate lists. They do not produce findings. Any variant you are tempted to act on needs confirmation in a clinical laboratory and interpretation by a clinician or genetic counselor.

Finally, check your work with counts and ratios rather than eyeballing rows:

bcftools stats sample.clean.vcf.gz > stats.txt
grep -E '^SN|^TSTV' stats.txt

For a single human WGS after PASS and genotype filtering, expect roughly 4.0–4.5M SNVs, 700k–900k indels, a genome-wide transition/transversion ratio near 2.0–2.1, and a het/hom-alt ratio around 1.5–2.0. A Ti/Tv below 1.8 genome-wide usually means false positives survived. A Ti/Tv above 2.3 genome-wide usually means you filtered into a biased subset. On exome-restricted calls Ti/Tv near 3.0 is normal because coding sites are enriched for transitions, so compare like with like.

Common problems

grep PASS returned fewer variants than expected. It matched header lines and INFO strings and missed . filters. Use bcftools view -f PASS,..

“The sequence does not match the reference” from bcftools norm. Wrong reference FASTA or wrong contig naming. Compare bcftools view -h file.vcf.gz | grep contig against grep '>' ref.fa. Renaming contigs with bcftools annotate --rename-chrs fixes naming but not a genuine build mismatch (GRCh37 coordinates against a GRCh38 FASTA), which requires a liftover.

Every record is PASS and there are 30 million of them. You have a gVCF with reference blocks, not a VCF. Look for END= in INFO and <NON_REF> in ALT. Run gatk GenotypeGVCFs first, or filter with bcftools view -e 'ALT="<NON_REF>"' as a rough pass.

FILTER column is empty for every record. The caller never applied filters. PASS filtering is a no-op here and you need to build your own criteria from QUAL, DP, and QD.

A variant you expect is absent from the PASS set. Check whether it is in the unfiltered VCF at all. If it is there with a non-PASS filter, look up that filter’s description in the header. If it is not there at all, the site may have had no coverage, which is a different problem from a failed filter. Check with samtools depth -r chr7:117559590-117559600 sample.bam.

Annotation columns misaligned after splitting multiallelics. You annotated before normalizing. Redo it in the order: norm, then annotate, then query.

Sample-level filters ignored. FORMAT/FT exists in some VCFs and bcftools view -f does not look at it. Filter it explicitly with -i 'FMT/FT="PASS"'.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Giulio Formenti, Arang Rhie, Brian P. Walenz, et al. Merfin: improved variant filtering, assembly evaluation and polishing via k-mer validation. Nature Methods, 2022. https://doi.org/10.1038/s41592-022-01445-y ↩

2. Jonathan R. Belyeu, Harrison Brand, Harold Wang, et al. De novo structural mutation rates and gamete-of-origin biases revealed through genome sequencing of 2,396 families. The American Journal of Human Genetics, 2021. https://doi.org/10.1016/j.ajhg.2021.02.012 ↩

3. Peter J. Freeman, Reece K. Hart, Liam J. Gretton, et al. VariantValidator: Accurate validation, mapping, and formatting of sequence variation descriptions. Human Mutation, 2017. https://doi.org/10.1002/humu.23348 ↩

4. Colby Chiang, Ryan M Layer, Gregory G Faust, et al. SpeedSeq: ultra-fast personal genome analysis and interpretation. Nature Methods, 2015. https://doi.org/10.1038/nmeth.3505 ↩

5. Ehsan Estaji, Jian-Feng Mao. VariantFlow: a selective-execution engine for efficient population genomic computation on large variant datasets. 2026. https://doi.org/10.64898/2026.09.01.748643 ↩

6. Jonathan L. King, Antti Sajantila, Bruce Budowle. mitoSAVE: Mitochondrial sequence analysis of variants in Excel. Forensic Science International: Genetics, 2014. https://doi.org/10.1016/j.fsigen.2014.05.013 ↩

7. Cong Feng, Xingwei Wang, Shishi Wu, et al. HAPPE: A Tool for Population Haplotype Analysis and Visualization in Editable Excel Tables. Frontiers in Plant Science, 2022. https://doi.org/10.3389/fpls.2022.927407 ↩

8. Jing Zhai, Nelly Pitteloud, Federico A. Santoni. GenMasterTable: a user-friendly desktop application for filtering, summarising, and visualising large-scale annotated genetic variants. BMC Bioinformatics, 2025. https://doi.org/10.1186/s12859-025-06238-6 ↩

9. Robert Vaser, Swarnaseetha Adusumalli, Sim Ngak Leng, et al. SIFT missense predictions for genomes. Nature Protocols, 2015. https://doi.org/10.1038/nprot.2015.123 ↩

10. Damian Smedley, Julius O B Jacobsen, Marten Jäger, et al. Next-generation diagnostics and disease-gene discovery with the Exomiser. Nature Protocols, 2015. https://doi.org/10.1038/nprot.2015.124 ↩