Molecules of life

What a cell is made of

Lesson 1 measured the cell: how big it is, how many molecules it holds, and how fast things happen inside it. This lesson opens the box and looks at what those molecules are. By the end you will know the main kinds of molecule in a cell and why they behave the way they do in water, and you will have written Python that reads the real DNA of syn3A.

Start with the whole cell. The widget below breaks one E. coli cell down in two ways: by mass, and by the number of molecules. Switch between the views and tap the parts.

What a cell is made of

One E. coli cell: about 1 µm³ in volume and about 1 pg (a trillionth of a gram) in mass. Tap any part to see its numbers.

By massBy count

Most of a cell is water; most of the rest is protein.

Water70%

Dry 30%

The dry 30%, enlarged

ProteinProtein55%

RNARNA20.4%

LipidsLipids9.1%

Small molecules and ionsSmall molecules4.0%GlycogenGlycogen2.5%

Cell wall and outer-membrane sugarsWall and sugars5.9%DNADNA3.1%

Water

not part of the dry mass70% of the whole cellabout 0.7 pg

The liquid that everything else is dissolved in or floats in. Water also takes part in many of the cell’s reactions.

Values for E. coli growing fast, rounded. Source: Neidhardt, Ingraham and Schaechter, Physiology of the Bacterial Cell (1990), as summarized in Milo and Phillips, Cell Biology by the Numbers. syn3A has no cell wall or outer membrane, so it has no peptidoglycan or lipopolysaccharide.

Some things to notice:

• By mass, a cell is mostly water: about 70%. More than half of the rest is protein.
• By count, small molecules win by a mile. Water molecules outnumber proteins about 10,000 to 1. Sugars, amino acids, salts and the other small molecules outnumber proteins about 100 to 1.
• The big molecules are few, but they do most of the work. One cell has only about 19,000 ribosomes and two copies of its DNA.

Atoms: the parts list

Every molecule is made of atoms joined by chemical bonds. There are over a hundred elements, but six of them make up almost all of a living cell: carbon (C), hydrogen (H), nitrogen (N), oxygen (O), phosphorus (P), and sulfur (S). Biologists remember them as CHNOPS. About half of a bacterium’s dry mass is carbon. A cell also needs small amounts of a few others, such as potassium, sodium, magnesium, and iron, mostly as ions dissolved in water.

How many bonds each atom makes

Each kind of atom forms a typical number of bonds with its neighbors, called its valence. That number decides which shapes a molecule can take, a bit like the number of ports on a connector.

Atom	Symbol	Bonds it usually makes	Mass of one atom (Da)
Hydrogen	H	1	1.008
Oxygen	O	2	15.999
Nitrogen	N	3	14.007
Carbon	C	4	12.011
Sulfur	S	2	32.06
Phosphorus	P	5, in phosphate	30.974

A double bond counts twice. In carbon dioxide, written O=C=O, the carbon has two double bonds: 2 + 2 = 4.

Math toolkit: formulas and molecular mass

A molecular formula lists the atoms in a molecule. Each element’s symbol is followed by how many of that atom there are, and no number means one. Glucose is $\text{C}_6\text{H}_{12}\text{O}_6$: 6 carbons, 12 hydrogens, and 6 oxygens. Water is $\text{H}_2\text{O}$: 2 hydrogens and 1 oxygen. Methane is $\text{CH}_4$: 1 carbon and 4 hydrogens.

A molecule’s mass is the sum of the masses of its atoms. Recall from Lesson 1 that molecular masses are measured in daltons (Da), and that a molecule of mass $M$ daltons weighs $M$ grams per mole.

Worked example: the mass of glucose. Multiply each atom’s mass by how many there are, then add.

1. Carbon: $6 \times 12.011 = 72.066$ Da
2. Hydrogen: $12 \times 1.008 = 12.096$ Da
3. Oxygen: $6 \times 15.999 = 95.994$ Da
4. Add them up: $72.066 + 12.096 + 95.994 = 180.156$ Da

So glucose is about 180 Da, and one mole of it weighs about 180 grams.

Worked example: how many molecules are in a teaspoon of sugar? A teaspoon of sugar weighs about 4 g. Table sugar is really sucrose, but treating it as glucose keeps the numbers simple.

