Finding sequences

Three ways to search a genome

The question for today

You have a long sequence and a short one.

Where does the short one occur in the long one?

  • Where would EcoRI cut this plasmid?
  • Where does this 20-nucleotide sequencing read come from?
  • Does this promoter motif appear upstream of my gene?
  • Which of these 10 million reads contain the adapter?

All of these are the same computational question

Why this is a computational problem

Sequence lengths

phiX174 5.4 kb
E. coli 4.6 Mb
Human 3.1 Gb

Query volumes

One Illumina run can produce billions of reads.

Each one needs to be placed somewhere in the genome.

An approach that takes one second per query is fine for a few queries.

But, that is 32 years for a billion queries.

Stating the problem

Given a text \(T\) of length \(n\) and a pattern \(P\) of length \(m\):

Report every position in \(T\) where \(P\) occurs.

Two things to notice:

  • We want all occurrences, not just the first one
  • \(m\) is usually much smaller than \(n\)

In the exercise, you will build three solutions and time them against each other

Approach 1

What does a linear search cost?

  • Setup: none
  • Memory: essentially none
  • Per search: up to \(n\) window positions, each costing up to \(m\) character comparisons, so the time complexity is \(O(nm)\)

For DNA the comparison usually does not match on the first or second character, so in practice it behaves like \(O(n)\).

The central idea

Do the work once.

Instead of re-reading the genome for every query, read it one time and build a structure that makes future queries fast.

Buying time with memory!

This is called indexing, and it’s analogous to the index at the back of a textbook.

Approach 2

A \(k\)-mer dictionary

A \(k\)-mer is any substring of length \(k\).

Walk over the sequence once, and record where every \(k\)-mer starts.

Index:    0 1 2 3 4 5 6 7 8 9 10
Sequence: A T C G T T C A T C G
{
  "ATC": [0, 7],   "TCG": [1, 8],   "CGT": [2],
  "GTT": [3],      "TTC": [4],      "TCA": [5],
  "CAT": [6],
}

Looking things up is very fast

Searching for "TCG" is now a dictionary lookup, not a search!

kmer_dict["TCG"]   # -> [1, 8]
  • Setup: one pass over the sequence, time complexity \(O(nk)\)
  • Per search: \(O(1)\) on average, regardless of how long the genome is
  • Memory: a key and a position for nearly every position in the sequence

The genome could be a thousand times longer and the lookup would take the same time. But, it will take up more memory!

Two catches

  1. The length is fixed. You can only search for 3-mers in a 3-mer dictionary. Searching for "ATCAG" means finding every "ATC" and checking whether "AG" follows. This is called \(k\)-mer expansion.

  2. It costs memory. Some implementations reserve a slot for every possible \(k\)-mer, and there are \(4^k\) of those. E.g., about \(1.1 \times 10^{12}\) for \(k = 20\).

    You can save memory by only storing the \(k\)-mers that actually occur in the sequence, but this is still several times the size of the sequence.

Approach 3

Suffixes

A suffix is any substring that runs to the end of the sequence.

ATTCGTTCAG has ten of them:

ATTCGTTCAG  0     TTCAG  5
 TTCGTTCAG  1      TCAG  6
  TCGTTCAG  2       CAG  7
   CGTTCAG  3        AG  8
    GTTCAG  4         G  9

Note: Every substring of a sequence is the beginning of some suffix.

Sort the suffixes

row pos suffix
0 8 AG
1 0 ATTCGTTCAG
2 7 CAG
3 3 CGTTCAG
4 9 G
5 4 GTTCAG
6 6 TCAG
7 2 TCGTTCAG
8 5 TTCAG
9 1 TTCGTTCAG

Sorted alphabetically

Look at the two suffixes starting with TC

  • They are adjacent. Sorting guarantees this.
  • Nothing that doesn’t start with TC can sort in between them

Suffix array

Storing all those suffix strings is costly (roughly \(n^2/2\) characters)

But we don’t need them. Given the sequence and a position, any suffix can be recovered with a slice.

So a suffix array is just the positions:

[8, 0, 7, 3, 9, 4, 6, 2, 5, 1]

\(n\) integers instead of \(n^2/2\) characters

Searching a suffix array

Because the suffixes are sorted, we can use binary search

  • Compare the target to the middle entry
  • Throw away the half it cannot be in
  • Repeat

Each step halves what is left, so the cost (time complexity) grows as \(\log_2 n\)

For a 5,000-nucleotide sequence, that’s 13 comparisons instead of 5,000.

Example: Finding TC

row pos suffix
0 8 AG
1 0 ATTCGTTCAG
2 7 CAG
3 3 CGTTCAG
4 9 G
5 4 GTTCAG
6 6 TCAG
7 2 TCGTTCAG
8 5 TTCAG
9 1 TTCGTTCAG
step window middle row suffix there verdict
1 rows 0–9 4 G G < TC, go down
2 rows 5–9 7 TCGTTCAG match

Only 2 comparisons. But we have only found one of the two matches.

Binary search is only part of the solution

It finds a match and stops. We threw away half the array at every step, so other matches can be discarded.

Searching a suffix array in 2 stages

Remember that matches are adjacent

  1. Binary search to land on any one matching row
  2. Walk outward — up, then down — until the neighbors stop matching
row 5   GTTCAG      no match  <- stop walking up
row 6   TCAG        match     <- position 6
row 7   TCGTTCAG    match     <- position 2  (found by binary search)
row 8   TTCAG       no match  <- stop walking down

Answer: positions 2 and 6. Cost: about \(O(m \log n)\).

Putting it all together

Three approaches

Setup Per search Memory Any query length?
Linear search none \(O(nm)\) none yes
\(k\)-mer dictionary \(O(nk)\) \(O(1)\) large no — fixed \(k\)
Suffix array sorting \(O(m \log n)\) \(n\) integers yes

There is no “best algorithm” in this table.

There is only a “best algorithm” for a given problem.

The indexing strategy

Build once. Search many.

For one search, linear search wins. Both the \(k\)-mer index and suffix array have to read the whole genome just to get started.

For many searches, the setup cost is recovered with every query. If doing many searches, the payoff can be huge.

This is what real tools do

  • BLAST indexes short words, finds exact seed matches, then extends them
  • Bowtie, BWA build an FM-index (based on the Burrows-Wheeler transform), a compressed cousin of the suffix array
  • minimap2 indexes a cleverly chosen subset of \(k\)-mers called minimizers

They all use the same strategy: A genome gets searched far more often than it changes, so pay the indexing cost up front.

In the exercise

You will build all three using Python functions

  1. linear_search
  2. get_kmer_dict and kmer_dict_search
  3. get_suffix_array and suffix_binary_search

Use all three to search the phiX174 genome for the recognition sites of restriction enzymes

Watch for

The timing result might surprise you. Try to work out why before you read the explanation.

Exercise

  1. Log into molab
  2. Open the exercise in molab
  3. Click “Fork” button

Acknowledgements

I used Anthropic’s Claude (Opus 5 model) to generate a draft of these slides.