Mutations and genetic variation

Why this matters for bioinformatics

No two human genomes are the same. The 1000 Genomes Project found that a typical person’s genome differs from the reference genome at 4.1 million to 5.0 million sites. Almost all of these differences are harmless. A handful can cause disease, and in cancer, mutations that arise in one tissue during a person’s life drive the growth of a tumour.

Finding variants, describing them unambiguously and working out what they do are among the most common jobs in bioinformatics. A clinical geneticist looks for the one variant behind a child’s rare disease, a cancer team catalogues the mutations in a tumour, a population geneticist reads human history from shared variants, and a public-health lab tracks new virus lineages from their mutations. All of them use the same file formats, naming rules and reasoning, which this topic introduces.

Where variants come from

A mutation is a change in a DNA sequence; the word describes the event. A variant is any position where a genome differs from a reference, and clinical genetics prefers this neutral word because it makes no claim about harm. The different versions found at one position are its alleles. Two more distinctions matter in every analysis:

• Germline variants are present in an egg or a sperm cell, so they are in every cell of the child and can be passed on. Variants inherited from parents and de novo mutations, which appear for the first time in a child, are both germline.
• Somatic variants arise in one cell of the body after conception and are found only in that cell’s descendants. They are not inherited, and they are what tumour sequencing looks for.

Copying three billion base pairs is remarkably accurate, but not perfect. Sequencing parents and children shows about 1.2 new mutations per 100 million bases per generation, which works out at roughly 70 de novo mutations in each child. Most come from the father, and each extra year of the father’s age adds about one to two.

How mutations happen

• Copying errors. DNA polymerase occasionally inserts the wrong base. Proofreading and mismatch repair catch most of these mistakes, but not all.
• Chemical damage. Bases are damaged by water, oxygen, ultraviolet light and chemicals such as those in tobacco smoke. Each source leaves its own pattern of changes, a mutational signature, that can be read in tumour genomes: ultraviolet light mainly causes C→T changes next to another pyrimidine.
• Slippage in repeats. In runs such as AAAAAA or CACACA, the two strands can slip out of register during copying, adding or removing repeat units. That is why small insertions and deletions cluster in repeats.
• Recombination between similar sequences. When two nearby copies of a sequence pair up by mistake, a crossover between them deletes or duplicates the DNA between them. The β-globin cluster is a classic case: haemoglobin Lepore is made by a gene fused from part of HBD and part of HBB, produced by exactly this kind of crossover.

The most common mutation: C to T at CpG

Cytosine can lose its amino group, a reaction called deamination. Ordinary cytosine becomes uracil, a base that does not belong in DNA, and the cell’s repair enzymes remove it. But in human DNA most cytosines that are followed by a guanine (a CpG site) carry a methyl group, and 5-methylcytosine deaminates to thymine: a normal DNA base, so the mistake is easy to miss. If it is not repaired before the DNA is copied, the C·G pair becomes a T·A pair.

As a result, C→T changes at CpG sites happen roughly ten times more often than changes elsewhere, and over millions of years they have worn CpG down: in the human genome the dinucleotide CG is found at about a fifth of the frequency its C and G content would predict. You will measure both effects in the code lab.

Types of variant: from one base to whole chromosomes

Variants are classified first by size, because size decides how they are detected, stored and interpreted.

Counts per genome from the 1000 Genomes Project, phase 3 (median autosomal sites per genome across five continental groups). People with more African ancestry carry the most variants, because human genetic diversity is greatest in Africa.

More than 99.9% of the variants in a genome are SNVs and short indels, yet structural variants change more bases in total. Copy-number variants (CNVs) are structural variants that change how many copies of a stretch of DNA a person has. Variants of different sizes need different methods: SNVs and short indels are called from short reads lined up against the reference, while large and complex structural variants are found much more reliably with long reads.

Two kinds of SNV: transitions and transversions

The four DNA bases come in two shapes: the purines A and G have two rings, the pyrimidines C and T have one. A transition swaps one purine for the other (A↔G) or one pyrimidine for the other (C↔T). A transversion swaps a purine for a pyrimidine or the reverse. There are twice as many possible transversions as transitions, yet transitions are more common, partly because of CpG deamination. In high-quality whole-genome data the transition-to-transversion ratio (Ti/Tv) is about 2.0 to 2.1. Random sequencing errors would give 0.5, so a low Ti/Tv is a classic sign that a set of variant calls contains many false ones.

