59  Regular expressions

A regular expression is a tiny program that describes a pattern in text. The re module is the machine that runs it.

import sandbox_widget

You already know how to look for a piece of text inside a string. The in operator tells you whether it is there, .find tells you where, and .count tells you how many times. If you want to know whether a DNA sequence contains a site for the restriction enzyme EcoRI, which cuts DNA wherever it sees GAATTC, that is all you need:

'GAATTC' in 'CCGAATTCATGG'

Now take the enzyme HinfI. It cuts at GANTC, where N means any base. 'GANTC' in dna will not help you, because it looks for a real letter N, and there is none in your DNA. You could check GAATC, GACTC, GAGTC and GATTC one at a time, but that gets old fast. The TATA box, a short sequence found in many promoters, has the consensus TATAWAW, where W means A or T. That is four strings to check. And if you want to find a poly-A tail, a run of As of any length, there is no list of strings you can write down at all.

What you need is a way to describe a pattern rather than spell out every string that fits it. That is what a regular expression is. It is a small language of its own, with its own grammar, and it does exactly one thing: it instructs a little machine to look for text that fits a description. A course called Instructing Machines could hardly leave it out.

What is a regular expression?

A regular expression (people often just say regex) is a description of a set of strings. The pattern GAATTC describes exactly one string. The HinfI site GANTC describes four DNA strings, one for each base in the middle. And the pattern for a poly-A tail describes infinitely many: A, AA, AAA, and so on forever. You write the description, and the re module checks text against it.

That lets you ask four kinds of questions, and you will meet all of them in this chapter:

  • Searching: is there anything in this text that fits the description, and where?
  • Validating: does the whole string fit? Is this string a DNA sequence at all?
  • Extracting: which pieces fit, and what are the parts inside them? Pull the chromosome, start and end out of chr7:114086327-114693772.
  • Replacing: swap everything that fits for something else, for instance to strip the numbers and spaces out of a sequence copied from a GenBank file.

The idea is older than most programming languages. The mathematician Stephen Kleene invented the notation in the 1950s, and in the late 1960s Ken Thompson built it into a text editor, and later into the command-line tool grep, whose name stands for globally search for a regular expression and print. Since then regular expressions have found their way into almost every programming language and text editor. What you learn here also works, with small differences in the details, in R, on the command line, and in the search box in VS Code (click the .* button).

Handing the job to an engine

Here is the part that makes regular expressions more than a convenience. When you write a for-loop that walks along a sequence looking for something, Python carries out every single step of it: take the next index, cut out a slice, compare it, add one to a counter, go round again. Python reads and carries out your instructions one at a time, and each of those small steps costs time.

When you call re.search, none of that walking happens in Python. The re module translates your pattern into a compact set of instructions for a regular expression engine: a small, specialized machine, written in the much faster language C, whose only job is to run patterns over text. Engines like it have been tuned and optimized by programmers for decades, because searching text is something computers are asked to do all the time. Your one line of Python hands the whole job over, and the engine does the walking.

Exercise 59-1

SOLO

The code below counts the EcoRI sites in a sequence of 1.8 million bases twice: first with a for-loop that compares every six-base slice to GAATTC, then with re.findall, which you will meet properly in ?sec-re-findall. The time module from the standard library lets you read the clock before and after each one. Do both ways find the same number of sites? And which do you think is faster, and by how much: twice as fast, ten times, a hundred times? Decide before you run it.

%%sandbox
import re
import time

dna = 'ATGCGAATTCTTAGCCAT' * 100000

start = time.time()
count = 0
for i in range(len(dna) - 5):
    if dna[i:i+6] == 'GAATTC':
        count = count + 1
print('Loop:', count, 'sites in', time.time() - start, 'seconds')

start = time.time()
count = len(re.findall('GAATTC', dna))
print('re:  ', count, 'sites in', time.time() - start, 'seconds')