1. Moles: divide the mass by the mass of one mole, $\frac{4\ \text{g}}{180\ \text{g/mol}} \approx 0.022$ mol.
2. Molecules: multiply by Avogadro’s number, $0.022 \times 6.022 \times 10^{23} \approx 1.3 \times 10^{22}$.

About $10^{22}$ molecules in one spoonful: ten thousand billion billion.

Now look at real molecules. Each model in Molecule Zoo is built from atom positions in PubChem, a public chemistry database. Tap any atom to see how many bonds it makes. Then press Play: count the atoms.

3D game · Molecule Zoo

Real molecules, atom by atom

Sound on

Small moleculesSugarsAmino acidsDNA and RNA basesEnergy and fat

GlucoseRibose

Each model is built from a real 3D structure. Click any atom to see how many bonds it makes.

Play: count the atoms

Glucose · 24 atoms

Reset view

Drag to turn · Scroll to zoom · Click an atom

Glucose

Sugars

C₆H₁₂O₆24 atoms · about 0.65 nm from end to end

6 carbon12 hydrogen6 oxygen

Mass: 6 × C (12.011) + 12 × H (1.008) + 6 × O (15.999) = 180.16 Da

The cell’s main fuel; breaking it down releases energy.

A dalton (Da) is about the mass of one hydrogen atom. 3D coordinates: computed model from PubChem CID 5793.

Bonds each atom makes:H hydrogen 1O oxygen 2N nitrogen 3C carbon 4S sulfur 2P phosphorus 5 (in phosphate)

Your first Python of the lesson turns a formula like C6H12O6 into a mass, using the atomic masses in the engine’s constants.

Python exercise

Molecular mass from a formula

import re # regular expressions: patterns for finding things in text from virtualcell.constants import ATOMIC_MASS # {"H": 1.008, "C": 12.011, ...} in daltons def molecular_mass(formula): """Mass of one molecule in daltons (Da), from a formula such as "C6H12O6".""" # Step 1: find the (element, count) pairs. This line does it: # pairs = re.findall(r"([A-Z][a-z]?)(\d*)", formula) # ([A-Z][a-z]?) is an element symbol; (\d*) is the digits after it, if any. # For "CH4" it gives [("C", ""), ("H", "4")]. # Step 2: each count is text. An empty count means 1 atom; int() turns the rest into numbers. # Step 3: add up count × ATOMIC_MASS[element] for all the pairs, and return the total. raise NotImplementedError("Write your code here") # Remove the # below to try it. Glucose should be about 180 Da. # print(molecular_mass("C6H12O6"))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfWhat is the molecular mass of water, H₂O?

Two hydrogens and one oxygen: $2 \times 1.008 + 15.999 = 18.015$ Da. So a mole of water weighs about 18 grams.

Check yourselfMethane is CH₄. Why does it have four hydrogens?

Carbon makes 4 bonds and hydrogen makes 1. A single carbon needs four partners, so four hydrogens complete it.

Check yourselfHow many atoms are in one ATP molecule, C₁₀H₁₆N₅O₁₃P₃?

Add the counts: $10 + 16 + 5 + 13 + 3 = 47$ atoms.

Water: where everything happens

Water is about 70% of a cell’s mass, and almost every reaction in the cell happens in it. Water has one property that shapes all of biology: it is polar.

In a water molecule, the oxygen pulls the shared electrons toward itself. That leaves the oxygen end slightly negative, written δ− (“delta minus”), and the two hydrogen ends slightly positive (δ+). Opposite charges attract, so the hydrogen of one water molecule is drawn to the oxygen of the next. This weak attraction is called a hydrogen bond.

Molecules that are polar or charged, such as sugars, salts, and many amino acids, can form hydrogen bonds with water, so they dissolve. They are called hydrophilic, “water-loving”. Oily molecules, made mostly of carbon and hydrogen, can’t form hydrogen bonds with water. They are hydrophobic, “water-fearing”.

How concentrated is water? A liter of water weighs about 1,000 g, and a mole of water weighs 18.015 g. So pure water is $\frac{1{,}000}{18.015} \approx 55.5$ moles per liter: 55.5 molar. A typical protein is at a micromolar or less, so water molecules outnumber it by tens of millions to one.

Math toolkit: energy, heat, and kT

Bonds hold molecules together, and heat shakes them apart. To compare the two, we need one unit of energy for both.

