Genome analysis scales everything you have learned to billions of bases. This chapter covers variant calling, mutation classification, and the computational challenges of working with whole genomes.

You will parse VCF files (the standard variant format), classify mutations by their effect on protein coding, and understand the workflow from raw sequencing reads to variant calls.

Genome analysis is where bioinformatics meets public health. Every genetic test, every ancestry report, and every pharmacogenomics decision starts with variant analysis.

VCF Parsing and Variant Annotation#

VCF (Variant Call Format) is the standard format for genetic variants. This script parses VCF files and extracts variant information for downstream analysis.

Vanilla version (ch09_vanilla_01.py):

#!/usr/bin/env python3


# Sample VCF data (simplified)
vcf_data = """\
##fileformat=VCFv4.2
##source=sample
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr1	11856378	rs123456	G	A	99	PASS	DP=50;AF=0.45
chr1	11856424	rs789012	T	C	85	PASS	DP=42;AF=0.38
chr7	55191822	rs1212127	T	G	95	PASS	DP=67;AF=0.52
chr17	41245466	rs80357906	A	G	99	PASS	DP=55;AF=0.01
chr12	25398284	rs121913529	C	T	99	PASS	DP=48;AF=0.02
"""

# Parse VCF
print("Parsing VCF file:\n")
print(f"{'CHROM':8s} {'POS':>10s} {'ID':>12s} {'REF':>4s} {'ALT':>4s} {'QUAL':>5s} {'AF':>6s}")
print("-" * 55)

variants = []
for line in vcf_data.strip().split("\n"):
    if line.startswith("#"):
        continue
    fields = line.split("\t")
    chrom, pos, vid, ref, alt, qual, filt, info = fields
    # Extract allele frequency
    af = 0.0
    for item in info.split(";"):
        if item.startswith("AF="):
            af = float(item.split("=")[1])
    variants.append({"chrom": chrom, "pos": int(pos), "id": vid, "ref": ref, "alt": alt, "qual": int(qual), "af": af})
    print(f"{chrom:8s} {pos:>10s} {vid:>12s} {ref:>4s} {alt:>4s} {qual:>5s} {af:>6.2f}")

# Classify variants
print(f"\nTotal variants: {len(variants)}")

snps = [v for v in variants if len(v["ref"]) == 1 and len(v["alt"]) == 1]
indels = [v for v in variants if len(v["ref"]) != len(v["alt"])]
print(f"SNPs: {len(snps)}")
print(f"Indels: {len(indels)}")

# Common vs rare
common = [v for v in variants if v["af"] > 0.05]
rare = [v for v in variants if v["af"] <= 0.05]
print(f"\nCommon (AF > 5%): {len(common)}")
for v in common:
    print(f"  {v['id']}: {v['ref']}>{v['alt']} AF={v['af']:.2f}")
print(f"Rare (AF <= 5%): {len(rare)}")
for v in rare:
    print(f"  {v['id']}: {v['ref']}>{v['alt']} AF={v['af']:.2f}")

We parsed VCF columns (chromosome, position, reference allele, alternate allele, quality), computed simple statistics, and classified variants by type (SNP, indel).

AI version (ch09_ai_01.py):

#!/usr/bin/env python3


import os

try:
    from openai import OpenAI
    client = OpenAI(api_key=os.environ.get("OPENAI_API_KEY", ""), base_url=os.environ.get("OPENAI_BASE_URL", "https://openrouter.ai/api/v1"))
    AI_AVAILABLE = True
except ImportError:
    AI_AVAILABLE = False
    print("Note: Install openai package for AI features\n")

def ask_ai(prompt: str) -> str:
    if not AI_AVAILABLE:
        return "(AI not available)"
    return client.chat.completions.create(model=os.environ.get("OPENAI_MODEL", "google/gemini-2.5-flash"), messages=[{"role": "user", "content": prompt}], temperature=0.3).choices[0].message.content

variants = [
    {"chrom": "chr7", "pos": 55191822, "id": "rs1212127", "ref": "T", "alt": "G", "af": 0.52, "gene": "BRAF"},
    {"chrom": "chr17", "pos": 41245466, "id": "rs80357906", "ref": "A", "alt": "G", "af": 0.01, "gene": "BRCA1"},
    {"chrom": "chr12", "pos": 25398284, "id": "rs121913529", "ref": "C", "alt": "T", "af": 0.02, "gene": "KRAS"},
]

print("Notable variants found:")
for v in variants:
    print(f"  {v['id']}: {v['ref']}>{v['alt']} in {v['gene']} (AF={v['af']:.2f})")

print("\n--- AI: What do these variants mean for health? ---\n")
var_text = "\n".join(f"  {v['id']}: {v['ref']}>{v['alt']} in {v['gene']}, allele frequency={v['af']}" for v in variants)
result = ask_ai(
    f"I found these genetic variants in a patient's genome:\n{var_text}\n\n"
    "Please explain:\n"
    "1. What does each variant do? Are they pathogenic or benign?\n"
    "2. What is the clinical significance of BRCA1 variants?\n"
    "3. What is pharmacogenomics? Could these variants affect drug response?\n"
    "4. How do genetic counselors use this information?\n\n"
    "Be sensitive -- this could be someone's real genetic data."
)
print(result)

The AI interprets variants: which ones are likely pathogenic, what the quality scores mean, and how to prioritize variants for clinical follow-up.

You can run any vanilla script without an API key. The AI scripts will print a message and run in offline mode if no key is set.

Mutation Classification#