Describing a variant: VCF and HGVS

A variant needs a name that points to exactly one change in exactly one reference sequence. Two systems do most of the work: the VCF file format, used by software, and HGVS nomenclature, used in papers and clinical reports.

VCF: one variant per line

Here is the sickle-cell variant as it appears in ClinVar’s VCF file, with the long INFO column shortened:

11  5227002  15333  T  A  .  .  ALLELEID=30372;CLNHGVS=NC_000011.10:g.5227002T>A;CLNSIG=Pathogenic;MC=SO:0001583|missense_variant;RS=334

Indels are written with the base before them included, the anchor base, so that REF and ALT are never empty. The four-base deletion c.126_129del is 11 5226762 CAAAG C: the C at 5,226,762 is kept and the four bases after it are deleted. A position with several alternative alleles can be written on one line (ALT = A,G) or split into one line per allele, as both files in the code lab do.

HGVS: names that say where and what

HGVS names start with a reference sequence and a prefix that says what kind of sequence it is:

• g. a genomic position: NC_000011.10:g.5227002T>A
• c. a position in the coding DNA of a transcript, where the A of the start codon is c.1: NM_000518.5:c.20A>T. Bases before the start codon get negative numbers (c.−78) and bases after the stop codon a star (c.*110). Intronic bases are numbered from the nearest exon: c.92+1 is the first base of intron 1, after exon base c.92, and c.93−21 is 21 bases before exon base c.93.
• p. the predicted change in the protein: NP_000509.1:p.(Glu7Val), with the parentheses saying it is predicted rather than measured. Stops are Ter (or *), so a nonsense variant is p.Gln40Ter, and a frameshift is p.Phe42LeufsTer19: phenylalanine 42 becomes leucine, and the new reading frame ends in a stop at position 19, counting the leucine as position 1.

Changes are written as A>T (substitution), del (deletion), dup (duplication), ins (insertion) and delins (deletion-insertion).

Genomic names are on NC_000011.10 (GRCh38); coding names on NM_000518.5. Each was checked against the ClinVar records in the code lab.

Notice that the bases in the g. and c. names are complements of each other. HBB lies on the minus strand, so the transcript reads the opposite strand to the one stored in the genome FASTA, and c. numbers run downwards along the chromosome. This is the same strand logic as in the previous topic, and it trips up people all the time.

What a variant does to a gene

Whether a variant matters depends above all on where it falls and what it changes. In a protein-coding gene:

• Synonymous variants change a codon but not its amino acid, such as CAT to CAC, both histidine. They are usually harmless, though a few disrupt splicing.
• Missense variants change one amino acid. The effect ranges from none to severe, depending on how different the new amino acid is and where it sits in the folded protein. Glu→Val at β-globin residue 6 replaces a charged surface residue with a hydrophobic one, as you saw in the previous topic.
• Nonsense variants turn a codon into a stop codon. Frameshifts, indels whose length is not a multiple of three, change every codon after them and usually meet a stop codon soon afterwards. Indels of three bases, or multiples of three, remove or add whole amino acids instead (in-frame indels).
• Start-lost and stop-lost variants destroy the start or stop codon. Translation then starts somewhere else, or runs on into the 3′ UTR.
• Splice-site variants change the GT at the start of an intron or the AG at its end, or other bases the spliceosome recognises, so exons are skipped or parts of introns are kept.
• Regulatory variants in promoters, enhancers and untranslated regions change how much protein is made rather than what it is like.

β-thalassaemia: one gene, every mechanism

People with β-thalassaemia make too little β-globin (β⁺) or none at all (β⁰). Almost 300 different HBB mutations are known to cause it, and between them they break almost every step in making a protein. That makes HBB the best single gene for learning how variants act:

