Genes, genomes and chromosomes

Why this matters for bioinformatics

A genome sequence on its own is millions or billions of letters with no punctuation. Analysis becomes possible when the letters are tied to positions and labels: this gene starts here, this exon ends there, this variant sits at that base. Those positions and labels come from two things you will use in almost every project: a reference genome and its annotation.

This topic explains how genomes are organised, from whole chromosomes down to genes, transcripts and exons, and how an annotation file describes them. You will see the human β-globin gene from the previous topic in its real setting on chromosome 11, explore DNA wrapped around histone proteins in 3D, and write Python that reads the annotation of a bacterial genome and of a human chromosome.

Genomes: the complete set of genetic material

A genome is all of an organism’s genetic material. For many viruses it is a single molecule of RNA. A bacterium usually has one circular chromosome, and many bacteria also carry small extra circles of DNA called plasmids. A human cell has 23 pairs of chromosomes in its nucleus, one set from each parent, plus many copies of a small circular genome in its mitochondria.

Virus and bacterium: counted from the NCBI annotations on the datasets page; you will reproduce the E. coli figures in the code lab. Human: from the human genome reference papers listed at the end of this topic.

The last row is the striking one. Bacterial genes sit almost end to end, with short gaps between them. In the human genome, protein-coding sequence is a small minority. The rest includes introns, untranslated regions, regulatory DNA such as promoters and enhancers, genes for RNAs that are never translated, and repeated sequences: nearly half of the human genome is derived from transposable elements, “jumping genes” that copied themselves around the genome over evolutionary time.

Chromosomes and how DNA is packaged

Laid end to end, the DNA in one human cell would stretch about 2 metres, yet it fits inside a nucleus around a hundredth of a millimetre across. It can do so because it is wound around proteins called histones. About 147 base pairs of DNA wrap around a core of eight histones, two each of H2A, H2B, H3 and H4, to form a nucleosome. Nucleosomes repeat along the DNA roughly every 200 base pairs, joined by short stretches of linker DNA, and this DNA–protein complex is called chromatin.

The anatomy of a chromosome

Each human chromosome is one long DNA molecule. Two landmarks matter for analysis. The centromere is where the two sister chromatids are held together and where the cell’s machinery attaches to pull them apart during division; it divides the chromosome into a short arm, called p, and a long arm, q. That is why genome locations are often given as bands such as 11p15.4, the band on the short arm of chromosome 11 that holds the β-globin genes. Telomeres cap the ends of each chromosome with long runs of a short repeated sequence, TTAGGG in humans, which stop the ends being treated as broken DNA.

Humans have 22 pairs of numbered chromosomes, called autosomes, and one pair of sex chromosomes, XX or XY. Centromeres and telomeres are made of long, nearly identical repeats, which made them impossible to assemble from short sequencing reads. They were only completed in 2022, when the Telomere-to-Telomere consortium published the first gapless human genome.

Bacteria have no nucleus and no histones. Their circular chromosome is compacted by other DNA-binding proteins into a region called the nucleoid. The laboratory strain E. coli K-12 MG1655, whose genome you will use in the code lab, carries no plasmids.

Genes, transcripts and proteins

A gene is a stretch of DNA that is transcribed into a functional RNA, together with the sequences that control when and where that happens. Most genes discussed in biology are protein-coding: their RNA is translated. Many others make RNAs that work as RNA, such as the tRNAs and rRNAs of translation and thousands of long non-coding RNAs whose roles are still being studied.

A gene, its RNAs and its proteins are related but not interchangeable. One gene can produce several transcripts, called isoforms, by starting transcription at different promoters, splicing exons together in different combinations or ending at different polyadenylation sites. Each transcript may encode a different protein, or none. When you count sequencing reads “per gene”, you are adding together every transcript of that gene.

Because of this, each level has its own identifiers. These are the identifiers of human β-globin, as they appear in the files you will use:

The number after the dot is a version: NM_000518.5 is the fifth version of that record, and the number changes whenever the sequence does. MANE Select marks the one transcript per protein-coding gene that NCBI and EMBL-EBI agreed to treat as the standard, so that clinical reports agree on how to number a variant.