Chemists measure bond energies in kilojoules per mole (kJ/mol): the energy needed to break a whole mole of those bonds. Breaking a mole of carbon–carbon bonds takes about 350 kJ, so a C–C bond is about 350 kJ/mol. A hydrogen bond is far weaker: about 5 to 20 kJ/mol, depending on the atoms involved.

Everything in a cell jiggles with heat. In water, a molecule is bumped by its neighbors trillions of times a second, and the typical energy of that jiggling is written kT: Boltzmann’s constant $k$ times the temperature $T$. For a mole of molecules, the same amount is $R \times T$, where $R$ is the gas constant:

$$kT \text{ per mole} = R \times T, \qquad R = 0.008314\ \text{kJ per mol per kelvin}$$

The temperature must be in kelvin: add 273.15 to the temperature in degrees Celsius. At body temperature, 37 °C, $T = 310.15$ K, so

$$kT = 0.008314 \times 310.15 \approx 2.58\ \text{kJ/mol}$$

Now divide each energy by kT to see how many “jiggles” it is worth:

What	Energy	In units of kT, at 37 °C
The jiggling of heat	2.58 kJ/mol	1 kT
A hydrogen bond	5–20 kJ/mol	about 2–8 kT
The energy released by using one ATP in a cell	about 50 kJ/mol	about 20 kT
A carbon–carbon bond	about 350 kJ/mol	about 136 kT

How often does heat break a bond? Physics gives a simple answer. The chance that random jiggling gathers an energy $E$ in one place is proportional to the Boltzmann factor:

$$\text{Boltzmann factor} = e^{-E/kT}$$

The letter $e$ stands for a number, like $\pi$ does: about 2.718. It is the base of natural growth and decay, and every calculator and programming language has it. In Python, math.exp(x) means $e^x$. You don’t need to know where $e$ comes from. You need two facts about $e^{-x}$:

• When $x$ is 0, $e^{0} = 1$. No energy is needed, so it always happens.
• Every extra kT of energy needed makes it about 2.7 times less likely, and every extra 2.3 kT makes it 10 times less likely.

Try it with the table:

• A hydrogen bond of 10 kJ/mol is $\frac{10}{2.58} \approx 3.9$ kT. Its Boltzmann factor is $e^{-3.9} \approx 0.02$: about 1 in 50.
• A carbon–carbon bond of 350 kJ/mol is about 136 kT. Its Boltzmann factor is $e^{-136}$.

To turn $e^{-x}$ into a power of ten, divide $x$ by 2.303 (because $e^{2.303} \approx 10$). So $e^{-136} = 10^{-136/2.303} \approx 10^{-59}$. A molecule jiggles about $10^{13}$ times a second, so heat would take about $10^{46}$ seconds to break that bond: vastly longer than the age of the universe.

Python exercise

How often does heat break a bond?

import math # math.exp(x) is e to the power of x from virtualcell.constants import BODY_TEMPERATURE_K # 37 °C is 310.15 K from virtualcell.constants import GAS_CONSTANT_KJ # R = 0.008314 kJ per mol per kelvin def boltzmann_factor(energy_kj_per_mol, temperature_k): """How likely heat alone is to supply an energy E: exp(-E / (R × T)).""" # Step 1: work out R × T, the thermal energy in kJ/mol. (It is about 2.58 at 37 °C.) # Step 2: divide the energy by R × T. Return e to the power of minus that. raise NotImplementedError("Write your code here") # Remove the # below to try it: a 10 kJ/mol hydrogen bond at body temperature. # print(boltzmann_factor(10, BODY_TEMPERATURE_K))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

The four families of big molecules

The big molecules of a cell fall into four families. Each is built from its own kind of building block.

Family	Building block	Kinds of block	Main jobs
Nucleic acids (DNA and RNA)	nucleotides	4	store and carry genetic information
Proteins	amino acids	20	almost everything: enzymes, pumps, motors, structure
Carbohydrates	sugars	many	fuel, energy stores, cell walls
Lipids	fatty acids and others	many	membranes and energy stores

In the first three families, breaking a link between two building blocks uses up one water molecule, a reaction called hydrolysis. That is how cells, and your own gut, take these molecules apart.

Nucleic acids: DNA and RNA

DNA (deoxyribonucleic acid) stores the cell’s genes. Its building block, a nucleotide, has three parts: a phosphate, a sugar called deoxyribose, and one of four bases: adenine (A), thymine (T), guanine (G), or cytosine (C). The phosphates and sugars link into a backbone, and the bases stick out to one side. The order of the bases along the strand is the genetic information.

