Bio applied variant calling and snp analysis

Variant Calling and SNP Analysis

When to Use

• Building a BAM→VCF pipeline with GATK HaplotypeCaller or bcftools and need read-group/BQSR/joint-genotyping steps in the right order
• Parsing a VCF file's INFO/FORMAT columns into structured records without a full library
• Decoding a genotype string (0/1, 1|2, ./.) into alleles, zygosity, and phasing
• Applying quality filters (QUAL, DP, GQ, FILTER) or checking allele balance to flag suspect heterozygous calls
• Testing whether a SNP's genotype counts deviate from Hardy-Weinberg equilibrium (population stratification, genotyping error, selection)

Version Compatibility

GATK ≥4.5, bcftools ≥1.19, bwa ≥0.7.17, samtools ≥1.19, Python ≥3.10 with numpy ≥1.26 and scipy ≥1.11.

Prerequisites

• gatk4, bcftools, bwa, samtools on PATH (conda: bioconda::gatk4 bioconda::bcftools bioconda::bwa bioconda::samtools)
• A reference FASTA indexed with samtools faidx and gatk CreateSequenceDictionary
• Familiarity with SAM/BAM basics and the VCF spec (##INFO/##FORMAT header lines)

Variant Type Reference

Type	Size	Example
SNV	1 bp	A→G
MNV	2+ bp equal length	AT→GC
Insertion	1–50 bp	A→ATCG
Deletion	1–50 bp	ATCG→A
Large SV	>50 bp	detected by split reads / depth
CNV	variable	depth-based (deletion=low depth, dup=high depth)

Goal: Go from aligned reads to a filtered, joint-genotyped VCF. Approach: Follow GATK Best Practices — mark duplicates, recalibrate base quality, call per-sample GVCFs, joint-genotype, then hard-filter (or VQSR for large cohorts).

# BWA alignment — always set read groups or MarkDuplicates/HaplotypeCaller will fail
bwa mem -R '@RG\tID:s1\tSM:s1\tPL:ILLUMINA\tLB:lib1' ref.fa R1.fq R2.fq | samtools sort -o sorted.bam
samtools index sorted.bam

# Mark duplicates
gatk MarkDuplicates -I sorted.bam -O dedup.bam -M metrics.txt

# BQSR (optional but recommended when known-sites VCF is available)
gatk BaseRecalibrator -I dedup.bam -R ref.fa --known-sites dbsnp.vcf -O recal.table
gatk ApplyBQSR -I dedup.bam -R ref.fa --bqsr-recal-file recal.table -O recal.bam

# Call variants per-sample (GVCF mode enables later cohort joint-genotyping)
gatk HaplotypeCaller -I recal.bam -R ref.fa -O raw.g.vcf.gz -ERC GVCF

# Joint genotyping across samples
gatk GenomicsDBImport --genomicsdb-workspace-path gdb -V raw.g.vcf.gz -L intervals.list
gatk GenotypeGVCFs -R ref.fa -V gendb://gdb -O genotyped.vcf.gz

# Hard filters (use VQSR instead for >=30 WGS or >=10 WES samples)
gatk VariantFiltration -R ref.fa -V genotyped.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 60.0" --filter-name "StrandBias" \
  --filter-expression "MQ < 40.0" --filter-name "LowMQ" \
  -O filtered.vcf.gz

# Fast single-sample alternative with bcftools
bcftools mpileup -f ref.fa recal.bam | bcftools call -mv -Oz -o calls.vcf.gz

Variant Calling Tool Comparison

Tool	Algorithm	Best For
GATK HaplotypeCaller	Local de novo assembly	WGS/WES cohorts, gold standard
bcftools mpileup/call	Pileup-based	Fast single-sample WGS
FreeBayes	Bayesian haplotype	Low-frequency variants
DeepVariant	Deep learning (CNN)	High accuracy WGS/WES
Strelka2	Statistical model	Tumor-normal somatic

VCF Format and Genotype Encodings

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Read Depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS    ID       REF  ALT  QUAL  FILTER  INFO          FORMAT  SAMPLE1
chr1    10000  rs123456 A    G    5000  PASS    DP=200;AF=0.5 GT:DP   0/1:100

GT	Meaning
0/0	Hom-ref
0/1	Het
1/1	Hom-alt
1/2	Het multi-allelic
./.	Missing
0|1	Phased het

Goal: Parse a VCF record and decode a sample's genotype into alleles, zygosity, and phasing. Approach: Split INFO on ;/= into a dict, split each sample's FORMAT-keyed values on :, then map GT allele indices back onto REF/ALT.

def classify_variant(ref: str, alt: str) -> str:
    """Classify a variant by comparing REF and ALT allele lengths."""
    if len(ref) == 1 and len(alt) == 1:
        return 'SNV'
    elif len(ref) > 1 and len(alt) > 1 and len(ref) == len(alt):
        return 'MNV'
    elif len(ref) < len(alt):
        return 'insertion'
    elif len(ref) > len(alt):
        return 'deletion'
    return 'complex'


def parse_vcf_record(line: str, header: list[str]) -> dict:
    """Parse one VCF data line into a structured record with decoded INFO and samples."""
    fields = line.rstrip('\n').split('\t')
    info = {}
    for item in fields[7].split(';'):
        if '=' in item:
            k, v = item.split('=', 1)
            info[k] = v
        else:
            info[item] = True  # flag field, e.g. "DB"

    record = {
        'chrom': fields[0], 'pos': int(fields[1]), 'id': fields[2],
        'ref': fields[3], 'alt': fields[4].split(','),
        'qual': float(fields[5]) if fields[5] != '.' else None,
        'filter': fields[6], 'info': info,
        'type': classify_variant(fields[3], fields[4].split(',')[0]),
    }

    if len(fields) > 9:
        fmt_keys = fields[8].split(':')
        record['samples'] = {}
        for i, sample_data in enumerate(fields[9:]):
            values = sample_data.split(':')
            record['samples'][header[9 + i]] = dict(zip(fmt_keys, values))
    return record


def decode_genotype(gt_string: str, ref: str, alts: list[str]) -> tuple[str, str, bool]:
    """Decode a GT string (e.g. '0/1', '1|2', './.') into (alleles, zygosity, is_phased)."""
    alleles = [ref] + alts
    sep = '|' if '|' in gt_string else '/'
    indices = gt_string.split(sep)
    decoded = [alleles[int(i)] if i != '.' else '.' for i in indices]

    if '.' in indices:
        zygosity = 'missing'
    elif len(set(indices)) == 1:
        zygosity = 'hom-ref' if indices[0] == '0' else 'hom-alt'
    else:
        zygosity = 'het'
    return '/'.join(decoded), zygosity, sep == '|'

Goal: Filter variants on QUAL/DP/GQ and flag suspiciously skewed heterozygous calls. Approach: Reject on any failing criterion, tracking the failure reason; compute allele balance from AD for heterozygotes (expect ~0.5).

from collections import defaultdict


def filter_variants(variants: list[dict], min_qual=30, min_dp=10, require_pass=False):
    """Apply QUAL/DP/FILTER thresholds; returns (passed, {reason: [failed_variants]})."""
    passed, failed = [], defaultdict(list)
    for v in variants:
        reasons = []
        if v['qual'] is not None and v['qual'] < min_qual:
            reasons.append(f"LowQUAL({v['qual']}<{min_qual})")
        dp = int(v['info'].get('DP', 0))
        if dp < min_dp:
            reasons.append(f"LowDP({dp}<{min_dp})")
        if require_pass and v['filter'] not in ('PASS', '.'):
            reasons.append(f"NotPASS({v['filter']})")
        if reasons:
            for r in reasons:
                failed[r].append(v)
        else:
            passed.append(v)
    return passed, failed


def allele_balance(ad_string: str) -> float | None:
    """Compute alt/(ref+alt) read-depth ratio from an AD field; ~0.5 expected for true hets."""
    counts = [int(x) for x in ad_string.split(',') if x != '.']
    if len(counts) < 2 or sum(counts) == 0:
        return None
    return counts[1] / sum(counts)

Goal: Test whether a SNP's genotype counts are consistent with Hardy-Weinberg equilibrium. Approach: Compute allele frequencies (p, q), derive expected genotype counts (p², 2pq, q²), and run a 1-df chi-squared test.

import numpy as np
from scipy.stats import chi2


def hardy_weinberg_test(genotypes: list[tuple]) -> dict:
    """Chi-squared HWE test. genotypes: list of (0,0)/(0,1)/(1,1) tuples; None entries are skipped."""
    valid = [g for g in genotypes if None not in g]
    n = len(valid)
    obs_hom_ref = sum(1 for g in valid if g == (0, 0))
    obs_het = sum(1 for g in valid if g in ((0, 1), (1, 0)))
    obs_hom_alt = sum(1 for g in valid if g == (1, 1))

    alt_count = sum(sum(g) for g in valid)
    q = alt_count / (2 * n)  # alt allele frequency
    p = 1 - q

    expected = [p**2 * n, 2 * p * q * n, q**2 * n]
    observed = [obs_hom_ref, obs_het, obs_hom_alt]
    chi2_stat = sum((o - e) ** 2 / e for o, e in zip(observed, expected) if e > 0)
    p_value = 1 - chi2.cdf(chi2_stat, df=1)

    return {'n': n, 'p': p, 'q': q, 'chi2': chi2_stat, 'p_value': p_value,
            'in_hwe': p_value > 0.05}

Pitfalls

• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based inclusive — off-by-one errors happen when mixing them
• Missing read groups: GATK MarkDuplicates and HaplotypeCaller require @RG tags — always pass -R to bwa mem
• VQSR vs hard filters: VQSR requires ≥30 WGS samples or ≥10 WES samples; use hard filters for smaller cohorts
• Multi-allelic sites: ALT may be comma-separated; split (or run bcftools norm -m-) before per-allele analysis
• Allele balance outliers: hets with AB far from 0.5 (e.g. <0.2 or >0.8) often indicate mapping artifacts or contamination
• HWE deviation isn't always error: population stratification, selection, or a genuine batch effect can also violate HWE — check before discarding a SNP
• Multiple testing: Apply Benjamini-Hochberg FDR when testing thousands of variants (e.g. genome-wide HWE or GWAS scans)

See Also

• bio-applied-snp-calling-pipeline — end-to-end SNP calling pipeline walkthrough
• bio-applied-gwas — genome-wide association testing on called variants
• bio-applied-population-genetics — allele frequencies and population structure
• bio-applied-clinical-genomics — clinical annotation and interpretation of called variants