| name | bio-vcf-statistics |
| description | Compute and interpret VCF quality-control metrics (Ti/Tv, het/hom, novel/known, missingness, HWE, contamination, relatedness) with bcftools stats, vcftools, plot-vcfstats, and identity tools (somalier, peddy, KING). Use when judging whether a callset is trustworthy, diagnosing a low Ti/Tv or outlier het/hom sample, deciding whether an HWE deviation is error or biology, screening a cohort for sample swaps/contamination/wrong-sex before analysis, or comparing call sets before and after filtering. Not for applying filters (see variant-calling/filtering-best-practices) or normalizing representation (see variant-calling/variant-normalization). |
| tool_type | cli |
| primary_tool | bcftools |
Version Compatibility
Reference examples tested with: bcftools 1.19+, vcftools 0.1.16+, somalier 0.2.19+, peddy 0.4.8+, cyvcf2 0.30+, numpy 1.26+
Before using code patterns, verify installed versions match. If versions differ:
- Python:
pip show <package> then help(module.function) to check signatures
- CLI:
<tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Note: bcftools gtcheck was rewritten around bcftools 1.10; the old -G/--GTs-only flag is gone. Modern gtcheck cross-checks all samples in one file when -g is omitted, and uses -E/--error-probability. Verify flags against the installed build.
VCF Statistics
"Is this callset any good, and are the samples who the manifest says they are?" -> Summarize variant counts and quality distributions, then read each QC metric as a signal whose expected value is context-dependent.
- CLI:
bcftools stats (+ plot-vcfstats), vcftools, somalier, peddy
- Python:
cyvcf2 for custom per-record statistics
The governing principle
No QC metric has a universal pass value. Each metric's expected range depends on the assay (WGS vs WES), the sample's ancestry, and the cohort it sits in, so QC is not threshold-checking but reading a deviation for its mechanistic meaning and acting on it. Three rules follow. First, always compare a sample against a matched cohort (same assay, same inferred ancestry), never against an absolute number. Second, order of operations matters: apply genotype-level filters (set ./. where GQ/DP/allele-balance fail) BEFORE computing cohort missingness, call rate, or HWE, because low-quality genotypes left in the matrix drive spurious missingness and HWE deviation. Third, sample-identity QC is a graph problem, not a per-sample check: build the all-pairs relatedness matrix and reconcile it against the declared pedigree/manifest, and do this BEFORE any association or burden analysis, because one undetected swap or contaminated sample can fabricate or erase a genome-wide-significant hit.
QC metric decision table
Expected values are stated conventions (they shift with capture kit, ancestry, reference build, and caller); treat them as starting points and confirm against a matched cohort.
| Metric | Expected (WGS) | Expected (WES) | A deviation MEANS | Action |
|---|
| Ti/Tv (overall) | ~2.0-2.1 | ~3.0-3.3 | Low -> false-positive transversions dilute the signal (random errors have Ti/Tv ~0.5); high -> over-filtering removed transversions | Tighten site filters if low; check filter for transversion bias if high |
| Ti/Tv (novel only) | near overall | near overall | Novel Ti/Tv well below the known-site value -> FP contamination in the novel fraction | Raise stringency; the novel set is where FPs concentrate |
| Het/hom ratio | ancestry-dependent (~1.5-1.6 EUR, ~2.0+ AFR) | same | High vs same-ancestry cohort -> contamination or reference bias; low -> inbreeding/consanguinity/ROH or chromosome loss | Stratify by ancestry first; then flag within-ancestry outliers |
| Novel fraction (vs dbSNP) | low for common; high for rare/singleton | same | Novel COMMON variants -> almost always artifacts; high novel rare -> expected or under-studied population | Stratify novel% by frequency; investigate common novels |
| Call rate (per sample) | >95-98% | >95-98% | Low -> low-coverage/low-quality sample | Drop worst samples, then recompute per-variant missingness (iterate) |
| Call rate (per variant) | >95-99% | >95-99% | Low -> site in a hard-to-genotype region | Drop worst variants after sample QC; rare variants tolerate more missing |
| Het allele-balance (hets) | centered on 0.5 | 0.5 | Shifted away from 0.5 + elevated het count -> contamination signature | Confirm with VerifyBamID2 (on the BAM) or CHARR (on the VCF/gVCF) -- the real test |
| Excess-het / HWE | in equilibrium within ancestry | same | EXCESS het -> collapsed paralog/CNV mapping artifact; het DEFICIT -> often real (Wahlund/inbreeding) | Filter on excess het only, within ancestry, in controls; do NOT blanket-filter HWE |
| Relatedness (kinship) | matches manifest | matches manifest | Unexpected high kinship -> sample swap/duplicate; a "replicate" that is not -> mislabel | Reconcile the all-pairs matrix against the pedigree with KING/somalier/peddy |
bcftools stats
Goal: Generate comprehensive variant statistics (counts, Ti/Tv, per-sample het/hom and singletons, indel and depth distributions).
Approach: Run bcftools stats and read the section-tagged output lines; add -s - for per-sample metrics.
bcftools stats input.vcf.gz > stats.txt
bcftools stats -s - input.vcf.gz > per_sample.txt
bcftools stats file1.vcf.gz file2.vcf.gz > cmp.txt
plot-vcfstats -p qc_plots/ stats.txt
Output sections: SN summary numbers, TSTV transition/transversion, SiS singletons, AF allele-frequency spectrum, QUAL quality distribution, IDD indel-length distribution, ST substitution types, DP depth distribution, PSC per-sample counts (hom-ref, het, hom-alt, transitions, transversions, missing), PSI per-sample indels.
bcftools stats input.vcf.gz | grep "^SN" | cut -f3-
bcftools stats input.vcf.gz | grep "^TSTV" | cut -f5
bcftools stats -s - input.vcf.gz | grep "^PSC"
Ti/Tv ratio
Transitions (purine<->purine A<->G, pyrimidine<->pyrimidine C<->T) are favored over transversions because CpG deamination (methylated C->T) is the single most common vertebrate point mutation and is a transition, and because transitions are more often synonymous. A random error spectrum gives Ti/Tv ~0.5, so a callset diluted with false positives drifts DOWNWARD. WES runs higher than WGS (~3.0-3.3 vs ~2.0-2.1) because coding regions are CpG- and constraint-enriched and transitions at CpGs plus codon degeneracy push the ratio up.
Ti/Tv is the single fastest gross-error smell test. A WES callset reporting Ti/Tv ~2.1 is telling the analyst the filtering is too loose. Stratify: compute Ti/Tv on the novel fraction separately (below) since that is where false positives concentrate. Ancestry and target design shift the exact number, so compare against a matched cohort, not the absolute.
Het/hom ratio and ancestry
The per-sample het:non-ref-hom ratio reflects heterozygosity relative to the reference, so it is strongly ancestry-dependent: African-ancestry genomes diverge more from GRCh (a mostly European-ancestry assembly) and carry more het calls (commonly ~2.0+), European ancestry sits ~1.5-1.6, and the value shifts across populations. A single global het/hom cutoff is therefore wrong: it would flag every AFR sample in a EUR-tuned pipeline.
Infer ancestry first (peddy/somalier project onto 1000 Genomes PCs), then flag outliers WITHIN each ancestry group. An elevated het/hom versus same-ancestry peers indicates contamination (foreign reads manufacture spurious hets) or reference bias; a depressed het/hom indicates inbreeding/consanguinity, a long run of homozygosity, or chromosome loss.
Novel/known ratio via dbSNP
Overlapping the callset with dbSNP gives an orthogonal false-positive signal independent of Ti/Tv. Annotate known/novel and stratify by frequency, because the diagnostic differs by allele frequency: novel COMMON variants are almost always artifacts (real common variants are already catalogued), while novel rare/singleton variants are expected and biologically real. A single novel% without frequency stratification is uninformative.
bcftools annotate -a dbsnp.vcf.gz -c ID input.vcf.gz -Oz -o annotated.vcf.gz
bcftools view -H annotated.vcf.gz | awk '{n++; if($3==".") novel++} END{print "novel:", novel/n}'
Missingness and call rate
Call rate is the fraction of sites with a non-missing genotype for a sample; missingness is its complement. Only ./. (no-call) counts as missing; 0/0 (confident hom-ref) does NOT, so treating no-call as hom-ref biases allele frequencies. GWAS convention drops samples below ~95-98% call rate and variants below ~95-99% (dataset-dependent conventions).
The iteration trap: sample-level and variant-level missingness are coupled, so applying both thresholds in one pass is wrong. Drop the worst samples, recompute per-variant missingness, drop the worst variants, and repeat. Critically, apply genotype-level filters (./. where GQ<20 or DP<8) BEFORE computing cohort missingness or HWE, or low-quality genotypes will drive both.
vcftools --gzvcf input.vcf.gz --missing-indv --out sample_miss
vcftools --gzvcf input.vcf.gz --missing-site --out site_miss
HWE filtering (excess-het only)
Hardy-Weinberg testing flags genotype frequencies that deviate from p^2:2pq:q^2, but a naive two-sided HWE gate throws away real biology. The discipline:
- Filter on EXCESS heterozygosity only. Heterozygote excess is the signature of a mapping artifact: a duplicated/collapsed region piles reads from two paralogous copies onto one locus, manufacturing spurious hets everywhere. Heterozygote DEFICIT, by contrast, is often real (Wahlund effect from substructure, inbreeding, selection, a true null allele). GATK's
ExcessHet and InbreedingCoeff and gnomAD's ExcessHet filter the excess side only.
- Compute within an ancestry-homogeneous subgroup. Pooling populations with different allele frequencies induces a Wahlund het-deficit that mimics genotyping error; HWE on a mixed cohort filters good variants.
- Use the EXACT test, not chi-square. The chi-square approximation is anticonservative for rare variants and small samples; the standard is the Wigginton-Cutler-Abecasis exact test (implemented in vcftools
--hardy and PLINK).
- In case/control studies, compute HWE in CONTROLS only. A true association at a strong-effect locus produces HWE deviation in cases; filtering it removes the hit.
vcftools --gzvcf controls.vcf.gz --hardy --out hwe
bcftools query -f '%CHROM\t%POS\t%INFO/ExcessHet\n' input.vcf.gz | awk '$3>54.69'
The threshold 54.69 is GATK's default ExcessHet cutoff (phred-scaled p ~= 3.4e-6, ~the 1000-sample z=-4.5 boundary); tune it to the cohort size.
Contamination signatures
Cross-sample contamination is visible in VCF statistics before any dedicated test: the het allele-balance distribution shifts away from 0.5 (foreign reads add minor-allele support at true hom sites and skew true hets), the het count and het/hom ratio rise, and the novel-fraction Ti/Tv drops. These are SIGNALS, not the measurement. The real test runs on the BAM/CRAM: VerifyBamID2 (Zhang et al. 2020) estimates the contamination fraction alpha ancestry-agnostically by modeling observed allele fractions against population frequencies, and CHARR estimates alpha directly from VCF-level reference-read counts at hom-alt sites. An alpha above ~0.02-0.03 is a red flag; somatic pipelines are sensitive to even 1%. GATK pipelines feed --contamination alpha from VerifyBamID2. See variant-calling/gatk-variant-calling for wiring contamination estimates into calling.
bcftools query -i 'GT="het"' -f '[%AD]\n' input.vcf.gz | \
awk -F',' '{ab=$2/($1+$2); s+=ab; n++} END{print "mean het AB:", s/n}'
Sample-swap, relatedness, and sex checks (mandatory cohort QC)
Sample swaps are among the most common errors in sequencing studies, so identity QC is not optional. Build the all-pairs relatedness matrix and reconcile it against the manifest; swaps appear as off-diagonal surprises (unexpected relatedness) or on-diagonal failures (a "replicate" that is not).
| Tool | Input | Detects | Notes |
|---|
bcftools gtcheck | VCF | Same-file swaps/duplicates (quick) | Native, no reference panel; discordance score, not a relatedness graph |
| KING (Manichaikul 2010) | PLINK bed | Robust kinship without allele-freq/ancestry assumptions | Reference method; kinship bands below |
| peddy (Pedersen 2017) | VCF + PED | Reported vs inferred sex, relationships, ancestry (PCA on 1000G) | Fast, VCF-only; ideal PED-vs-VCF reconciliation |
| somalier (Pedersen 2020) | BAM/CRAM/VCF | Relatedness/ancestry/sex at scale; cross-checks RNA-seq vs WGS | Tiny per-sample sketches; tens of thousands of samples in seconds |
KING kinship coefficient bands: >0.354 duplicate/MZ twin, [0.177, 0.354] first-degree (parent-child, full sib), [0.0884, 0.177] second-degree, [0.0442, 0.0884] third-degree. A pair not expected to be related at ~0.5 is a swap or duplicate.
bcftools gtcheck input.vcf.gz > gtcheck.txt
bcftools gtcheck -g reference.vcf.gz query.vcf.gz
python -m peddy -p 4 --plot --prefix cohort_qc input.vcf.gz cohort.ped
somalier extract -d extracted/ --sites sites.vcf.gz -f ref.fa input.vcf.gz
somalier relate --ped cohort.ped extracted/*.somalier
somalier ancestry --labels 1kg-labels.tsv 1kg/*.somalier ++ extracted/*.somalier
vcftools --gzvcf input.vcf.gz --relatedness2 --out kin
Stratified evaluation
A caller with 99% overall accuracy may drop to 70% in difficult regions, so single-number accuracy hides where a callset fails. Evaluate stratified by region class using GIAB stratification BED files (github.com/genome-in-a-bottle/genome-stratifications).
bcftools stats -R easy_regions.bed input.vcf.gz > easy.txt
bcftools stats -R difficult_regions.bed input.vcf.gz > difficult.txt
Key strata and their failure modes: homopolymer runs (systematic indel errors, Illumina/Ion Torrent), tandem repeats / low-complexity (alignment ambiguity inflates FP and FN), segmental duplications (paralogous mapping -> false hets), high-GC >70% / low-GC <25% (coverage-bias missingness), MHC/centromeric (extreme polymorphism or repetitiveness). A drop in Ti/Tv within difficult regions confirms elevated false positives there.
Quick counts with query
bcftools view -H input.vcf.gz | wc -l
bcftools view -v snps -H input.vcf.gz | wc -l
bcftools view -v indels -H input.vcf.gz | wc -l
bcftools view -f PASS -H input.vcf.gz | wc -l
bcftools query -f '%QUAL\n' input.vcf.gz | awk '{s+=$1;n++} END{print "mean QUAL:", s/n}'
See examples/vcf_stats.py for a cyvcf2 script computing counts, Ti/Tv, and mean QUAL in one pass; the usage guide covers per-sample genotype distributions and the allele-frequency spectrum.
Quick Reference
| Task | Command |
|---|
| Full stats | bcftools stats input.vcf.gz |
| Per-sample het/hom/missing | bcftools stats -s - input.vcf.gz | grep "^PSC" |
| Ti/Tv ratio | bcftools stats input.vcf.gz | grep "^TSTV" | cut -f5 |
| Per-sample missingness | vcftools --gzvcf in.vcf.gz --missing-indv |
| Exact HWE (controls) | vcftools --gzvcf controls.vcf.gz --hardy |
| KING kinship | vcftools --gzvcf in.vcf.gz --relatedness2 |
| Sex/ancestry/relatedness | python -m peddy -p 4 --plot --prefix qc in.vcf.gz in.ped |
| Scalable identity QC | somalier extract ... then somalier relate --ped ... |
| Plot stats | plot-vcfstats -p dir stats.txt |
Common Errors
| Symptom | Cause | Fix |
|---|
| Every AFR sample flagged as high het/hom | Single global cutoff across ancestries | Infer ancestry, flag outliers within each group |
| Good variants filtered by HWE | Two-sided HWE on a mixed cohort or on cases | Excess-het only, within ancestry, controls-only |
| Spurious HWE/missingness | Cohort metrics computed before genotype filtering | Set ./. on low GQ/DP first, then recompute |
bcftools gtcheck -G 1 errors | -G removed after the 1.10 rewrite | Drop -G; cross-check runs by default without -g |
| Low Ti/Tv only in novel set | False positives concentrate in novel variants | Raise stringency; recheck against dbSNP overlap |
plot-vcfstats not found / no plots | matplotlib missing or not on PATH | pip install matplotlib; check which plot-vcfstats |
Related Skills
- variant-calling/filtering-best-practices - Apply the site and genotype filters these metrics motivate
- variant-calling/gatk-variant-calling - VerifyBamID2 contamination estimate feeding
--contamination
- variant-calling/variant-normalization - Normalize before comparing or annotating call sets
- variant-calling/vcf-basics - View, query, and understand VCF fields
- variant-calling/vcf-manipulation - Compare and merge call sets
- variant-calling/joint-calling - Cohort genotyping where population QC applies
- alignment-files/bam-statistics - Upstream alignment QC that drives variant statistics
References
- Danecek P, Bonfield JK, Liddle J, et al. Twelve years of SAMtools and BCFtools. 2021 GigaScience 10:giab008. (bcftools stats/gtcheck)
- Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. 2011 Bioinformatics 27:2156-2158. (vcftools QC)
- Wigginton JE, Cutler DJ, Abecasis GR. A note on exact tests of Hardy-Weinberg equilibrium. 2005 American Journal of Human Genetics 76:887-893. (exact HWE test)
- Manichaikul A, Mychaleckyj JC, Rich SS, et al. Robust relationship inference in genome-wide association studies. 2010 Bioinformatics 26:2867-2873. (KING kinship)
- Pedersen BS, Quinlan AR. Who's Who? Detecting and Resolving Sample Anomalies in Human DNA Sequencing Studies with Peddy. 2017 American Journal of Human Genetics 100:406-413. (peddy)
- Pedersen BS, Bhetariya PJ, Brown J, et al. Somalier: rapid relatedness estimation for cancer and germline studies using efficient genome sketches. 2020 Genome Medicine 12:62. (somalier)
- Zhang F, Flickinger M, Gagliano Taliun SA, et al. Ancestry-agnostic estimation of DNA sample contamination from sequence reads. 2020 Genome Research 30:185-194. (VerifyBamID2)