A DNA molecule is two strands wound around each other: the double helix. The bases of one strand pair with the bases of the other, held by hydrogen bonds, and they follow two rules:

• A pairs with T, held by 2 hydrogen bonds.
• G pairs with C, held by 3 hydrogen bonds.

So if you know one strand, you know the other. The partner strand is called the complement.

The two strands run in opposite directions. Each strand has a 5′ end (“five prime”) and a 3′ end (“three prime”), named after the positions of atoms on the sugar, and sequences are always written from 5′ to 3′. Because the strands run in opposite directions, the partner of a strand, written in its own 5′ to 3′ direction, is the reverse complement: pair each base, then read the result backward.

Worked example: the reverse complement of 5′-ATGC-3′.

1. Pair each base: A–T, T–A, G–C, C–G. The partner strand is 3′-TACG-5′.
2. Read it from its 5′ end, which is on the right: 5′-GCAT-3′.

The double helix is about 2 nm wide. Each base pair adds 0.34 nm to its length, and it makes one full turn about every 10.5 base pairs. syn3A’s genome has 543,379 base pairs, so stretched out straight it would be

$$543{,}379 \times 0.34\ \text{nm} \approx 185{,}000\ \text{nm} = 185\ \mu\text{m}$$

That is about 460 times the width of the cell it has to fit inside.

RNA (ribonucleic acid) is DNA’s working copy. It uses a slightly different sugar, ribose, and the base uracil (U) in place of thymine. U pairs with A. RNA is usually a single strand. Cells copy genes into RNA, and ribosomes read the RNA to build proteins: that is the subject of Lesson 3.

Build some DNA yourself. In Zip the DNA, a template strand is waiting. Add the partner bases one at a time. Later levels make you think about direction, build RNA, and end with the first bases of syn3A’s first gene.

3D game · Zip the DNA

Build the partner strand

Streak 0Stars 0/15

Sound on

1. Pairing rules2. Against the clock3. Read it 5′ to 3′4. Make RNA5. Real genome

1. Pairing rules

The template strand is already in place. Build its partner, one base at a time, from left to right.

A pairs with T, held by 2 hydrogen bonds. G pairs with C, held by 3.

Tap a base below. On a keyboard, press A, T, G or C.

Start

Template

3′········5′

5′3′

New strand building from the left →

AadenineTthymineGguanineCcytosine

Drawn to scale: the helix is about 2 nm wide, rises 0.34 nm per base pair and turns once every 10.5 base pairs.

Because a G–C pair has three hydrogen bonds and an A–T pair only two, the share of G and C in a stretch of DNA matters. It is called the G–C content:

$$\text{G–C content} = \frac{\text{number of G} + \text{number of C}}{\text{length}}$$

Genomes differ a lot. E. coli’s genome is about 51% G–C. syn3A’s is much lower, as you will measure yourself later in this lesson.

Heat can separate DNA’s two strands. It breaks the hydrogen bonds between the bases, but not the strong covalent bonds of the backbones. This is called melting. In Melt the DNA, heat a helix whose left half is A–T-rich and whose right half is G–C-rich, and watch which half comes apart first.

3D game · Melt the DNA

How much heat pulls the strands apart?

Challenges 0/4

Sound on

1. Melt only the A–T half2. Melt it all3. Zip it back up4. Design it

Challenge: Heat the DNA until the A–T half comes apart (over 80% open) while the G–C half stays zipped (under 20% open). Hold it there for 2 seconds.

A–T half: 0% open

G–C half: 0% open

37 °C · kT = 2.58 kJ/mol

A–T-rich half (the start of syn3A’s dnaA gene) · 0% open

ATTAGCATATCGGCTAATATATCGGCATTAATTATATATA

G–C-rich half · 0% open

GCCGGCGCCGGCCGCGTAGCCGGCGCCGCGGCCGATGCCG

Top strand 5′ → 3′, left to right. Open pairs are drawn apart.

Temperature: 37 °C (310 K)2037 (body)6080100

Cool to 37 °C

The temperatures in this model are rough. Real DNA melts somewhere between about 60 °C and 95 °C, depending on its G–C content and the salt around it.

Base pairs open over the last 30 seconds

100%0%

A–T half G–C half Temperature (20–100 °C)