Gene families, clusters and pseudogenes

Genes are often duplicated during evolution, and the copies can diverge. The human β-globin genes are a classic example: five working genes on chromosome 11 that encode related globin chains, each used at a different stage of life. Explore the real annotation below. Every drawing in it is read from the same GFF3 file you will open in the code lab.

Things to notice:

• Every gene here is on the minus strand. They are read from right to left, so each gene begins at its larger coordinate.
• The order matches development. Reading along the minus strand, the genes appear in the order they are switched on: HBE1 (ε-globin) in the embryo, HBG2 and HBG1 (γ-globin) in the fetus, and HBD and HBB (δ- and β-globin) after birth. The locus control region, a cluster of regulatory elements upstream of HBE1, is needed for high expression of the whole cluster.
• The working globin genes share one structure: three exons and two introns, and a 444-base coding sequence for a 147-amino-acid protein. HBB and HBD encode proteins that differ at only 10 of those 147 positions.
• HBBP1 is a pseudogene: a gene-like sequence, related to the others, that no longer makes a working globin. NCBI annotates it as transcribed, but with no coding sequence. BGLT3 is a long non-coding RNA gene.

Reading a genome annotation

An annotation lists features, such as genes, transcripts, exons and coding sequences, with their positions. The most common format is GFF3 (General Feature Format version 3): plain text, one feature per line, nine columns separated by tabs. Here is the β-globin mRNA from the file you just explored, column by column:

Features are linked by their attributes. Each line can have an ID, and a line’s Parent points to the ID of the feature it belongs to. The HBB gene (ID=gene-HBB) is the parent of the mRNA (ID=rna-NM_000518.5), which is the parent of three exon lines and three CDS lines. In bacterial annotations from NCBI, a CDS usually points straight at its gene, because there is no splicing to describe.

The phase column

A coding sequence split across exons does not always break between codons. The phase of a CDS line is the number of bases to skip at its start to reach the first complete codon. HBB’s first CDS piece is 92 bases long: 30 whole codons plus 2 bases of codon 31. The third base of codon 31 is in the next exon, so the second CDS line has phase 1. Intron 1 sits inside codon 31.

Coordinates and strands

GFF3 counts from 1 and includes both ends, as do VCF and SAM. BED files, and Python, count from 0 and leave out the end. The same six bases are therefore “4 to 9” in GFF3, “3 9” in BED and seq[3:9] in Python. Even dbSNP’s machine-readable SPDI notation counts from 0: it writes the sickle-cell position, 5,227,002, as 5227001.

Reference genomes, assemblies and versions

A reference genome is a representative sequence for a species, assembled from the DNA of one or a few individuals. It is a coordinate system, not a description of any one person: a typical human genome differs from the reference at 4 to 5 million positions. Reads from a new sample are aligned to the reference, and variants are reported as differences from it.

References improve over time, and each improvement is a new assembly with its own coordinates. The human genome has had several:

• GRCh37 (also called hg19), released in 2009 and still found in many older datasets and clinical pipelines.
• GRCh38 (hg38), released in 2013 and the current standard, with patch releases such as GRCh38.p14 that add sequence without moving the main coordinates.
• T2T-CHM13, published in 2022: the first complete, gapless sequence of the 22 autosomes and the X chromosome, including the centromeres and the short arms of five chromosomes that are missing from GRCh38.

Because the sequences differ, so do the coordinates: the sickle-cell base is at 5,227,002 on chromosome 11 in GRCh38 but at 5,248,232 in GRCh37. RefSeq makes the version explicit: GRCh38’s chromosome 11 is NC_000011.10, GRCh37’s is NC_000011.9. Annotations are versioned too: genes and transcripts are added and corrected in each NCBI annotation release or Ensembl release, and a transcript can gain a new version when its sequence changes.

Code lab: explore real genome annotations

Time to work with real annotation files. Each box below is a small Python program 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. Each box runs on its own, and your edits are saved in this browser until you press Reset.

