domains / data / vcf-quality-control / vcf-statistics

vcf-statistics

Skill card

combined

The two skills cover the same task (VCF statistics and QC) but with different emphases and non-overlapping material. The OpenClaw copy is a broad command reference (bcftools stats, gtcheck, query, cyvcf2) with specific expected Ti/Tv ranges and plotting steps. The GPTomics copy is a decision-oriented QC guide that adds a governing principle (cohort-relative interpretation, order of operations, identity as a graph problem), a detailed decision table for interpreting deviations, version compatibility notes, and additional tools (vcftools, somalier, peddy, KING). Neither is a strict reproduction of the other; crowning one would hide the other's unique interpretive framework and tooling.

Repositories

Other copies

Combined skill

# VCF Statistics and Quality Control

## 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

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.

## 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. QC is not threshold-checking but reading a deviation for its mechanistic meaning and acting on it.

1. Always compare a sample against a matched cohort (same assay, same inferred ancestry), never against an absolute number.
2. Order of operations matters: apply genotype-level filters (set `./.` where GQ/DP/allele-balance fail) BEFORE computing cohort missingness, call rate, or HWE.
3. Sample-identity QC is a graph problem: build the all-pairs relatedness matrix and reconcile it against the declared pedigree/manifest BEFORE any association or burden analysis.

## QC Metric Decision Table

| Metric | Expected (WGS) | Expected (WES) | A deviation MEANS | Action |
|--------|---------------|---------------|-------------------|--------|
| Ti/Tv (overall) | ~2.0-2.1 | ~3.0-3.3 | Low -> false-positive transversions; high -> over-filtering | 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 known-site value -> FP contamination | Raise stringency |
| Het/hom ratio | ancestry-dependent | same | High vs same-ancestry cohort -> contamination; low -> inbreeding/ROH | Stratify by ancestry first |
| Novel fraction (vs dbSNP) | low for common | same | Novel COMMON variants -> almost always artifacts | Investigate common novels |
| Call rate (per sample) | >95-98% | >95-98% | Low -> low-coverage/low-quality sample | Drop worst samples, then recompute |
| Call rate (per variant) | >95-99% | >95-99% | Low -> hard-to-genotype region | Drop worst variants after sample QC |
| Het allele-balance | centered on 0.5 | 0.5 | Shifted away from 0.5 -> contamination signature | Confirm with VerifyBamID2 or CHARR |
| Excess-het / HWE | in equilibrium | same | EXCESS het -> mapping artifact; het DEFICIT -> often real | Filter on excess het only, within ancestry |
| Relatedness (kinship) | matches manifest | matches manifest | Unexpected high kinship -> sample swap/duplicate | Reconcile with KING/somalier/peddy |

## bcftools stats

### Basic Statistics
```bash
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
```

### View Key Metrics
```bash
bcftools stats input.vcf.gz | grep "^SN"
bcftools stats input.vcf.gz | grep "^TSTV"
```

Output sections:
- `SN` - Summary numbers
- `TSTV` - Transitions/transversions
- `SiS` - Singleton stats
- `AF` - Allele frequency distribution
- `QUAL` - Quality distribution
- `IDD` - Indel distribution
- `ST` - Substitution types
- `DP` - Depth distribution

### Region-Specific Stats
```bash
bcftools stats -r chr1:1000000-2000000 input.vcf.gz > region_stats.txt
bcftools stats -R exome.bed input.vcf.gz > exome_stats.txt
```

### Plotting Statistics
```bash
bcftools stats input.vcf.gz > stats.txt
plot-vcfstats -p output_dir stats.txt
```

## bcftools gtcheck

### Check Sample Identity
```bash
bcftools gtcheck -g reference.vcf.gz query.vcf.gz
```

### Detect Sample Swaps
```bash
bcftools gtcheck input.vcf.gz > relatedness.txt
```

### Check Against Reference Panel
```bash
bcftools gtcheck -g 1000genomes.vcf.gz unknown_sample.vcf.gz
```

## Quick Statistics with Query

### Count Variants
```bash
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
```

### Quality and Depth Distribution
```bash
bcftools query -f '%QUAL\n' input.vcf.gz | awk '{sum+=$1; count++} END {print "Mean QUAL:", sum/count}'
bcftools query -f '%INFO/DP\n' input.vcf.gz | awk '{sum+=$1; count++} END {print "Mean DP:", sum/count}'
```

### Allele Frequency Spectrum
```bash
bcftools query -f '%INFO/AF\n' input.vcf.gz | awk '{
    if($1<0.01) rare++
    else if($1<0.05) low++
    else if($1<0.5) common++
    else freq++
} END {
    print "Rare (<1%):", rare
    print "Low (1-5%):", low
    print "Common (5-50%):", common
    print "Frequent (>50%):", freq
}'
```

## Sample Statistics

### List and Count Samples
```bash
bcftools query -l input.vcf.gz
bcftools query -l input.vcf.gz | wc -l
```

### Per-Sample Variant Counts
```bash
for sample in $(bcftools query -l input.vcf.gz); do
    count=$(bcftools view -s "$sample" -H input.vcf.gz | bcftools view -c 1 -H | wc -l)
    echo "$sample: $count"
done
```

### Missing Genotypes per Sample
```bash
bcftools stats -s - input.vcf.gz | grep "^PSC"
```

## cyvcf2 Statistics

### Basic Counts
```python
from cyvcf2 import VCF

stats = {'snps': 0, 'indels': 0, 'other': 0}

for variant in VCF('input.vcf.gz'):
    if variant.is_snp:
        stats['snps'] += 1
    elif variant.is_indel:
        stats['indels'] += 1
    else:
        stats['other'] += 1

print(f'SNPs: {stats["snps"]}')
print(f'Indels: {stats["indels"]}')
print(f'Other: {stats["other"]}')
```

### Quality Statistics
```python
from cyvcf2 import VCF
import numpy as np

quals = []
for variant in VCF('input.vcf.gz'):
    if variant.QUAL:
        quals.append(variant.QUAL)

quals = np.array(quals)
print(f'Mean QUAL: {np.mean(quals):.1f}')
print(f'Median QUAL: {np.median(quals):.1f}')
print(f'Min QUAL: {np.min(quals):.1f}')
print(f'Max QUAL: {np.max(quals):.1f}')
```