On my laptop, the engine is about forty times faster than the loop. Your numbers will be different, and they will change a little every time you run the cell, but the engine wins by a wide margin. On a whole genome that is the difference between waiting a second and waiting a minute.

One fair warning before you get carried away: string methods such as .find and .count are written in C too, so for exact text like GAATTC they are just as fast as re. The engine earns its keep when you are looking for a pattern rather than a fixed string, which is what the rest of this chapter is about. And speed is never an excuse for a pattern nobody can read (see ?sec-re-when-not).

Characters with special powers

Most characters in a pattern just match themselves: a G matches a G. A handful of characters are special. They are called metacharacters, and they are what makes a regular expression more than a string.

The first one is the dot. A . matches any single character. That is exactly the N in the HinfI site GANTC:

%%sandbox
import re

re.search('GA.TC', 'TTGACTCAA')

Exercise 59-6

SOLO

For each of the four searches below, decide whether it finds a match and, if it does, what text it matches. Write your four answers down before you run the code.

%%sandbox
import re

print(re.search('GA.TC', 'GATTCGAGTC'))
print(re.search('GA.TC', 'GATC'))
print(re.search('GA.TC', 'GA.TC'))
print(re.search('GA..C', 'AAGAATTCAA'))

The second one is the one people get wrong: the dot has to match exactly one character. It cannot match nothing.

Character classes

The dot is often too generous. In the TATA box consensus TATAWAW, the W is not any base, it is A or T. Square brackets let you list the characters allowed at a single position. [AT] matches one character that is either an A or a T. This is called a character class:

%%sandbox
import re

promoter = 'GCGCTATAAAAGGCGC'
re.search('TATA[AT]A[AT]', promoter)

A character class always matches exactly one character, however many characters you list inside the brackets. A few more tricks go inside the brackets:

  • [A-Z] is a range: any uppercase letter. [0-9] is any digit.
  • [^GC] with a ^ first is the opposite: any character that is not G or C.

Some classes are so common that they have their own shorthand. \d matches any digit, like [0-9]. \s matches any whitespace (a space, a tab or a newline). \w matches any word character: a letter, a digit or an underscore.

pattern matches one
G G
. any character
[AT] A or T
[A-Z] uppercase letter
[^GC] anything but G, C
\d digit
\s whitespace
\w letter, digit or _

Exercise 59-7

SOLO

Predict what each line prints, then run the code to check.

%%sandbox
import re

print(re.search('[^ACGT]', 'ACGTACGNACGT'))
print(re.search(r'\d', 'chromosome 7'))
print(re.search(r'\s', 'FOXP2 gene'))
print(re.search('[GC][GC][GC]', 'ATGCATGGCATT'))

Why the little r?

You probably noticed that the patterns with a backslash in them were written as r'\d' instead of '\d'. In the chapter on files, you saw that a backslash already means something in an ordinary Python string: '\n' is a single newline character, not a backslash followed by an n. Regular expressions use backslashes too, and the two uses get in each other’s way. Putting an r in front of the string makes it a raw string, where Python leaves every backslash alone and passes it on untouched to the re module.

Some patterns happen to work without the r and some silently do something else, and you do not want to have to remember which. Always write your patterns as raw strings. Just trust me on that one.

When you mean the dot itself

