# Bio Applied Variant Calling And Snp Analysis > Run GATK/bcftools BAM-to-VCF calling, parse VCF fields, decode genotypes (GT/AD/DP/GQ), hard-filter variants, test Hardy-Weinberg equilibrium. Use for SNP/indel calling, VCF/GVCF parsing, zygosity decoding, or HWE checks. - Skill: `pavel-kravchenko/bio-applied-variant-calling-and-snp-analysis` (Agent Skill) - Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-variant-calling-and-snp-analysis` - Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-variant-calling-and-snp-analysis/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/pavel-kravchenko/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). ```bash # 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 ```text ##fileformat=VCFv4.2 ##INFO= ##FORMAT= #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. ```python 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). ```python 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. ```python 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