Tug of war, per base pair, at 37 °C

A–T half holds with 26.6 kT

G–C half holds with 28.2 kT

The white line is the pull of freedom: an open pair can wiggle in many more ways, worth about 24.1 kT at any temperature. When a bar falls short of the line, those pairs open.

Heat makes kT bigger, so the same bonds are worth fewer kT. One hydrogen bond is 3.8 kT here. One backbone bond is about 350 kJ/mol, or 136 kT, so heat never breaks the strands.

How this model works

Each base pair is either open or closed. A closed pair is held by its hydrogen bonds (9.74 kJ/mol each: 2 for A–T, 3 for G–C) and by stacking onto each closed neighbor (46.7 kJ/mol). An open pair gains freedom worth 0.2 kJ/mol for every kelvin of temperature.

Many times a second the computer picks a pair, or a short run of pairs, and tries flipping it. If the flip lowers the energy it is kept. If it raises the energy by ΔF, it is kept with probability e^−ΔF/kT. This is the Metropolis Monte Carlo method. Because stacking rewards closed neighbors, whole stretches open and close together, like a zipper.

Real DNA melts over a wider range of temperatures than this model, and its stacking energy depends on which bases sit next to each other. The model keeps the main ideas: G–C holds tighter than A–T, and neighbors open and close together.

Check yourselfWhat are the complement and the reverse complement of 5′-AACG-3′?

Pair each base, keeping the order: the complement is 3′-TTGC-5′. Read that from its 5′ end and you get the reverse complement, 5′-CGTT-3′.

Check yourselfHow many hydrogen bonds hold the base pairs A–T, T–A, G–C and C–G together?

$2 + 2 + 3 + 3 = 10$. Each A–T pair has 2 and each G–C pair has 3.

Check yourselfWhat is the G–C content of ATGCGC?

There are 4 G or C bases out of 6, so $\frac{4}{6} \approx 0.67$, or 67%.

Math toolkit: counting possibilities

How much information can a sequence hold? Count the possible sequences. A single base can be any of 4 letters. Two bases can be $4 \times 4 = 16$ different pairs. Each extra base multiplies the count by 4 again, so a sequence of $n$ bases can be any of

$$4^n \text{ different sequences}$$

• 10 bases: $4^{10} = 1{,}048{,}576$, about a million.
• 20 bases: $4^{20} \approx 1.1 \times 10^{12}$, about a trillion.

For long sequences the numbers are too big to write out, so use logarithms from Lesson 1. $4^n$ is 4 multiplied by itself $n$ times, and multiplying numbers adds their logs, so its log is $\log_{10}(4)$ added up $n$ times:

$$\log_{10}(4^n) = n \times \log_{10}(4) = n \times 0.602$$

For syn3A’s genome, $543{,}379 \times 0.602 \approx 327{,}000$. So there are about $10^{327{,}000}$ possible sequences that long, and syn3A’s is one of them. Proteins use a 20-letter alphabet, and $\log_{10}(20) \approx 1.30$, so even a short protein of 100 amino acids could be any of about $10^{130}$ sequences. For comparison, there are about $10^{80}$ atoms in the observable universe.

Proteins

Proteins do most of the work in a cell. Each is a chain of amino acids, typically a few hundred long: about 300 on average, as you saw in Lesson 1. Every amino acid has the same backbone part, which joins to its neighbors through a peptide bond, plus one of 20 different side chains. The side chains make the amino acids different: some are oily, some are polar, some carry a positive or negative charge, and two contain sulfur. One of those, cysteine, can link to another cysteine through its sulfur atom.

A chain of amino acids doesn’t stay floppy. In water it folds into a particular shape, and the shape decides what the protein can do. The main force behind folding is the hydrophobic effect again. The chain folds so that most of its oily side chains end up buried in the middle, away from the water, and the polar ones face outward. Hydrogen bonds and other weak bonds then hold the shape in place.

Fold it turns folding into a puzzle, using a famous simplified model by the biophysicist Ken Dill. Each amino acid is either H (hydrophobic, oily) or P (polar), and the chain lies on a grid. Fold it so the H beads touch. The best folds are rare: the 20-bead chain in the last level has about 42 million possible shapes, and only 2 of them are best.

Game · Fold it

First fold

Sound on

1. First fold2. Hide the core3. Oil drop4. Stack it up5. Zip it up6. The benchmark