• Transcription. c.−78A>G changes the TATA box, 28 bases upstream of the transcription start, and c.−138C>T a CACCC box. Less mRNA is made.
• Splicing. c.92+1G>A destroys the GT at the start of intron 1, so the normal splice cannot happen. c.93−21G>A (IVS-I-110) creates a new AG 19 bases inside intron 1, and the spliceosome uses it most of the time, keeping 19 intron bases in the mRNA. c.316−197C>T (IVS-II-654) creates a new splice donor deep in intron 2, and 73 intron bases are spliced in as an extra exon.
• Processing the mRNA end. c.*110T>C changes the poly(A) signal AATAAA to AACAAA, so the mRNA is cut and polyadenylated poorly.
• Translation. c.1A>G removes the start codon. c.118C>T (codon 39) makes a premature stop, and c.126_129del a frameshift.
• Both at once. HbE, c.79G>A, is a missense change, Glu→Lys, and also creates the sequence GGTAAG, which looks like a splice donor. Some transcripts are spliced at the wrong place, so less normal β-globin is made.

Nonsense-mediated decay: why many stops make no protein at all

Cells check new mRNAs for premature stop codons. During splicing, a protein marker is left on the mRNA at each exon–exon junction, and the first ribosome to translate the mRNA removes the markers as it passes. If the ribosome stops more than about 50 to 55 bases before the last junction, markers remain downstream, and the mRNA is destroyed by nonsense-mediated decay (NMD). In HBB the last junction is between c.315 and c.316, so the codon 39 stop (c.118–120) triggers NMD: no β-globin is made, carriers are mostly healthy and the disease is recessive. A stop in exon 3, such as c.364G>T, escapes NMD. The cell then makes a shortened, unstable β-globin that damages red blood cells, and one copy is enough to cause disease: dominant β-thalassaemia. The rule has exceptions, such as stops very close to the start codon, but it is the first thing variant-effect software checks.

Variants in populations

For each autosomal position you carry two alleles, one from each parent: your genotype. In a VCF file 0 means the reference allele and 1 the first alternative, so 0/0 is homozygous for the reference, 0/1 heterozygous and 1/1 homozygous for the alternative. A vertical bar, as in 0|1, means the genotype is phased: the file records which allele came on which chromosome copy.

The allele frequency is the fraction of all copies in a population that carry an allele: (number of heterozygotes + 2 × number of homozygotes) / (2 × number of people). If people pair up at random with respect to a variant, genotype frequencies follow from allele frequencies. With alleles of frequency p and q = 1 − p, the expected proportions are p² (two reference copies), 2pq (one of each) and q² (two alternative copies): the Hardy–Weinberg proportions. Departures from them can point to selection, to how samples were chosen, or to genotyping errors.

Rare and common

Most variants are rare, but most of the variants in any one genome are common. A new mutation starts in one person, and only a few reach high frequency. In the 1000 Genomes data around HBB, 72 of the 141 SNVs have a frequency below 0.1%, while 25 are carried by at least 5% of chromosomes. Rare variants also tend to be confined to one part of the world: 86% of the variants found in the 1000 Genomes Project were seen in only one continental group.

Frequency is one of the strongest clues in interpretation. A variant carried by 30% of healthy people cannot cause a rare, severe childhood disease, so clinical pipelines filter variants against population databases such as gnomAD, whose 2020 release alone summarised the variants of 141,456 people. Frequencies must be compared within matching ancestry groups, because a variant common in one population can be absent from another.

The sickle-cell allele and malaria

People with one copy of the sickle-cell allele (sickle-cell trait) are usually healthy and are partly protected against severe malaria, which Anthony Allison reported in 1954. People with two copies have sickle-cell disease. Where malaria has been common, the protection given to carriers has kept the allele at high frequency despite the harm it does to people with two copies, an example of balancing selection. In the code lab you will count the allele in 2,504 unrelated people from 26 populations.

Code lab: real variant data

Each box below is Python that runs in your browser: press Run (or Ctrl + Enter). The first run downloads Python, which takes a few seconds; after that, runs are quick. Your edits are saved in this browser until you press Reset.

You will work with two real VCF files. One holds ClinVar’s interpretations of 2,105 variants in the β-globin region; the other holds the genotypes of all 3,202 people in the high-coverage 1000 Genomes data at 176 variants around HBB. Boxes marked Exercise are checked automatically, with a hint and a worked solution.

import gzip
from collections import Counter

def parse_info(text):
    """Turn 'RS=334;CLNSIG=Pathogenic;...' into a dictionary."""
    info = {}
    for item in text.split(";"):
        key, has_value, value = item.partition("=")
        info[key] = value if has_value else True        # a key without "=" is a flag
    return info