The first five boxes use the E. coli K-12 genome and its NCBI annotation; the last two, and the mini project, use the human β-globin region of chromosome 11. Boxes marked Exercise are checked automatically, and each has a hint and a worked solution.

import gzip
from collections import Counter

with gzip.open("NC_000913.3.gff3.gz", "rt") as handle:
    lines = handle.read().splitlines()

for line in lines[:7]:
    print(line[:100])

features = [line.split("\t") for line in lines if line and not line.startswith("#")]
print()
print(f"{len(features):,} features, each with {len(features[0])} columns")
for feature_type, count in Counter(fields[2] for fields in features).most_common():
    print(f"  {feature_type:<24}{count:>6,}")

import gzip
from collections import Counter

def parse_attributes(text):
    """Turn 'ID=gene-b0344;Name=lacZ;...' into a dictionary."""
    return dict(item.split("=", 1) for item in text.split(";") if item)

genes = []
with gzip.open("NC_000913.3.gff3.gz", "rt") as handle:
    for line in handle:
        fields = line.rstrip("\n").split("\t")
        if len(fields) == 9 and fields[2] == "gene":
            info = parse_attributes(fields[8])
            genes.append({"name": info["Name"], "start": int(fields[3]), "end": int(fields[4]),
                          "strand": fields[6], "biotype": info["gene_biotype"]})

print(len(genes), "genes:", dict(Counter(gene["biotype"] for gene in genes)))
print()
for gene in genes:
    if gene["name"].startswith("lac"):
        print(gene)

import gzip
import statistics
import matplotlib.pyplot as plt

def parse_attributes(text):
    """Turn 'ID=gene-b0344;Name=lacZ;...' into a dictionary."""
    return dict(item.split("=", 1) for item in text.split(";") if item)

genes = []
with gzip.open("NC_000913.3.gff3.gz", "rt") as handle:
    for line in handle:
        fields = line.rstrip("\n").split("\t")
        if len(fields) == 9 and fields[2] == "gene":
            info = parse_attributes(fields[8])
            genes.append({"name": info["Name"], "start": int(fields[3]), "end": int(fields[4]),
                          "strand": fields[6], "biotype": info["gene_biotype"]})

coding = [gene for gene in genes if gene["biotype"] == "protein_coding"]

lengths = []        # TODO: the length of each gene in coding, in the same order
median = None       # TODO: the median of lengths
longest = None      # TODO: the gene dictionary (from coding) with the greatest length

print(len(lengths), "protein-coding genes; median length", median, "bp")
print("Longest:", longest)

plt.figure(figsize=(7, 3.2))
plt.hist(lengths, bins=80, color="#0fa391")
plt.xlabel("Gene length (bp)")
plt.ylabel("Number of genes")
plt.show()

import gzip

GENOME_LENGTH = 4_641_652        # from the FASTA file, and the region line of the GFF3

intervals = []
with gzip.open("NC_000913.3.gff3.gz", "rt") as handle:
    for line in handle:
        fields = line.rstrip("\n").split("\t")
        if len(fields) == 9 and fields[2] == "CDS" and "pseudo=true" not in fields[8]:
            intervals.append((int(fields[3]), int(fields[4])))

def merge(intervals):
    """Merge overlapping (start, end) intervals. Both ends are included."""
    merged = []
    for start, end in sorted(intervals):
        # TODO: if this interval overlaps the last one in merged, extend that one;
        #       otherwise, append it as a new interval.
        pass
    return merged

merged = merge(intervals)
coding_bases = sum(end - start + 1 for start, end in merged)
print(len(intervals), "CDS pieces merge into", len(merged), "stretches")
print(f"{coding_bases:,} of {GENOME_LENGTH:,} bases code for protein: {coding_bases / GENOME_LENGTH:.1%}")

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 translate(dna):
    """Translate codon by codon, stopping at the first stop codon."""
    protein = ""
    for i in range(0, len(dna) - 2, 3):
        amino_acid = codon_table[dna[i:i + 3]]
        if amino_acid == "*":
            break
        protein += amino_acid
    return protein