Fold the chain so the oily H beads touch each other and hide from the water. Tap a highlighted square to place the next bead; the chain can’t cross itself. Each glowing dashed line is a contact: two H beads side by side that aren’t neighbors along the chain.

6 beads · H means hydrophobic (oily) and P means polar. This is Ken Dill’s HP model from 1985, the classic simple model of folding. Start with six beads.

NextHPPHH

Contacts 0 / best 2H sides touching water 4

UndoRestartShow water

Drag along squares to place several beads, and drag back to take them off. Arrow keys work too; Backspace undoes.Drag your finger from the last bead to place several beads. Drag back to take them off.

Lipids and membranes

Lipids are the oily molecules of the cell. The ones in membranes are mostly phospholipids: a water-loving head, which contains a phosphate, attached to two oily tails. Each tail is a fatty acid, like the palmitic acid in Molecule Zoo. A molecule with a water-loving part and an oily part is called amphipathic.

Put phospholipids in water and nothing has to arrange them. The hydrophobic effect pushes the tails together, out of the water, and the heads face the water. Lipids with one tail are shaped like cones, and they gather into small balls called micelles. Phospholipids, with two tails, are shaped more like cylinders, and they pack side by side into a flat double layer: a lipid bilayer, about 4 to 5 nm thick. A bilayer closes up into a hollow bubble, so that no oily edge touches the water. That bubble is a cell membrane.

Watch it happen in Hide the tails, a simulation of lipids drifting in water. Nothing in it tells the lipids where to go. There is only an attraction between oily tails, and the random jiggling of heat.

Simulation game · Hide the tails

Lipids that build themselves

Challenges 0/3

Sound on

1. Hide the tails2. Build a bilayer3. Melt it

Challenge: These lipids have one tail each. Watch them, and wait until fewer than 15% of the tail beads touch water.

Tails touching water–

Starting…

• Water-loving head
• Oily tail
• Tail touching water

A tail bead counts as touching water when a water molecule could fit beside it, or when its lipid is floating alone. Each bead is about 1 nm across; a real bilayer is about 4–5 nm thick (scale bar: 5 nm). Drag on the box to stir.To stir, turn on Stir with finger, then drag.

Lipid type

One tailTwo tails

Speed

1×3×

Hydrophobic pull: 3.0 kTnonestrong

ShakeStir with finger

kT is the energy of random thermal motion. A weaker pull has the same effect as more heat.

How many lipids does syn3A’s membrane need? A bilayer holds about 5 million lipids per square micrometer, counting both layers. syn3A is a sphere about 0.4 µm across, so its radius is 0.2 µm and its surface area is

$$4\pi r^2 = 4 \times 3.14 \times (0.2\ \mu\text{m})^2 \approx 0.50\ \mu\text{m}^2$$

So its membrane holds about $0.50 \times 5$ million $\approx$ 2.5 million lipids. Like other mycoplasmas, syn3A has no cell wall. This one membrane, 4 to 5 nm thick, is all that separates the cell from the world outside.

Sugars and small molecules

The rest of the cell’s molecules are small: sugars, amino acids, nucleotides, vitamins, and salts. Small molecules that take part in the cell’s chemistry are called metabolites. In E. coli, the metabolites add up to about 300 millimolar, and the most common one is glutamate, an amino acid, at about 100 mM.

One small molecule deserves special mention: ATP (adenosine triphosphate), the cell’s energy currency. It is an adenine base and a ribose sugar with a chain of three phosphates attached. Splitting off the last phosphate releases about 50 kJ/mol in a living cell, roughly 20 kT: enough to drive a reaction that wouldn’t happen on its own. Cells use ATP to power almost everything, and they remake it constantly from the food they take in. You will follow that energy in Lesson 5.

Reading DNA in code

Biologists store DNA sequences in text files, and two formats are everywhere: FASTA and GenBank. You will write Python that reads both, and use it to measure syn3A’s real genome.

FASTA files

A FASTA file is the simplest format: a header line that starts with >, followed by the sequence, split across many lines. Here are the first lines of syn3A’s genome, as stored in the GenBank database under accession number CP016816.2:

>CP016816.2 Synthetic bacterium JCVI-Syn3A, complete genome
ATGAACGTAAACGATATTTTAAAAGAACTTAAACTAAGTTTAATGGCTAATAAAAATATTGATGAATCCG
TGTATAACGACTATATAAAGACAATAAATATTCATAAAAAGGGGTTTTCTGATTATATTGTTGTTGTTAA
ATCACAATTTGGTTTGTTAGCTATAAAACAGTTTCGTCAAACTATTGAAAATGAGATAAAAAATATTTTA

The whole file is about 7,800 lines long. A FASTA file can hold many records, one after another, and the first word of each header is the record’s ID.

Start with two small functions you will use again and again, then write the FASTA reader. The last check of the first exercise runs your code on syn3A’s whole genome, which the exercise downloads from this site.

Python exercise

G–C content of a sequence

def gc_content(sequence): """Fraction of the bases in a DNA sequence that are G or C, from 0 to 1.""" # Step 1: make the letters uppercase, so "g" counts the same as "G". # Step 2: count the Gs and the Cs. A string's .count("G") counts its Gs. # Step 3: divide by the number of bases, len(sequence). # An empty sequence has no bases: return 0.0 for it. raise NotImplementedError("Write your code here") # Remove the # below to try it. Half of these bases are G or C. # print(gc_content("ATGC")) # Then try a real genome: all 543,379 bases of syn3A. # from virtualcell.genomes import load_syn3a_genome # print(gc_content(load_syn3a_genome()))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Python exercise

Reverse complement

PAIR = {"A": "T", "T": "A", "G": "C", "C": "G"} # each base and its partner def reverse_complement(sequence): """The other strand of a DNA sequence, read in its own direction.""" # Step 1: make the sequence uppercase, so lowercase input works too. # Step 2: swap each base for its partner. PAIR["A"] is "T". # Step 3: reverse the order. In Python, text[::-1] is text backwards. # Return one string. "".join(letters) turns a list of letters into a string. raise NotImplementedError("Write your code here") # Remove the # below to try it. You should see GCAT. # print(reverse_complement("ATGC"))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Python exercise

Read a FASTA file

from virtualcell.genomes import read_data_text # reads the course's data files def parse_fasta(text): """Read FASTA text into a dict of {record id: sequence}.""" # A FASTA file looks like this: # >CP016816.2 Synthetic bacterium JCVI-Syn3A, complete genome # ATGAACGTAAACGATATTTTAAAAGAACTTAAACTAAGTTTAATGGCTAATAAAAATATTGATGAATCCG # TGTATAACGACTATATAAAGACAATAAATATTCATAAAAAGGGGTTTTCTGATTATATTGTTGTTGTTAA # A line that starts with ">" begins a new record. Its first word, without the ">", is the id. # The lines after it, up to the next ">", are the sequence: join them into one uppercase string. # Skip blank lines, and strip the spaces and line breaks off each line. raise NotImplementedError("Write your code here") # Remove the #s below to read syn3A's real genome file. # records = parse_fasta(read_data_text("syn3a.fasta")) # for record_id, sequence in records.items(): # print(record_id, len(sequence), sequence[:20])

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

GenBank files

A GenBank file holds the same sequence plus what is known about it: where each gene starts and stops, what its protein is called, and who sequenced it. The sequence comes last, after a line that says ORIGIN, in numbered blocks of ten letters. Here is part of the record for syn3A’s first gene, dnaA. Its protein starts each round of copying the chromosome.

     gene            1..1353
                     /gene="dnaA"
                     /locus_tag="JCVISYN3A_0001"
     CDS             1..1353
                     /gene="dnaA"
                     /product="Chromosomal replication initiator protein"
                     /translation="MNVNDILKELKLSLMANKNIDESVYNDYIKTINIHKKGFSDYIV
                     VVKSQFGLLAIKQFRQTIENEIKNILKEPVNISFTYEQEYKKQLEKDELINKDHSDII
ORIGIN
        1 atgaacgtaa acgatatttt aaaagaactt aaactaagtt taatggctaa taaaaatatt
       61 gatgaatccg tgtataacga ctatataaag acaataaata ttcataaaaa ggggttttct

The record says the gene runs from base 1 to base 1,353. The /translation line is the protein the gene makes, written with one letter per amino acid. In Lesson 3 you will make that protein from the DNA yourself.

Python exercise

Challenge: read a GenBank record

