DNA, RNA and the central dogma

Why this matters for bioinformatics

Almost every file you will meet as a bioinformatician is a description of DNA, RNA or protein. A FASTA file is a strand written as letters. A FASTQ file is millions of short pieces of strands read by a sequencer. A VCF file lists the places where one genome’s letters differ from another’s.

To use those files correctly you need a working picture of the molecules behind the letters. Many analysis bugs come from forgetting which strand a sequence is on, which direction it is written in, or where a gene really starts. This topic builds that picture, using real structures and real genomes rather than cartoons alone.

Nucleotides: the building blocks

DNA (deoxyribonucleic acid) and RNA (ribonucleic acid) are polymers: long chains built from repeating units called nucleotides. Every nucleotide has three parts:

1. A phosphate group, which carries a negative charge. It is why DNA moves towards the positive electrode in gel electrophoresis.
2. A five-carbon sugar: deoxyribose in DNA, ribose in RNA. Its carbons are numbered 1′ to 5′ (“one-prime” to “five-prime”). The primes distinguish sugar carbons from the atoms of the base.
3. A nitrogenous base, attached to the 1′ carbon. The base is the part that carries information.

There are five common bases. Adenine (A) and guanine (G) are purines, with two fused rings. Cytosine (C), thymine (T) and uracil (U) are pyrimidines, with a single ring. DNA uses A, C, G and T; RNA uses A, C, G and U.

Strands have a direction: 5′ and 3′

Nucleotides join through a phosphodiester bond: the phosphate on the 5′ carbon of one nucleotide links to the hydroxyl on the 3′ carbon of the previous one. Repeating this makes a sugar–phosphate backbone with the bases sticking out to one side.

Because every link runs from a 3′ carbon to a 5′ carbon, a strand has two different ends. The 5′ end finishes with a free phosphate; the 3′ end finishes with a free hydroxyl. Enzymes that make DNA and RNA can only add new nucleotides to the 3′ end, so strands always grow in the 5′ → 3′ direction.

In a double helix the two strands run in opposite directions: they are antiparallel. Convention says we always write a sequence from 5′ to 3′. So the partner of 5′-ATGC-3′ is written 5′-GCAT-3′: complement each base, then reverse the order. This is the reverse complement, and you will compute it constantly.

Base pairing and the double helix

The two strands are held together by hydrogen bonds between bases. The pairing is specific: A pairs with T through two hydrogen bonds, and G pairs with C through three. These are called Watson–Crick base pairs.

Base pairing explains a pattern Erwin Chargaff measured around 1950, before the structure was known. In double-stranded DNA the amount of A equals the amount of T, and G equals C. If a DNA sample is 30% A, it must also be 30% T, which leaves 40% for G and C together: 20% each.

Try it below. Type a sequence and watch the partner strand, the hydrogen bonds and the GC content change.

The shape: a right-handed double helix

In 1953 James Watson and Francis Crick proposed that the two antiparallel strands wind around a common axis, with the bases paired inside and the backbones outside. Their model relied on X-ray diffraction images taken by Rosalind Franklin and Raymond Gosling, published alongside it in the same issue of Nature. The common form in cells, B-DNA, is about 2 nm wide, rises 0.34 nm per base pair and makes one full turn about every 10.5 base pairs.

RNA: same alphabet, different chemistry

RNA differs from DNA in three ways that matter:

• Sugar. Ribose has a hydroxyl (–OH) on its 2′ carbon; deoxyribose (“de-oxy”) has only a hydrogen there. The 2′-OH makes RNA more chemically reactive and much less stable, which suits a short-lived working copy.
• Base. RNA uses uracil (U) where DNA uses thymine (T). U pairs with A just as T does.
• Shape. RNA is usually made as a single strand. It folds back on itself, pairing short complementary stretches into stems and loops, and these folds give many RNAs their function.

Cells make several kinds of RNA. Messenger RNA (mRNA) carries the instructions for a protein. Transfer RNA (tRNA) brings amino acids to the ribosome. Ribosomal RNA (rRNA) forms the core of the ribosome and catalyses the joining of amino acids. Many other RNAs, such as microRNAs and long non-coding RNAs, regulate genes without ever being translated.

Below is a real tRNA. Its 76 nucleotides fold into an L shape held together by base pairs between distant parts of the chain, including G·U “wobble” pairs and unusual pairs that DNA never makes. Notice how much of the molecule is double-stranded even though it is one chain.

How does a single strand decide which parts to pair? The simplest model looks for the arrangement with the most base pairs. Try it with the sequences below: complementary stretches separated by a loop fold into a hairpin.

The central dogma: DNA → RNA → protein

In 1958 Francis Crick summarised how genetic information flows, and restated it in 1970 as the central dogma of molecular biology. Information is copied from DNA to DNA (replication), from DNA to RNA (transcription) and from RNA to protein (translation). Some viruses add special routes, but information never flows back out of a protein into nucleic acid.

Transcription