Mutations affect proteins in different ways: missense mutations change one amino acid, nonsense mutations create premature stops, and synonymous mutations do not change the protein at all.

Vanilla version (ch09_vanilla_02.py):

#!/usr/bin/env python3


CODON_TABLE = {
    "TTT": "F", "TTC": "F", "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L", "CTA": "L", "CTG": "L",
    "ATT": "I", "ATC": "I", "ATA": "I", "ATG": "M", "GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V",
    "TCT": "S", "TCC": "S", "TCA": "S", "TCG": "S", "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P",
    "ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T", "GCT": "A", "GCC": "A", "GCA": "A", "GCG": "A",
    "TAT": "Y", "TAC": "Y", "TAA": "*", "TAG": "*", "CAT": "H", "CAC": "H", "CAA": "Q", "CAG": "Q",
    "AAT": "N", "AAC": "N", "AAA": "K", "AAG": "K", "GAT": "D", "GAC": "D", "GAA": "E", "GAG": "E",
    "TGT": "C", "TGC": "C", "TGA": "*", "TGG": "W", "CGT": "R", "CGC": "R", "CGA": "R", "CGG": "R",
    "AGT": "S", "AGC": "S", "AGA": "R", "AGG": "R", "GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G",
}


def annotate_variant(codon: str, position: int, new_base: str) -> str:
    """Predict the effect of a single-base change in a codon."""
    ref_aa = CODON_TABLE.get(codon, "?")
    mutated_codon = codon[:position] + new_base + codon[position + 1:]
    alt_aa = CODON_TABLE.get(mutated_codon, "?")

    if ref_aa == alt_aa:
        return "synonymous"
    elif alt_aa == "*":
        return "nonsense (stop gain)"
    elif ref_aa == "*":
        return "stop loss"
    else:
        return "missense"


# Known pathogenic variants
variants = [
    {"gene": "BRAF", "codon": "GTG", "pos": 0, "alt": "T", "note": "V600E -- most common BRAF mutation in cancer"},
    {"gene": "KRAS", "codon": "GGT", "pos": 1, "alt": "A", "note": "G12D -- common in pancreatic cancer"},
    {"gene": "HBB", "codon": "GAG", "pos": 1, "alt": "T", "note": "E6V -- causes sickle cell disease"},
]

print("Variant annotation:")
print("=" * 60)
for v in variants:
    ref_aa = CODON_TABLE.get(v["codon"], "?")
    effect = annotate_variant(v["codon"], v["pos"], v["alt"])
    mut_codon = v["codon"][:v["pos"]] + v["alt"] + v["codon"][v["pos"] + 1:]
    alt_aa = CODON_TABLE.get(mut_codon, "?")
    print(f"\n  {v['gene']}: {v['codon']}({ref_aa}) -> {mut_codon}({alt_aa})")
    print(f"  Effect: {effect}")
    print(f"  Note: {v['note']}")

We built a codon table and classified mutations as missense, nonsense, or synonymous based on how they change the protein sequence.

AI version (ch09_ai_02.py):

#!/usr/bin/env python3


import os

try:
    from openai import OpenAI
    client = OpenAI(api_key=os.environ.get("OPENAI_API_KEY", ""), base_url=os.environ.get("OPENAI_BASE_URL", "https://openrouter.ai/api/v1"))
    AI_AVAILABLE = True
except ImportError:
    AI_AVAILABLE = False
    print("Note: Install openai package for AI features\n")

def ask_ai(prompt: str) -> str:
    if not AI_AVAILABLE:
        return "(AI not available)"
    return client.chat.completions.create(model=os.environ.get("OPENAI_MODEL", "google/gemini-2.5-flash"), messages=[{"role": "user", "content": prompt}], temperature=0.3).choices[0].message.content

variants = [
    {"gene": "BRAF", "change": "V600E", "effect": "missense", "disease": "melanoma"},
    {"gene": "KRAS", "change": "G12D", "effect": "missense", "disease": "pancreatic cancer"},
    {"gene": "HBB", "change": "E6V", "effect": "missense", "disease": "sickle cell disease"},
]

print("Pathogenic variants:")
for v in variants:
    print(f"  {v['gene']} {v['change']}: {v['effect']} -- {v['disease']}")

print("\n--- AI: From mutation to disease mechanism ---\n")
result = ask_ai(
    f"These three mutations all change one amino acid:\n"
    f"  BRAF V600E: valine -> glutamic acid at position 600\n"
    f"  KRAS G12D: glycine -> aspartic acid at position 12\n"
    f"  HBB E6V: glutamic acid -> valine at position 6\n\n"
    "For each mutation:\n"
    "1. How does ONE amino acid change cause disease?\n"
    "2. What is the molecular mechanism?\n"
    "3. Is there a drug that targets this specific mutation?\n"
    "4. How common is this mutation in the population?\n\n"
    "Then explain: why are some single-letter changes devastating "
    "while others have no effect? Use an analogy."
)
print(result)

The AI explains the clinical significance of each mutation type: why nonsense mutations are usually more harmful, why some missense mutations are benign, and how to use mutation classification in genetic counseling.

The AI version always includes the same core logic as the vanilla version. The AI calls are added at the end, so you can compare the two approaches side by side.

Chapter Summary#

This chapter covered genome-scale variant analysis. VCF parsing and mutation classification.

The vanilla scripts gave you hands-on experience with the core techniques. The AI scripts showed how LLMs can interpret, explain, and extend your analysis. Together, they prepare you for real-world bioinformatics work where code and AI work side by side.