Reading Genomes

Genetics and Heredity Fundamentals 2 showed that a gene is a stretch of DNA written in four letters, and Genetics and Heredity Fundamentals 3 showed how a cell reads those letters three at a time. In this chapter we do the reading ourselves. A DNA sequence is just a string, so Python's string tools, dictionaries and counters are enough to measure it, search it and compare it with another one. All the programs here are Python 3.

What Sequencing Gives You

Sequencing is reading the order of the bases in a piece of DNA. Current machines do not read a whole chromosome in one go. They read many short pieces, called reads, and software then fits the overlapping pieces back together into longer sequences, a step called assembly (check). The reads can be hundreds of letters long, or much longer with some newer methods (check). The result is stored as text.

The most common plain-text format is FASTA (check). Each record starts with a line that begins with > followed by a name and a free-text description. The following lines hold the sequence, usually wrapped at a fixed width. The parser below keeps the first word of each name, joins the wrapped lines and forces upper case, so that a and A are treated alike.

The two records come back as a dictionary from name to sequence, each already joined into a single string. A real genome file is the same idea at a huge scale: the human genome is roughly three billion letters (check), which is why real tools read files line by line rather than loading everything at once.

Counting Letters and GC Content

The simplest question about a sequence is what it is made of. collections.Counter tallies the letters. A particularly useful summary is the GC content: the fraction of letters that are G or C. Because G pairs with C by three hydrogen bonds and A with T by two (Genetics and Heredity Fundamentals 2), GC content is related to how tightly the two strands hold together, and it differs between species and between regions of a genome (check).

The same program also builds the reverse complement. The two strands of the double helix run in opposite directions, so the partner strand, read in its own direction, is the complement written backwards. str.maketrans and str.translate swap letters in one pass, and [::-1] reverses the result.

There are 11 G or C letters out of 20, so the GC content is 0.55. Taking the reverse complement twice returns the original, a handy sanity check for the function. Reverse complements matter because a gene can sit on either strand, so a search for a pattern usually has to look at both.

Finding Patterns

Searching a sequence for a short pattern, or motif, is the everyday job. Python gives two tools. str.find(pattern, start) returns the position of the first match at or after start, or -1. The re module finds every match, with a catch: a normal re.finditer does not report matches that overlap, so in AGGGGA it finds GG only twice (positions 1 and 3). Wrapping the pattern in a zero-width lookahead, (?=GG), reports every position where the pattern begins, including overlapping ones.

A famous motif is the start codon ATG (Genetics and Heredity Fundamentals 3). A stretch that begins with ATG and continues codon by codon until a stop codon (TAA, TAG or TGA) with no earlier stop is called an open reading frame, or ORF, a candidate for a protein-coding gene (check). The function below tries every ATG, steps along in threes, and records the stretch if it reaches a stop and has at least a minimum number of codons before it. It only looks at the strand it is given. A full gene finder would repeat the search on the reverse complement and be far more careful (check).

Go through the output in order. The first line shows find returning 0, then 12 once we start searching from position 1, then -1 for a missing pattern. The second line shows all three ATG positions (0, 12 and 20). The third and fourth lines show the overlap difference: [1, 3] against [1, 2, 3]. The last line holds the two ORFs: the whole sequence from position 0 to 24 (eight codons, ending in the stop TGA), and a shorter one that starts at 12, in the same reading frame, and ends at the same stop. The ATG at 20 gives no ORF because only four letters remain, so no stop codon can be reached in its frame.

Two of 20 positions differ, so the sequences are 90% identical. The list gives the position (counting from 0) and the two letters. The first change, C to T at position 5, is in the third position of codon 2 (GCC becomes GCT), and both codons code for the same amino acid, alanine (check), a kind of change called silent (Genetics and Heredity Fundamentals 3). The second change, G to C at position 14, turns the codon ATG (methionine) into ATC (isoleucine), which does change the amino acid (check).

The Hamming distance has a weakness. If one letter is deleted, every letter after it shifts by one place, and the two strings look different almost everywhere. The edit distance (also called Levenshtein distance) counts the fewest single-letter insertions, deletions and substitutions needed to turn one string into the other, so it copes with shifts (check). It is computed by filling a table, row by row, keeping only the previous row in memory.

The first line is the textbook example, "kitten" to "sitting", which needs 3 edits. In the second line, the distance between our two DNA strings is 2, matching the Hamming distance. The last two numbers compare a deletion: ACGTACGT with the T removed becomes ACGACGT, so the edit distance is 1, while the Hamming distance (over the shorter length) is 4. Real aligners use refinements of this table, with different costs for different kinds of change (check).

k-mers and Sliding Windows

A k-mer is any run of k consecutive letters. Counting all the k-mers of a sequence is a fast way to summarise it, and it is the basis of many tools (check). A sequence of length n has n - k + 1 k-mers, which overlap. A sliding window is the same trick applied to a measurement: slide a fixed-width window along the sequence and compute something, such as GC content, at each step. It shows where along a sequence a property rises and falls.

The 2-mer AT occurs 5 times in the 14-letter string and the 3-mers ATA and TAT tie at 3. The window output reads as a profile: GC content is 0 at the start, climbs to 1.0 in the middle where the sequence is all GC repeats, and drops back to 0 at the end. In real genomes, dips and peaks in GC content can point to regions with a different history or function, though the meaning depends on the genome (check).

Challenges