header, records = [], []
with gzip.open("clinvar_20261019_chr11_5220001-5300000.vcf.gz", "rt") as handle:
    for line in handle:
        if line.startswith("#"):
            header.append(line.rstrip("\n"))
            continue
        chrom, pos, vid, ref, alt, qual, filt, info = line.rstrip("\n").split("\t")
        records.append({"chrom": chrom, "pos": int(pos), "id": vid, "ref": ref, "alt": alt, "info": parse_info(info)})

print(len(header), "header lines. The last one names the columns:")
print(header[-1])
print(len(records), "variants, by type:")
for kind, n in Counter(record["info"]["CLNVC"] for record in records).most_common():
    print(f"  {kind:<26}{n:>5}")

sickle = next(r for r in records if r["pos"] == 5227002 and r["alt"] == "A")
print()
print("The sickle-cell variant:", sickle["chrom"], sickle["pos"], sickle["id"], sickle["ref"], sickle["alt"])
for key in ("RS", "CLNHGVS", "CLNSIG", "CLNREVSTAT", "MC"):
    print(f"  {key:<11}{sickle['info'][key]}")

import gzip

def is_transition(ref, alt):
    """True for a change between two purines (A, G) or two pyrimidines (C, T)."""
    return False     # TODO

snvs = []
with gzip.open("1kGP_high_coverage_chr11_5225001-5230000.vcf.gz", "rt") as handle:
    for line in handle:
        if line.startswith("#"):
            continue
        chrom, pos, vid, ref, alt = line.split("\t", 5)[:5]     # stop splitting before the 3,202 genotypes
        if len(ref) == 1 and len(alt) == 1:                     # SNVs only: skip indels and structural variants
            snvs.append((int(pos), ref, alt))

transitions = sum(1 for pos, ref, alt in snvs if is_transition(ref, alt))
transversions = len(snvs) - transitions
print(len(snvs), "SNVs:", transitions, "transitions and", transversions, "transversions")
print(f"Ti/Tv = {transitions / transversions:.2f}")

import gzip
def read_fasta(path):
    """Return (header, sequence) for a single-record FASTA file."""
    with open(path) as handle:
        header = handle.readline().strip()
        return header, "".join(line.strip() for line in handle)

_, region = read_fasta("NC_000011.10_5220001-5300000.fasta")
OFFSET = 5220001                           # the FASTA file starts at chromosome position 5,220,001
start, end = 5225001, 5230000              # the window covered by the 1000 Genomes file

cpg = set()                                # chromosome positions of the C and the G of every CpG
for pos in range(start, end):
    if region[pos - OFFSET:pos - OFFSET + 2] == "CG":
        cpg.update({pos, pos + 1})

snvs = []
with gzip.open("1kGP_high_coverage_chr11_5225001-5230000.vcf.gz", "rt") as handle:
    for line in handle:
        if line.startswith("#"):
            continue
        chrom, pos, vid, ref, alt = line.split("\t", 5)[:5]     # stop splitting before the 3,202 genotypes
        if len(ref) == 1 and len(alt) == 1:                     # SNVs only: skip indels and structural variants
            snvs.append((int(pos), ref, alt))

at_cpg = [snv for snv in snvs if snv[0] in cpg]
other_snvs = len(snvs) - len(at_cpg)
other_bases = (end - start + 1) - len(cpg)
print(f"CpG bases:   {len(at_cpg)} SNVs at {len(cpg)} bases, {len(at_cpg) / len(cpg):.1%} of them")
print(f"Other bases: {other_snvs} SNVs at {other_bases:,} bases, {other_snvs / other_bases:.1%} of them")
print(f"A CpG base is {(len(at_cpg) / len(cpg)) / (other_snvs / other_bases):.1f} times as likely to carry an SNV")
deaminated = [snv for snv in at_cpg if (snv[1], snv[2]) in {("C", "T"), ("G", "A")}]
print(len(deaminated), "of the", len(at_cpg), "CpG variants are C>T, or G>A where the C is on the minus strand")

bases = "TCAG"
amino_acids = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
codon_table = {a + b + c: amino_acids[16 * i + 4 * j + k]
               for i, a in enumerate(bases)
               for j, b in enumerate(bases)
               for k, c in enumerate(bases)}

THREE = {"A": "Ala", "R": "Arg", "N": "Asn", "D": "Asp", "C": "Cys", "Q": "Gln", "E": "Glu", "G": "Gly", "H": "His",
         "I": "Ile", "L": "Leu", "K": "Lys", "M": "Met", "F": "Phe", "P": "Pro", "S": "Ser", "T": "Thr", "W": "Trp",
         "Y": "Tyr", "V": "Val", "*": "Ter"}