RNA polymerase opens the double helix and reads one strand, the template strand, from 3′ to 5′. It builds an RNA from 5′ to 3′ that is complementary to the template. That RNA therefore has the same sequence as the other strand, the coding strand, except that it has U where the coding strand has T. Sequence databases show genes as coding strand, 5′ to 3′.

Translation

A ribosome reads the mRNA three bases at a time. Each triplet is a codon. Four bases make 4 × 4 × 4 = 64 codons: 61 specify amino acids and 3 (UAA, UAG, UGA) mean “stop”. Translation usually begins at AUG, which also codes for methionine. Most amino acids have several codons, so the code is redundant (or degenerate): a change in the third base of a codon often leaves the protein unchanged.

Where you start reading matters. A sequence can be read in three reading frames on each strand, and only one of them, usually, gives the real protein.

One letter changes everything. In sickle-cell disease, base 20 of the HBB coding sequence is T instead of A. Codon 7 changes from GAG to GTG, and the protein gets valine instead of glutamate: a missense change written c.20A>T (p.Glu7Val). Older papers call it Glu6Val, because they count from the second amino acid: the first methionine is removed from the mature protein. Vernon Ingram identified this single amino-acid difference in 1957, the first time an inherited disease was traced to one amino-acid change in a protein.

Genes are more than coding sequence

In plants, animals and fungi, most genes are split. The pieces that remain in the mature mRNA are exons; the pieces in between, which are removed by splicing, are introns. The mRNA also keeps untranslated regions (UTRs) at both ends, and gains a 5′ cap and a poly(A) tail. HBB, for example, has three exons and two introns.

Code lab: your first genome analysis

Time to work with real data. Each box below is a small Python program that runs in your browser. Press Run (or Ctrl + Enter) to run it. The first run downloads Python, about 12 MB, which takes a few seconds; after that, runs are instant. Each box runs on its own, so you can do them in any order, and your edits are saved in this browser until you press Reset.

Boxes marked Exercise are checked automatically when you run them. If you get stuck, there is a hint and a worked solution.

with open("NC_045512.2.fasta") as handle:
    header = handle.readline().strip()               # the first line
    sequence = "".join(line.strip() for line in handle)  # all the rest, joined

print("Header:", header)
print("Length:", len(sequence), "bases")
print("First 60 bases:", sequence[:60])
print("Last 60 bases: ", sequence[-60:])

from collections import Counter

with open("NC_045512.2.fasta") as handle:
    handle.readline()                                    # skip the header
    sequence = "".join(line.strip() for line in handle)

counts = Counter(sequence)
for base in "ACGT":
    print(f"{base}: {counts[base]:>6,}  ({counts[base] / len(sequence):.1%})")

gc = (counts["G"] + counts["C"]) / len(sequence)
print(f"GC content: {gc:.2%}")

def reverse_complement(seq):
    """Return the reverse complement of a DNA sequence written 5'->3'."""
    pairs = {"A": "T", "T": "A", "G": "C", "C": "G"}
    # 1. Complement every base using the dictionary.
    # 2. Reverse the result so it reads 5'->3'.
    return ""   # replace this line


print(reverse_complement("ATGC"))     # should print GCAT
print(reverse_complement("GAATTC"))   # the EcoRI site: should print GAATTC

# The standard genetic code, built from the usual T, C, A, G ordering.
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")
cds = mrna[50:494]          # bases 51-494, as given in the GenBank record

def translate(dna):
    """Translate codon by codon, stopping at the first stop codon."""
    protein = ""
    # Step through dna three bases at a time: range(0, len(dna) - 2, 3)
    # Look each codon up in codon_table; stop (break) when you reach "*".
    return protein

protein = translate(cds)
print(len(cds), "bases ->", len(protein), "amino acids")
print(protein)

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

# The standard genetic code, built from the usual T, C, A, G ordering.
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)

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

_, mrna = read_fasta("NM_000518.5.fasta")
cds = mrna[50:494]

position = 20                                    # c.20 is base 20 of the CDS
print("Base at c.20:", cds[position - 1])       # should be A

mutant = cds          # TODO: replace base 20 with T (build a new string with slices)

normal_protein = translate(cds)
mutant_protein = translate(mutant)

codon_number = None   # TODO: which codon (counting from 1) contains base 20?
before = None         # TODO: that amino acid in normal_protein
after = None          # TODO: that amino acid in mutant_protein
print(codon_number, before, "->", after)

import gzip
from collections import Counter
import matplotlib.pyplot as plt

def read_sequence(path):
    opener = gzip.open if path.endswith(".gz") else open
    with opener(path, "rt") as handle:
        handle.readline()
        return "".join(line.strip() for line in handle)

genomes = {
    "E. coli K-12 (double-stranded DNA)": read_sequence("NC_000913.3.fasta.gz"),
    "SARS-CoV-2 (single-stranded RNA)": read_sequence("NC_045512.2.fasta"),
}