from virtualcell.genomes import read_data_text # reads the course's data files # The first 1,500 bases of syn3A's genome, as a GenBank record. Have a look: record = read_data_text("syn3a_dnaA.gb") lines = record.splitlines() print("\n".join(lines[:4])) # the top: name, length, description... print("...") print("\n".join(lines[-3:])) # ...and the end: the last lines of sequence, then // # print(record) # remove the # to see the whole record (95 lines) def genbank_sequence(text): """The DNA sequence from a GenBank record, as one uppercase string.""" # The sequence starts on the line after the one that begins with "ORIGIN", # and ends at the line "//". Each sequence line looks like this: # 1 atgaacgtaa acgatatttt aaaagaactt aaactaagtt taatggctaa taaaaatatt # Keep only the letters (ch.isalpha() is True for a letter ch), # join all the lines, and make the result uppercase. raise NotImplementedError("Write your code here") # Remove the # below to try it. You should see 1500, as on the LOCUS line. # print(len(genbank_sequence(record)))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Code milestone: chemistry and sequences

Two new modules join the engine in this lesson:

• virtualcell/chemistry.py holds molecular_mass and boltzmann_factor.
• virtualcell/sequences.py holds gc_content, reverse_complement, parse_fasta, and genbank_sequence.

Two more pieces are written for you. virtualcell/genomes.py reads the course’s real data files, including syn3A’s genome, and constants.py now has atomic masses and the gas constant.

To finish the milestone, copy each function you wrote into its module, replacing the stub, then run the checks from the engine folder:

python -m pytest

When they pass, the engine can read a genome and analyze it:

from virtualcell.genomes import load_syn3a_genome
from virtualcell.sequences import gc_content, reverse_complement

genome = load_syn3a_genome()
print(len(genome))                      # 543379
print(round(gc_content(genome), 4))     # 0.2428
print(reverse_complement(genome[:20]))  # AAAATATCGTTTACGTTCAT

Numbers to remember

Quantity	Value
A cell’s mass	about 70% water; more than half of the rest is protein
The elements of life	C, H, N, O, P, S
Bonds each atom makes	H 1, O 2, N 3, C 4
Glucose, $\text{C}_6\text{H}_{12}\text{O}_6$	about 180 Da
Pure water	55.5 M
kT at 37 °C	2.58 kJ/mol
Bond strengths	hydrogen bond 2–8 kT; carbon–carbon bond about 136 kT
Using one ATP	about 20 kT
Base pairs	A–T, 2 hydrogen bonds; G–C, 3
DNA double helix	2 nm wide, 0.34 nm per base pair, 10.5 base pairs per turn
Lipid bilayer	4–5 nm thick
syn3A’s genome	543,379 base pairs, about 24% G–C

Next lesson: the central dogma, or how a gene becomes a protein. You will find syn3A’s genes in its GenBank record, copy one into RNA, and translate it into its protein.

Sources and further reading

• Milo R, Phillips R. Cell Biology by the Numbers. Garland Science, 2015. Free online at book.bionumbers.org.
• Neidhardt FC, Ingraham JL, Schaechter M. Physiology of the Bacterial Cell: A Molecular Approach. Sinauer, 1990. The composition of E. coli in the cell recipe.
• Alberts B, et al. Molecular Biology of the Cell, 6th edition. Garland Science, 2014. Chemical bonds, water, and the number of lipids in a bilayer.
• Bennett BD, et al. Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nature Chemical Biology 5:593–599 (2009).
• Breuer M, et al. Essential metabolism for a minimal cell. eLife 8:e36842 (2019).
• Glass JI. Synthetic bacterium JCVI-Syn3A, complete genome. GenBank CP016816.2.
• Kim S, et al. PubChem 2023 update. Nucleic Acids Research 51:D1373–D1380 (2023). The 3D structures in Molecule Zoo.
• Watson JD, Crick FHC. Molecular structure of nucleic acids. Nature 171:737–738 (1953).
• Lau KF, Dill KA. A lattice statistical mechanics model of the conformational and sequence spaces of proteins. Macromolecules 22:3986–3997 (1989). The HP model in Fold it.
• Jumper J, et al. Highly accurate protein structure prediction with AlphaFold. Nature 596:583–589 (2021).
• Cooke IR, Kremer K, Deserno M. Tunable generic model for fluid bilayer membranes. Physical Review E 72:011506 (2005). The kind of lipid model in Hide the tails.
• Blattner FR, et al. The complete genome sequence of Escherichia coli K-12. Science 277:1453–1462 (1997).