def read_fasta(path):
    """Return (header, sequence) for a single-record FASTA file."""
    with open(path) as handle:
        header = handle.readline().strip()
        return header, "".join(line.strip() for line in handle)

_, mrna = read_fasta("NM_000518.5.fasta")
cds = mrna[50:494]           # the coding sequence, c.1 to c.444: start codon to stop codon

def consequence(cds, position, alt):
    """Effect of changing coding position c.<position> (1-based) to alt: (kind, protein change)."""
    number = (position - 1) // 3 + 1             # the codon number, which is also the amino acid number
    codon = cds[3 * (number - 1):3 * number]
    # TODO: 1. make the new codon by replacing the base at index (position - 1) % 3 with alt
    #       2. translate both codons, then return one of:
    #          ("synonymous", "p.His3="), ("start_lost", "p.Met1?"), ("stop_lost", "p.Ter148GluextTer?"),
    #          ("nonsense", "p.Gln40Ter") or ("missense", "p.Glu7Val"), built from THREE and number
    return "unknown", ""

examples = [("HbS", 20, "T"), ("HbC", 19, "A"), ("HbE", 79, "A"), ("Codon 39", 118, "T"),
            ("Codon 3", 9, "C"), ("Start codon", 1, "G"), ("Stop codon", 442, "G")]
for name, position, alt in examples:
    kind, change = consequence(cds, position, alt)
    print(f"{name:<12} c.{position}{cds[position - 1]}>{alt:<4} {kind:<11} {change}")

bases = "TCAG"
amino_acids = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
codon_table = {a + b + c: amino_acids[16 * i + 4 * j + k]
               for i, a in enumerate(bases)
               for j, b in enumerate(bases)
               for k, c in enumerate(bases)}

def read_fasta(path):
    """Return (header, sequence) for a single-record FASTA file."""
    with open(path) as handle:
        header = handle.readline().strip()
        return header, "".join(line.strip() for line in handle)

_, mrna = read_fasta("NM_000518.5.fasta")
coding = mrna[50:]            # from c.1 to the end of the mRNA, so translation can run past the normal stop
LAST_JUNCTION = 315           # exon 2 ends at c.315; exon 3, the last exon, starts at c.316

def translate_until_stop(seq):
    protein = ""
    for i in range(0, len(seq) - 2, 3):
        protein += codon_table[seq[i:i + 3]]
        if protein[-1] == "*":
            break
    return protein

normal = translate_until_stop(coding)
variants = {                  # c.N is index N - 1 of coding
    "c.118C>T  (codon 39, nonsense)": coding[:117] + "T" + coding[118:],
    "c.126_129del  (codons 41/42, -TTCT)": coding[:125] + coding[129:],
    "c.27dup  (codons 8/9, +G)": coding[:27] + "G" + coding[27:],
    "c.364G>T  (codon 121, nonsense in exon 3)": coding[:363] + "T" + coding[364:],
}
print("Normal β-globin:", len(normal) - 1, "amino acids, then a stop")
for name, seq in variants.items():
    protein = translate_until_stop(seq)
    first = next(i for i, (a, b) in enumerate(zip(normal, protein)) if a != b)
    stop_end = 3 * len(protein)                   # position of the new stop codon's last base, counted from its c.1
    junction = LAST_JUNCTION + len(seq) - len(coding)     # an indel before the junction moves it too
    if stop_end < junction:
        where = f"{junction - stop_end} bases before the last exon junction"
    else:
        where = "in the last exon"
    fate = "decay: no protein" if junction - stop_end > 50 else "no decay: a short protein is made"
    print(f"{name}\n    changes residue {first + 1}, stops after {len(protein) - 1} amino acids; stop codon {where}, so {fate}")

import gzip
from collections import defaultdict

population, super_population = {}, {}       # sample -> population code; population -> continental group
with open("integrated_call_samples_v3.20130502.ALL.panel") as handle:
    next(handle)                             # skip the header line
    for line in handle:
        sample, pop, super_pop = line.split("\t")[:3]
        population[sample] = pop
        super_population[pop] = super_pop