def read_fasta(path):
    """Return (header, sequence) for a single-record FASTA file, gzipped or not."""
    opener = gzip.open if path.endswith(".gz") else open
    with opener(path, "rt") as handle:
        header = handle.readline().strip()
        return header, "".join(line.strip() for line in handle)

def reverse_complement(seq):
    return seq.translate(str.maketrans("ACGT", "TGCA"))[::-1]

_, genome = read_fasta("NC_000913.3.fasta.gz")
start, end, strand = 363231, 366305, "-"        # lacZ, from the GFF3 file

gene = None      # TODO: slice out bases start..end (1-based, both included), then
                 #       reverse-complement them, because lacZ is on the minus strand

protein = translate(gene or "")
print(len(gene or ""), "bases, starting", (gene or "")[:3], "->", len(protein), "amino acids")
print(protein[:50], "...")

_, uniprot = read_fasta("P00722.fasta")
print("Matches UniProt P00722:", protein == uniprot)

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)

def parse_attributes(text):
    """Turn 'ID=gene-b0344;Name=lacZ;...' into a dictionary."""
    return dict(item.split("=", 1) for item in text.split(";") if item)

header, region = read_fasta("NC_000011.10_5220001-5300000.fasta")
print(header)
offset = int(header.split(":")[1].split("-")[0])     # chromosome position of region[0]
print(f"{len(region):,} bases; chromosome position p is region[p - {offset:,}]")
print()

genes = []
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] in ("gene", "pseudogene"):
            info = parse_attributes(fields[8])
            genes.append((int(fields[3]), int(fields[4]), fields[6], info["Name"], info["gene_biotype"]))

print("Genes in the order the minus strand reads them, right to left:")
for start, end, strand, name, biotype in sorted(genes, reverse=True):
    print(f"  {name:<9}{start:>11,}{end:>11,}  {strand}  {end - start + 1:>6,} bp  {biotype}")

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

# HBB's first CDS line (exon 1) in the GFF3 file
cds_start, cds_end, strand = 5226930, 5227021, "-"

c1_position = None      # TODO: chromosome position of c.1 (on the minus strand, which end is the start?)
c20_position = None     # TODO: chromosome position of c.20
plus_base = None        # TODO: the base at c20_position, as stored in region (the plus strand)
gene_base = None        # TODO: the same base read from HBB's strand: its complement

print("c.1 is at", c1_position, "and c.20 is at", c20_position)
print("Plus strand:", plus_base, "  HBB's strand:", gene_base)

Mini project: splice HBB from chromosome 11

Do what the cell does, in code. HBB’s coding sequence is split across three exons on the minus strand. Find its pieces in the annotation, cut them out of the chromosome, join them and translate the result, then compare it with the β-globin protein in UniProt.

1. Collect the CDS lines of the transcript NM_000518.5 (given).
2. Cut each piece out of the region file, converting chromosome coordinates with the offset, and join them in chromosome order.
3. Reverse-complement the joined sequence, because HBB is on the minus strand.
4. Translate it and check the protein against UniProt P68871.

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 translate(dna):
    """Translate codon by codon, stopping at the first stop codon."""
    protein = ""
    for i in range(0, len(dna) - 2, 3):
        amino_acid = codon_table[dna[i:i + 3]]
        if amino_acid == "*":
            break
        protein += amino_acid
    return protein

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)

def reverse_complement(seq):
    return seq.translate(str.maketrans("ACGT", "TGCA"))[::-1]

header, region = read_fasta("NC_000011.10_5220001-5300000.fasta")
offset = 5220001                    # the region file starts at chromosome position 5,220,001

# 1. Collect the coding pieces (CDS lines) of the HBB transcript NM_000518.5.
pieces = []
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] == "CDS" and "Parent=rna-NM_000518.5;" in fields[8]:
            pieces.append((int(fields[3]), int(fields[4]), fields[6]))
print("CDS pieces:", pieces)

# 2. Cut out each piece (chromosome position p is region[p - offset]) and join them,
#    in chromosome order (smallest start first).
joined = None       # TODO