What if you want to match a character that happens to be a metacharacter, like a real dot? You put a backslash in front of it. \. matches a dot and nothing else, \( matches an opening parenthesis, and so on. This is called escaping the character.

Exercise 59-8

SOLO

ImportantSpot the bug

The pattern in the cell below is wrong on purpose. Your job is to find out how.

The code below is supposed to report which lines mention the FOXP2 transcript NM_148900.4. Only the first line does. The third line mentions NM_14890014, a completely different accession number. Run the code and see that it reports the third line anyway. Explain why, then fix the pattern so that only the first line matches.

%%sandbox
import re

lines = ['Reads aligned to NM_148900.4',
         'Reads aligned to NM_148900.3',
         'Reads aligned to NM_14890014']
for line in lines:
    print(line, re.search('NM_148900.4', line))

How many?

So far, every piece of a pattern matched exactly one character. Quantifiers change that. A quantifier goes right after the thing it applies to and says how many times it may repeat:

  • + means one or more times,
  • * means zero or more times,
  • ? means zero or one time (i.e., it is optional),
  • {3} means exactly three times, and {2,5} means between two and five times.
pattern means
A+ one or more A
A* zero or more A
A? zero or one A
A{3} exactly three A
A{2,5} two to five A
A{4,} four or more A

Here is the poly-A tail that no list of strings could describe:

%%sandbox
import re

mrna = 'AUGGCCUAAGCGAAAAAAAAAAAA'
re.search('A+', mrna)

Hmm. It matched a single A, the one in AUGGCC. That is correct, though: A+ means one or more As, and the very first A in the string is one or more As. The search always reports the leftmost place where the pattern fits. If you want a real tail, say what you mean by one, for instance at least six As:

%%sandbox
import re

mrna = 'AUGGCCUAAGCGAAAAAAAAAAAA'
re.search('A{6,}', mrna)

When a quantifier has a choice, it takes as many characters as it can. Once it found six As, it kept going and took all twelve. Quantifiers are greedy.

Exercise 59-9

SOLO

Predict exactly what text each search matches before you run the code. The quantifier applies only to the single character right in front of it.

%%sandbox
import re

print(re.search('CA*G', 'TTCAAAGTT'))
print(re.search('CA*G', 'TTCGTT'))
print(re.search('CA+G', 'TTCGTT'))
print(re.search('colou?r', 'the color of the gel'))
print(re.search('[ACGT]{3}', 'NNATGNN'))

Exercise 59-10

SOLO

Greed can surprise you. The pattern below is meant to find a stretch that starts with ATG and ends with the stop codon TAA. The sequence contains two TAAs. Which one do you think the match ends at? Decide, then run it.

%%sandbox
import re

re.search('ATG.+TAA', 'ATGAAATAAGGGTAA')

The .+ gobbled up as much as it could and only gave characters back when it had to, so the match ends at the last TAA. If you put a ? after a quantifier, it becomes lazy and takes as few characters as it can instead. Try changing the pattern to 'ATG.+?TAA' and run it again.

(Before you get ideas: finding real open reading frames takes more than this, because the stop codon has to be in the same reading frame as the start codon. You will do it properly in the ORF project.)

Anchors

Two metacharacters match no character at all. They match a position. ^ matches the start of the string and $ matches the end. Together they let you insist that a pattern describes the whole string and not just some piece of it.

Exercise 59-11

SOLO

The pattern ^[ACGT]+$ reads from the start, one or more of A, C, G or T, then the end. Decide which of the four strings it matches before you run the code.

%%sandbox
import re

for seq in ['ACGTTGCA', 'ACGTNGCA', 'acgt', '']:
    print(repr(seq), re.search('^[ACGT]+$', seq))

Note that ^ means two different things: outside square brackets it is start of the string, and as the first character inside them it is not. Same symbol, different job.

Exercise 59-12

SOLO

Write a function is_dna that takes a string and returns True if it consists only of the uppercase letters A, C, G and T, and False otherwise. A function that returns True or False should return exactly those values and not a match object or None, so you need one small extra step after the search. These lines should all print True:

print(is_dna('ACGTTGCA') == True)
print(is_dna('ACGTNGCA') == False)
print(is_dna('acgt') == False)
%%sandbox
import re

# Your code here

Finding all of them

re.search stops at the first match. To get every match, use re.findall, which returns a list of the matched strings:

%%sandbox
import re

dna = 'TTGACTCAAGATTCCGAGTCA'
re.findall('GA.TC', dna)

That is three HinfI sites, and you can see which base sits in the middle of each. A list of strings does not tell you where they are, though. For that, use re.finditer. It gives you the match objects one at a time, so it works with a for-loop just like a list or a file does:

%%sandbox
import re

dna = 'TTGACTCAAGATTCCGAGTCA'
for match in re.finditer('GA.TC', dna):
    print(match.group(), 'at', match.start())

Exercise 59-13

SOLO

How many times do you think ATA occurs in ATATATA? Count them yourself by sliding your finger along the string. Then predict what re.findall returns, and run the code.

%%sandbox
import re

re.findall('ATA', 'ATATATA')

You should have counted three, at index 0, 2 and 4. re.findall found two. Once the regular expression machine has found a match, it continues searching right after that match, so matches never overlap. The .count method on strings works the same way. For restriction sites that rarely matters, but for short repeated motifs it can quietly give you the wrong count.

Exercise 59-14

AI: Explainer

Ask the assistant how to make re.findall (or re.finditer) also find overlapping matches. It will most likely give you a pattern containing something you have not seen, such as (?=...). Do not accept it until you can say in your own words what each character in the pattern does, so keep asking about the parts you cannot explain. Then let the machine settle it: run the suggested code on 'ATATATA' and on 'AAAAA' with the pattern AA, and compare both results with the count you get by sliding your finger along the string. Write in your logbook whether the assistant’s explanation agreed with what the code actually did.

Groups: pulling out the pieces

Often you do not just want to know that a pattern matched; you want to pull pieces out of the text it matched. Wrap a piece of the pattern in parentheses, and the match object remembers what that piece matched. The parenthesized pieces are called groups, and they are numbered from 1 in the order of their opening parentheses.

Genomic coordinates are often written like chr7:114086327-114693772: the chromosome, a colon, the start position, a dash and the end position. Here is a pattern that takes it apart:

%%sandbox
import re

region = 'chr7:114086327-114693772'
match = re.search(r'(\w+):(\d+)-(\d+)', region)
print(match.group(1))
print(match.group(2))
print(match.group(3))

match.group() with no argument (or with 0) still gives you the whole match. match.groups() gives you all the groups at once as a tuple.

Exercise 59-15

SOLO

Using region and the pattern from above, compute the length of the region by subtracting the start from the end. Before you run your code, decide what type match.group(2) returns. If your first attempt gives you an error, read it, and remember what you know about converting between types.

%%sandbox
import re

region = 'chr7:114086327-114693772'
match = re.search(r'(\w+):(\d+)-(\d+)', region)

# Your code here

When the pattern has groups, re.findall changes what it returns. Instead of a list of whole matches, you get a list of tuples with one string per group. That fits nicely with the tuple unpacking you learned in Chapter 22:

%%sandbox
import re

text = 'Regions: chr7:114086327-114693772, chrX:1000-2000 and chr2:500-900.'
for chrom, start, end in re.findall(r'(\w+):(\d+)-(\d+)', text):
    print(chrom, int(end) - int(start))

Exercise 59-16

SOLO

This is the DEFINITION line from the GenBank file for a FOXP2 transcript:

Homo sapiens forkhead box P2 (FOXP2), transcript variant 4, mRNA.

Write a pattern with one group that pulls out the gene symbol between the parentheses, so that match.group(1) is 'FOXP2'. Parentheses are metacharacters, so the ones you want to match need escaping (see ?sec-re-escaping) while the one you want as a group does not. Then add a second group that pulls out the variant number.

%%sandbox
import re

definition = 'Homo sapiens forkhead box P2 (FOXP2), transcript variant 4, mRNA.'

# Your code here

Replacing and splitting

re.sub finds every match of a pattern and replaces it with something else. It takes three arguments: the pattern, the replacement and the string. This is very handy for cleaning up data.

Here are the first three lines of the sequence part of that same GenBank file, exactly as they appear in it. Each line starts with a number, and the bases come in blocks of ten, separated by spaces:

%%sandbox
import re

block = '''
        1 acagacatga aagctaaccg aggacttgag agactcaaac tggtgctttt gtctctctct
       61 ctctgtcttt ctctctctca cacacacact cacacactca cacacatgca cacacacaca
      121 tacacacaca caaaaatgaa gcacttactt tagaaagatt atggtaagca tgctggctca
'''
seq = re.sub(r'[\d\s]', '', block).upper()
print(seq)
print(len(seq))

The pattern [\d\s] matches any single digit or whitespace character, and each one is replaced by the empty string, i.e., deleted. Then .upper() makes it look like the DNA you are used to.

Exercise 59-17

SOLO

Look at the second line of the block. It contains a long stretch of CA repeated over and over. Short repeats like this are called microsatellites, and they are used for, among other things, DNA fingerprinting. A quantifier can apply to a whole group: (CA)+ means one or more copies of CA. Starting from the cleaned seq, use re.finditer with the pattern (CA){5,} to print every run of at least five CA repeats, with its start position and its length. Predict how many runs you will find by looking at the block before you run your code.

Why did you have to clean the sequence first? What would re.finditer find in the raw block?

Finally, try re.findall with the same pattern, and predict what it returns before you run it. Remember what re.findall does when the pattern contains a group (?sec-re-groups).

%%sandbox
import re

block = '''
        1 acagacatga aagctaaccg aggacttgag agactcaaac tggtgctttt gtctctctct
       61 ctctgtcttt ctctctctca cacacacact cacacactca cacacatgca cacacacaca
      121 tacacacaca caaaaatgaa gcacttactt tagaaagatt atggtaagca tgctggctca
'''
seq = re.sub(r'[\d\s]', '', block).upper()

# Your code here

re.split is the regular-expression version of the string method .split. Where .split cuts the string at one fixed separator, re.split cuts it wherever the pattern matches. That helps when a colleague’s gene list was typed by hand:

%%sandbox
import re

genes = 'BRCA1, TP53;FOXP2   MYC,KRAS'
print(genes.split(','))
print(re.split(r'[,;\s]+', genes))

Peek behind the curtain: compiled patterns

When you call re.search, the re module first translates your pattern into an internal form it can run quickly. You can do that translation step yourself with re.compile, and you get a pattern object back. Like everything else in Python, it is an object with methods:

%%sandbox
import re

hinfi = re.compile('GA.TC')
print(type(hinfi))
print(hinfi.findall('TTGACTCAAGATTCCGAGTCA'))

re.findall('GA.TC', dna) says “Hey function, find this pattern in this string!”. hinfi.findall(dna) says “Hey pattern, find yourself in this string!”. It is the same function-versus-method story as capitalize. They do exactly the same thing. A compiled pattern is mostly a convenience: you give it a good name once and reuse it in many places. And print(type(...)) shows you that the match objects you have been using all along are objects of the class re.Match.

When not to use them

There is an old programmer joke: “Some people, when confronted with a problem, think ‘I know, I’ll use regular expressions.’ Now they have two problems.” Regular expressions are dense. A pattern that took you ten minutes to write can take you twenty to read a month later, and every character in it matters. So here is a rule of thumb. If you are looking for exact text, use in, .find, .count, .startswith or .split. They say what they do, and nobody needs to decode them. Reach for re when the thing you are looking for is a pattern: alternative letters at a position, repeats of unknown length, or pieces you want to pull out of a line.

Tip

Regular expressions could shortcut parts of the upcoming ORF project. Solve the project with the loops and slicing it asks for. The project is about decomposing a problem into functions, and a one-line pattern that you cannot debug teaches you nothing about that.

Regular expressions and the assistant

Regular expressions are one of the things AI assistants are best at writing and one of the things people are worst at reading. That is an uncomfortable combination. A pattern looks the same whether it is right or wrong, so the only way to know is to test it on strings for which you decided the answer beforehand.

Exercise 59-18

AI: Explainer

Here is a pattern of the kind you will find in other people’s code for reading FASTA files:

r'^>(\S+)\s+(.*)$'

Give the assistant this pattern and nothing else, and ask it to translate it into English one piece at a time. Then, with the assistant closed, write down three header lines you predict it matches and two you predict it does not, and for the ones that match, what you expect in group 1 and group 2. Run them through re.search and compare. Your predictions settle it, not the explanation: if a result surprises you, find the piece of the translation that misled you or that you misread.

%%sandbox
import re

pattern = r'^>(\S+)\s+(.*)$'
for header in ['>NM_148900.4 Homo sapiens forkhead box P2',
               # your test headers here
              ]:
    match = re.search(pattern, header)
    if match:
        print(repr(header), '->', match.groups())
    else:
        print(repr(header), '-> no match')

Exercise 59-19

AI: Comparer

Ask the assistant for two versions of the is_dna function from ?sec-re-anchors: one that uses re and one that uses no regular expression at all. Before you run either of them, decide what each one returns for '', the empty string, and whether the two versions agree. Then run both on 'ACGT', 'ACGN', 'acgt' and ''. Many reasonable versions disagree on the empty string, and neither is wrong until you decide whether an empty string is DNA. Write down which version you would rather read in six months, and what you think the empty string should give.

Exercise 59-20

AI: Drafter

Restriction sites are usually written in IUPAC code, where N is any base, W is A or T, S is G or C, R is A or G and Y is C or T. You want a function find_sites(seq, site) that takes a DNA sequence and an IUPAC site such as 'GANTC', and returns a list of the start positions of every match. To do that, it must translate the site into a regular expression, 'GANTC' into 'GA[ACGT]TC' and so on.

Do it in the Drafter order. First, with the assistant closed, finish the tests below: work out by hand the right answer for each one and add at least two more, including one with W or S and one where there is no match at all. Then give the assistant the name, the arguments, the return value and the IUPAC table, but not your tests, and ask for a draft. Read every line, and ask about any line you cannot explain. Paste the draft above the tests and run the cell. The tests decide whether the draft is right, not how convincing it looks.

%%sandbox
import re

# Paste the draft here


assert find_sites('GAATCGACTC', 'GANTC') == [0, 5]
assert find_sites('TTGAATTCAA', 'GAATTC') == [2]
# Your tests here

print('All tests passed')

General exercises

Exercise 59-21

SOLO

What does re.findall('[ACGT]{3}', 'ATGGCCTAA') return? Decide, then run it. Can you think of a use for it?

%%sandbox
import re

re.findall('[ACGT]{3}', 'ATGGCCTAA')

Exercise 59-22

SOLO

Write a function is_protein that returns True if a string consists only of the one-letter codes for the twenty standard amino acids (ACDEFGHIKLMNPQRSTVWY) and False otherwise. Which of these do you expect to be proteins: 'MKTAYIAK', 'MKTXYIAK', 'ACGT'? Decide before you test your function on them.

%%sandbox
import re

# Your code here

Exercise 59-23

SOLO

GC-rich stretches can make PCR primers stick in the wrong places. Write a function gc_runs that takes a DNA sequence and returns a list of every run of at least four bases that are all G or C. For 'ATGCGCCATTAGGGCATAGCA' it should return ['GCGCC', 'GGGC'].

%%sandbox
import re

# Your code here

Exercise 59-24

SOLO

Write a function parse_headers that takes a list of FASTA header lines like the ones below and returns a dictionary that maps each accession (the text between > and the first space) to the rest of the line.

headers = ['>NM_148900.4 Homo sapiens forkhead box P2 (FOXP2), transcript variant 4, mRNA',
           '>NM_000546.6 Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA',
           '>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome']

For the headers above, parse_headers(headers)['NM_000546.6'] should be 'Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA'.

%%sandbox
import re

# Your code here

Exercise 59-25

SOLO

Change your parse_headers so that the dictionary maps each accession to the gene symbol in parentheses instead, and to None if the header has no gene symbol. Which of the three headers above gives None? Decide before you run it.

%%sandbox
import re

# Your code here