with gzip.open("1kGP_high_coverage_chr11_5225001-5230000.vcf.gz", "rt") as handle:
    for line in handle:
        if line.startswith("#CHROM"):
            samples = line.rstrip("\n").split("\t")[9:]
        elif not line.startswith("#"):
            fields = line.rstrip("\n").split("\t")
            if fields[1] == "5227002" and fields[4] == "A":          # rs334, the sickle-cell allele
                genotypes = dict(zip(samples, fields[9:]))           # sample -> genotype such as "0|1"
print(len(samples), "people in the VCF;", len(population), "unrelated people in", len(super_population), "populations")

counts = defaultdict(lambda: [0, 0, 0])     # population -> [people with 0 copies, 1 copy, 2 copies]
# TODO: for each unrelated person (each sample in population), count the copies of the allele "1"
#       in their genotype, such as "0|1", and add 1 to counts[population[sample]][copies]

def frequency(c):
    """Allele frequency from genotype counts [0 copies, 1 copy, 2 copies]."""
    return 0.0     # TODO: copies of the allele / all copies, remembering that everyone carries two

for pop, c in sorted(counts.items(), key=lambda item: -frequency(item[1])):
    if frequency(c) > 0:
        print(f"{pop} ({super_population[pop]})  {sum(c):>4} people  {c[1]:>3} with one copy  {c[2]} with two  frequency {frequency(c):.1%}")
print("No copies at all in", sum(1 for c in counts.values() if frequency(c) == 0), "populations")

afr = [sum(counts[pop][k] for pop in counts if super_population[pop] == "AFR") for k in range(3)]
q, n = frequency(afr), sum(afr)
print()
print(f"African-ancestry populations: {n} people, allele frequency {q:.4f}")
print(f"Hardy-Weinberg expects {n * 2 * q * (1 - q):.1f} people with one copy and {n * q * q:.1f} with two;")
print(f"observed: {afr[1]} and {afr[2]}")

import numpy as np
from Bio.PDB import PDBParser

structure = PDBParser(QUIET=True).get_structure("p53", "1TSR.pdb")[0]
p53 = structure["B"]                    # the copy bound to its recognition site
dna = np.array([atom.coord for chain in "EF" for residue in structure[chain] if residue.id[0] == " " for atom in residue])
print(len(dna), "DNA atoms;", sum(1 for residue in p53 if residue.id[0] == " "), "amino acids in chain B")

def closest(residue):
    """Shortest distance from any atom of the residue to any DNA atom, in ångströms."""
    return min(np.linalg.norm(dna - atom.coord, axis=1).min() for atom in residue)

for number in (175, 245, 248, 249, 273, 282):
    residue = p53[number]
    distance = closest(residue)
    kind = "touches the DNA" if distance < 4 else "holds the binding surface in shape"
    print(f"{residue.get_resname()}{number}: {distance:4.1f} Å  {kind}")

touching = [f"{r.get_resname()}{r.id[1]}" for r in p53 if r.id[0] == " " and closest(r) < 3.5]
print()
print("All residues within 3.5 Å of the DNA:", ", ".join(touching))

Mini project: annotate every HBB SNV in ClinVar

Variant annotation tools map each variant onto the transcripts it falls in and name its consequence. Build a small annotator for HBB and test it against the consequence terms ClinVar gives.

1. Read the exons of NM_000518.5 from the annotation and ClinVar’s HBB SNVs (given).
2. Write annotate(pos, alt). For a position in an exon, convert it to a c. number, remembering that HBB is on the minus strand; for a position in an intron, decide whether it hits a splice site.
3. Compare your answers with ClinVar’s MC field for every SNV (given).

import gzip
bases = "TCAG"
amino_acids = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
codon_table = {a + b + c: amino_acids[16 * i + 4 * j + k]
               for i, a in enumerate(bases)
               for j, b in enumerate(bases)
               for k, c in enumerate(bases)}

def read_fasta(path):
    """Return (header, sequence) for a single-record FASTA file."""
    with open(path) as handle:
        header = handle.readline().strip()
        return header, "".join(line.strip() for line in handle)

COMPLEMENT = str.maketrans("ACGT", "TGCA")
_, mrna = read_fasta("NM_000518.5.fasta")
cds = mrna[50:494]
UTR5 = 50                                    # the 5′ UTR of NM_000518.5 is 50 bases long