# 3. HBB is on the minus strand: reverse-complement to read it 5' to 3'.
cds = None          # TODO

# 4. Translate and compare with UniProt.
protein = translate(cds or "")
_, uniprot = read_fasta("P68871.fasta")
print(len(cds or ""), "bases ->", len(protein), "amino acids")
print(protein)
print("Matches UniProt P68871:", protein == uniprot)

Summary and further reading

Key takeaways

• A genome is all of an organism’s genetic material. Bacteria usually have one circular chromosome; humans have 23 pairs of chromosomes plus mitochondrial DNA.
• In eukaryotes DNA is packaged as chromatin: 147 bp wrap around eight histones in each nucleosome. Chromosomes have a centromere and telomeres.
• A gene is a transcribed region plus its regulatory DNA. One gene can make several transcripts, and genes, transcripts and proteins each have versioned identifiers.
• Gene density varies enormously: about 86% of the E. coli genome codes for protein, less than 2% of the human genome.
• GFF3 describes each feature in nine columns and links features with ID and Parent. Its coordinates are 1-based and inclusive; BED and Python are 0-based and exclude the end.
• Minus-strand features begin at their end coordinate and must be reverse-complemented. Every coordinate belongs to one assembly, such as GRCh38.

Key terms

Data used on this page

• Human β-globin locus on chromosome 11 (GRCh38) NCBI RefSeq NC_000011.10 (GRCh38.p14) and its NCBI annotation
• Escherichia coli K-12 MG1655 genome and annotation NCBI RefSeq NC_000913.3 and UniProt P00722 and P02920
• Human β-globin (HBB) gene and protein NCBI RefSeq NM_000518.5 and UniProt P68871
• DNA, RNA and nucleosome 3D structures RCSB Protein Data Bank entries 1BNA, 2DCG, 1EHZ and 1KX5
• SARS-CoV-2 reference genome NCBI RefSeq NC_045512.2

References and further reading

1. Luger K, et al. Crystal structure of the nucleosome core particle at 2.8 Å resolution. Nature 389, 251–260 (1997). doi:10.1038/38444
2. Davey CA, et al. Solvent mediated interactions in the structure of the nucleosome core particle at 1.9 Å resolution. Journal of Molecular Biology 319, 1097–1113 (2002). doi:10.1016/S0022-2836(02)00386-8
3. Kornberg RD. Chromatin structure: a repeating unit of histones and DNA. Science 184, 868–871 (1974). doi:10.1126/science.184.4139.868
4. Blattner FR, et al. The complete genome sequence of Escherichia coli K-12. Science 277, 1453–1462 (1997). doi:10.1126/science.277.5331.1453
5. Jacob F, Monod J. Genetic regulatory mechanisms in the synthesis of proteins. Journal of Molecular Biology 3, 318–356 (1961). doi:10.1016/S0022-2836(61)80072-7
6. International Human Genome Sequencing Consortium. Initial sequencing and analysis of the human genome. Nature 409, 860–921 (2001). doi:10.1038/35057062
7. Schneider VA, et al. Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly. Genome Research 27, 849–864 (2017). doi:10.1101/gr.213611.116
8. Nurk S, et al. The complete sequence of a human genome. Science 376, 44–53 (2022). doi:10.1126/science.abj6987
9. The 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 526, 68–74 (2015). doi:10.1038/nature15393
10. Morales J, et al. A joint NCBI and EMBL-EBI transcript set for clinical genomics and research (MANE). Nature 604, 310–315 (2022). doi:10.1038/s41586-022-04558-8
11. Grosveld F, et al. Position-independent, high-level expression of the human β-globin gene in transgenic mice. Cell 51, 975–985 (1987). doi:10.1016/0092-8674(87)90584-8
12. Stamatoyannopoulos G. Control of globin gene expression during development and erythroid differentiation. Experimental Hematology 33, 259–271 (2005). doi:10.1016/j.exphem.2004.11.007
13. The GFF3 specification, maintained by the Sequence Ontology project. github.com