fig, ax = plt.subplots(figsize=(7, 3.2))
for offset, (name, seq) in enumerate(genomes.items()):
    counts = Counter(seq)
    print(f"{name}: {len(seq):,} bases, A/T = {counts['A'] / counts['T']:.3f}, G/C = {counts['G'] / counts['C']:.3f}")
    fractions = [counts[base] / len(seq) for base in "ACGT"]
    ax.bar([i + offset * 0.4 for i in range(4)], fractions, width=0.4, label=name)

ax.set_xticks([i + 0.2 for i in range(4)], list("ACGT"))
ax.set_ylabel("Fraction of genome")
ax.legend(frameon=False, fontsize=8)
plt.show()

from Bio import SeqIO

record = SeqIO.read("NM_000518.5.gb", "genbank")
print(record.id, "-", record.description)
print("Length:", len(record.seq), "bases\n")

for feature in record.features:
    if feature.type == "exon":
        print("exon", feature.location)
    elif feature.type == "CDS":
        cds = feature.extract(record.seq)
        protein = cds.translate(to_stop=True)
        print("CDS ", feature.location, "->", len(cds), "bases")
        print("Protein:", protein)
        print("Matches the annotation:", str(protein) == feature.qualifiers["translation"][0])

Mini project: find and translate the spike gene

The spike protein lets SARS-CoV-2 enter human cells, and it is the target of COVID-19 vaccines. Use the genome and its GFF3 annotation to pull out the spike gene, translate it and locate one of its best-known features: the furin cleavage site, the amino acids PRRAR, which is cut by a human enzyme as part of activating the protein.

1. Find the coordinates of the gene named S in the GFF3 file (given).
2. Extract its DNA from the genome, converting 1-based, inclusive coordinates to a Python slice.
3. Translate it, leaving out the stop codon.
4. Find where PRRAR begins, numbering amino acids from 1.

# The standard genetic code, built from the usual T, C, A, G ordering.
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)

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

_, genome = read_fasta("NC_045512.2.fasta")

# 1. Find the spike (S) gene in the annotation. GFF3 columns are tab-separated:
#    seqid, source, type, start, end, score, strand, phase, attributes
start = end = None
with open("NC_045512.2.gff3") as handle:
    for line in handle:
        fields = line.rstrip("\n").split("\t")
        if len(fields) == 9 and fields[2] == "gene" and "Name=S;" in fields[8] + ";":
            start, end = int(fields[3]), int(fields[4])
print("Spike gene:", start, "to", end)

# 2. Extract the gene's DNA (1-based inclusive -> 0-based slice).
spike_dna = None        # TODO

# 3. Translate it, without the stop codon.
spike_protein = None    # TODO

# 4. Where does the furin site PRRAR start (counting amino acids from 1)?
furin_position = None   # TODO

print(len(spike_protein or ""), "amino acids; PRRAR at", furin_position)

Summary and further reading

Key takeaways

• A nucleotide is a phosphate, a sugar (deoxyribose in DNA, ribose in RNA) and a base. Bases carry the information.
• Strands have direction. They grow 5′ → 3′, pair antiparallel, and are always written 5′ → 3′; the partner strand is the reverse complement.
• A pairs with T (or U) through two hydrogen bonds, G with C through three. B-DNA is a right-handed helix about 2 nm wide with about 10.5 base pairs per turn.
• RNA is usually single-stranded and folds into functional shapes, such as the L-shaped tRNA.
• Genes are transcribed into RNA, and mRNA is translated codon by codon into protein. One base change can change the protein, as in sickle-cell disease.
• Coordinates need care: annotation files are usually 1-based and inclusive, Python slices are 0-based and exclude the end.

Key terms

Data used on this page

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

References and further reading

1. Watson JD, Crick FHC. Molecular structure of nucleic acids: a structure for deoxyribose nucleic acid. Nature 171, 737–738 (1953). doi:10.1038/171737a0
2. Franklin RE, Gosling RG. Molecular configuration in sodium thymonucleate. Nature 171, 740–741 (1953). doi:10.1038/171740a0
3. Crick F. Central dogma of molecular biology. Nature 227, 561–563 (1970). doi:10.1038/227561a0
4. Ingram VM. Gene mutations in human haemoglobin: the chemical difference between normal and sickle cell haemoglobin. Nature 180, 326–328 (1957). doi:10.1038/180326a0
5. Drew HR, et al. Structure of a B-DNA dodecamer: conformation and dynamics. PNAS 78, 2179–2183 (1981). doi:10.1073/pnas.78.4.2179
6. Wang AH-J, et al. Molecular structure of a left-handed double helical DNA fragment at atomic resolution. Nature 282, 680–686 (1979). doi:10.1038/282680a0
7. Shi H, Moore PB. The crystal structure of yeast phenylalanine tRNA at 1.93 Å resolution: a classic structure revisited. RNA 6, 1091–1105 (2000). doi:10.1017/S1355838200000364
8. PDB-101 Molecule of the Month: DNA. A short, illustrated introduction from the RCSB Protein Data Bank. pdb101.rcsb.org