# 1. The exons of NM_000518.5, from the annotation, in transcript order (HBB is on the minus strand)
EXONS = []
with open("NC_000011.10_5220001-5300000.gff3") as handle:
    for line in handle:
        fields = line.rstrip("\n").split("\t")
        if len(fields) == 9 and fields[2] == "exon" and "Parent=rna-NM_000518.5;" in fields[8]:
            EXONS.append((int(fields[3]), int(fields[4])))
EXONS.sort(reverse=True)
print("Exons (start, end):", EXONS)

def coding_effect(position, alt):
    """Sequence Ontology term for changing coding position c.<position> to the base alt."""
    number = (position - 1) // 3 + 1
    codon = cds[3 * number - 3:3 * number]
    i = (position - 1) % 3
    before, after = codon_table[codon], codon_table[codon[:i] + alt + codon[i + 1:]]
    if before == after:
        return "synonymous_variant"
    if number == 1:
        return "initiator_codon_variant"
    if before == "*":
        return "stop_lost"
    if after == "*":
        return "nonsense"
    return "missense_variant"

# 2. ClinVar's HBB SNVs, with the consequence terms ClinVar gives in its MC field
records = []
with gzip.open("clinvar_20261019_chr11_5220001-5300000.vcf.gz", "rt") as handle:
    for line in handle:
        if line.startswith("#"):
            continue
        chrom, pos, vid, ref, alt, qual, filt, info = line.rstrip("\n").split("\t")
        info = dict(item.split("=", 1) for item in info.split(";"))
        if len(ref) == 1 and alt in ("A", "C", "G", "T") and info.get("GENEINFO", "").startswith("HBB:") and "MC" in info:
            records.append((int(pos), ref, alt, {term.split("|")[1] for term in info["MC"].split(",")}))
print(len(records), "SNVs to annotate")

# 3. Your annotator
def annotate(pos, alt):
    """Consequence term for changing chromosome position pos (plus strand) to alt, for NM_000518.5."""
    # TODO: 1. If pos is in an exon, find its c. number. Walk the exons in transcript order, adding up
    #          their lengths; on the minus strand, pos is (end - pos) bases into its exon. Subtract UTR5.
    #          c < 1 -> "5_prime_UTR_variant"; c > len(cds) -> "3_prime_UTR_variant"; otherwise
    #          coding_effect(c, alt.translate(COMPLEMENT)), because the transcript reads the minus strand.
    #       2. Otherwise pos is in an intron: "splice_donor_variant" for its first two bases (next to the
    #          exon before it in transcript order), "splice_acceptor_variant" for its last two bases,
    #          and "intron_variant" for the rest.
    return None

# 4. Compare with ClinVar
agree, disagree = 0, []
for pos, ref, alt, terms in records:
    mine = annotate(pos, alt)
    if mine in terms:
        agree += 1
    else:
        disagree.append((pos, ref, alt, mine, sorted(terms)))
print(f"Your annotation agrees with ClinVar for {agree} of {len(records)} SNVs")
for row in disagree[:5]:
    print("  disagrees:", row)

Summary and further reading

Key takeaways

• A variant is any difference from a reference. Germline variants are in every cell and can be inherited; somatic variants arise in one cell lineage, as in cancer.
• About 70 new mutations arise in each child. C→T changes at methylated CpG sites are the most common kind, which is why CpG is rare in human DNA.
• Variants range from single bases (SNVs) and short indels to structural variants of 50 bases or more. Transitions outnumber transversions, with Ti/Tv about 2 in whole genomes.
• VCF records variants on the plus strand with 1-based positions and an anchor base for indels; HGVS names them at the genome, coding and protein levels, and minus-strand genes complement the bases.
• Variants can change a protein (missense, nonsense, frameshift), how its mRNA is spliced or processed, or how much is made. Premature stops more than about 50 bases before the last exon junction usually trigger nonsense-mediated decay.
• Allele frequencies and Hardy–Weinberg proportions describe variants in populations; a variant common in healthy people is unlikely to cause a rare disease.

Key terms

Data used on this page

• Human β-globin variants: ClinVar and 1000 Genomes NCBI ClinVar and the International Genome Sample Resource (1000 Genomes)
• Human β-globin locus on chromosome 11 (GRCh38) NCBI RefSeq NC_000011.10 (GRCh38.p14) and its NCBI annotation
• Human β-globin (HBB) gene and protein NCBI RefSeq NM_000518.5 and UniProt P68871
• Protein–DNA complexes: p53, the Lac repressor and TBP RCSB Protein Data Bank entries 1TSR, 1EFA and 1CDW

References and further reading

1. 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 526, 68–74 (2015). doi:10.1038/nature15393
2. Sudmant PH, et al. An integrated map of structural variation in 2,504 human genomes. Nature 526, 75–81 (2015). doi:10.1038/nature15394
3. Byrska-Bishop M, et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell 185, 3426–3440 (2022). doi:10.1016/j.cell.2022.08.004
4. Landrum MJ, et al. ClinVar: updates to support classifications of both germline and somatic variants. Nucleic Acids Research 53, D1313–D1321 (2025). doi:10.1093/nar/gkae1090
5. Karczewski KJ, et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 (2020). doi:10.1038/s41586-020-2308-7
6. Danecek P, et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011). doi:10.1093/bioinformatics/btr330
7. Tan A, Abecasis GR, Kang HM. Unified representation of genetic variants. Bioinformatics 31, 2202–2204 (2015). doi:10.1093/bioinformatics/btv112
8. Hart RK, et al. HGVS Nomenclature 2024: improvements to community engagement, usability, and computability. Genome Medicine 16, 149 (2024). doi:10.1186/s13073-024-01421-5
9. Richards S, et al. Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the ACMG and AMP. Genetics in Medicine 17, 405–424 (2015). doi:10.1038/gim.2015.30
10. Kong A, et al. Rate of de novo mutations and the importance of father’s age to disease risk. Nature 488, 471–475 (2012). doi:10.1038/nature11396
11. Jónsson H, et al. Parental influence on human germline de novo mutations in 1,548 trios from Iceland. Nature 549, 519–522 (2017). doi:10.1038/nature24018
12. Duncan BK, Miller JH. Mutagenic deamination of cytosine residues in DNA. Nature 287, 560–561 (1980). doi:10.1038/287560a0
13. Alexandrov LB, et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013). doi:10.1038/nature12477
14. DePristo MA, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43, 491–498 (2011). doi:10.1038/ng.806
15. Thein SL. The molecular basis of β-thalassemia. Cold Spring Harbor Perspectives in Medicine 3, a011700 (2013). doi:10.1101/cshperspect.a011700
16. Hall GW, Thein S. Nonsense codon mutations in the terminal exon of the β-globin gene are not associated with a reduction in β-mRNA accumulation: a mechanism for the phenotype of dominant β-thalassemia. Blood 83, 2031–2037 (1994). doi:10.1182/blood.V83.8.2031.2031
17. Nagy E, Maquat LE. A rule for termination-codon position within intron-containing genes: when nonsense affects RNA abundance. Trends in Biochemical Sciences 23, 198–199 (1998). doi:10.1016/S0968-0004(98)01208-0
18. Romão L, et al. Nonsense mutations in the human β-globin gene lead to unexpected levels of cytoplasmic mRNA accumulation. Blood 96, 2895–2901 (2000). doi:10.1182/blood.V96.8.2895
19. Baglioni C. The fusion of two peptide chains in hemoglobin Lepore and its interpretation as a genetic deletion. PNAS 48, 1880–1886 (1962). doi:10.1073/pnas.48.11.1880
20. Allison AC. Protection afforded by sickle-cell trait against subtertian malarial infection. British Medical Journal 1, 290–294 (1954). doi:10.1136/bmj.1.4857.290
21. Cho Y, Gorina S, Jeffrey PD, Pavletich NP. Crystal structure of a p53 tumor suppressor-DNA complex: understanding tumorigenic mutations. Science 265, 346–355 (1994). doi:10.1126/science.8023157
22. Olivier M, Hollstein M, Hainaut P. TP53 mutations in human cancers: origins, consequences, and clinical use. Cold Spring Harbor Perspectives in Biology 2, a001008 (2010). doi:10.1101/cshperspect.a001008
23. The Huntington’s Disease Collaborative Research Group. A novel gene containing a trinucleotide repeat that is expanded and unstable on Huntington’s disease chromosomes. Cell 72, 971–983 (1993). doi:10.1016/0092-8674(93)90585-E
24. HGVS Nomenclature: the official recommendations, with examples of every kind of variant. hgvs-nomenclature.org