Lecture 7: Sequence Alignment¶
CBIO (CSCI) 4835/6835: Introduction to Computational Biology
Overview and Objectives¶
In our last lecture, we covered the basics of molecular biology and the role of sequence analysis. In this lecture, we'll dive deeper into how sequence analysis is performed and the role of algorithms in addressing sequence analysis. By the end of this lecture, you should be able to:
- Define the notion of algorithmic complexity and how it relates to sequence alignment and analysis
- Describe and define the abstract problems of shortest common superstring (SCS) and longest common substring (LCS), and how they specifically relate to sequence analysis
- Recall different methods of scoring sequence alignments and their advantages and drawbacks
- Describe the different distance metrics and methods of scoring sequence alignments
- Explain why local or global sequence alignments are preferred in certain situations
Part 0: range¶
This mysterious range function has showed up a few times so far. What does it do?
r = range(10)
print(r)
Not terribly useful output information, to be fair.
r = range(10)
l = list(r) # cast it as a list
print(l)
range(i) generates a list of numbers from 0 (inclusive) to i (exclusive). This is very useful for looping!
for i in range(10):
print(i, end = " ")
You can also provide a second argument to range, which specifies a starting point for the counting (other than 0). That starting point is still inclusive, and the ending point still exclusive.
for i in range(5, 10):
print(i, end = " ")
Finally, you can also provide a third argument, which specifies the interval between numbers in the output. So far, that interval has been 1: start at 0, to go i, by ones. You can change that "by ones" to whatever you want.
for i in range(5, 10, 2): # read: from 5, to 10, by 2
print(i)
You can get really crazy with this third one, if you want: you can go backwards by putting in a negative interval.
for i in range(10, 0, -2): # from 10, to 0, by -2
print(i)
Same rules apply, though: the starting point is inclusive (hence why we see a 10), and the ending point is exclusive (hence why we don't see a 0).
range is particularly useful as a way of looping through a list of items by index.
list_of_interesting_things = [93, 17, 5583, 47, 2359875, 4, 381]
for item in list_of_interesting_things:
print(item, end = " ")
This is how we've seen loops so far: the loop variable (here it's item) is a literal item in the list. But what if, in addition to the item, I needed to know where in the list that item was (i.e., the item's list index)?
list_length = len(list_of_interesting_things)
for index in range(list_length): # use range of the list length!
item = list_of_interesting_things[index] # pull out the item AT that index
print("Item " + str(item) + " at index " + str(index))
Part 1: Complexity¶
Big "Oh" Notation¶
From computer science comes this notion: how the runtime of an algorithm changes with respect to its input size.
$\mathcal{O}(n)$ - the "$\mathcal{O}$" is short for "order of the function", and the value inside the parentheses is always with respect to $n$, interpreted to be the variable representing the size of the input data.
Limits¶
Big-oh notation is a representation of limits, and most often we are interested in "worst-case" runtime. Let's start with the example from the last lecture.
a = [1, 2, 3, 4, 5]
for element in a:
print(element)
How many steps, or iterations, does this loop require to run?
Alright, back to complexity:
a = range(100)
for element in a:
print(element, end = " ")
How many iterations does this loop require?
For iterating once over any list using a single for loop, how many iterations does this require?
Algorithms which take $n$ iterations to run, where $n$ is the number of elements in our data set, are referred to as running in $\mathcal{O}(n)$ time.
This is roughly interpreted to mean that, for $n$ data points, $n$ processing steps are required.
Important to note: we never actually specify how much time a single processing step is. It could be a femtosecond, or an hour. Ultimately, it doesn't matter. What does matter when something is $\mathcal{O}(n)$ is that, if we add one more data point ($n + 1$), then however long a single processing step is, the algorithm should take only that much longer to run.
How about this code? What is its big-oh?
a = range(100)
b = range(1, 101)
for i in a:
print(a[i] * b[i], end = " ")
Still $\mathcal{O}(n)$. The important part is not (directly) the number of lists, but rather how we operate on them: again, we're using only 1 for loop, so our runtime is directly proportional to how long the lists are.
How about this code?
a = range(100)
x = []
for i in a:
x.append(i ** 2)
for j in a:
x.append(j ** 2)
Trick question! One loop, as we've seen, is $\mathcal{O}(n)$. Now we've written a second loop that is also $\mathcal{O}(n)$, so literally speaking the runtime is $2*\mathcal{O}(n)$, but what happens to the 2 in the limit as $n \rightarrow \infty$?
The 2 is insignificant, so the overall big-oh for this code is still $\mathcal{O}(n)$.
How about this code?
a = range(100)
for element_i in a:
for element_j in a:
print(element_i * element_j, end = " ")
Nested for loops are brutal--the inner loop runs in its entirety for every single iteration of the outer loop. In the limit, for a list of length $n$, there are $\mathcal{O}(n^2)$ iterations.
One more tricky one:
xeno = 100
while xeno > 1:
xeno /= 2
print(xeno, end = " ")
Maybe another example from the same complexity class:
xeno = 100000
while xeno > 1:
xeno /= 10
print(xeno, end = " ")
What does this "look" like?
# I'm just plotting the iteration number against the value of "xeno".
%matplotlib inline
import matplotlib.pyplot as plt
x = []
y = []
xeno = 10000
i = 1
while xeno > 1:
x.append(i)
y.append(xeno)
xeno /= 10
i += 1
plt.plot(x, y)
In the first one, on each iteration, we're dividing the remaining space by 2, halving again and again and again.
In the second one, on each iteration, we're dividing the space by 10.
$\mathcal{O}(\log n)$. We use the default (base 10) because, in the limit, constants don't matter.
Part 2: SCS and LCS¶
Recall from the last lecture what SCS (shortest common superstring) was:
- The shortest common superstring, given sequences $X$ and $Y$, is the shortest possible sequence that contains all the sequences $X$ and $Y$.
For example, let's say we have $X$ = ABACBDCAB and $Y$ = BDCABA. What would be the shortest common superstring?
Here is one alignment: BDCABA (second string) and ABACBDCAB (first string). The ABA is where the two strings overlap. The full alignment, BDCABACBDCAB, has a length of 12.
Can we do better?
ABACBDCAB and BDCABA, which gives a full alignment of ABACBDCABA, which has a length of only 10. So this alignment would be the SCS.
(When do we need to use SCS?)
Longest Common Substring (LCS)¶
In a related, but different, problem: longest common substring asks:
- Given sequences $X$ and $Y$, the longest common substring is the constituent of the sequences $X$ and $Y$ that is as long as possible.
Let's go back to our sequences from before: $X$ = ABACBDCAB and $Y$ = BDCABA. What would be the longest common substring?
The easiest substrings are the single characters A, B, C, and D, which both $X$ and $Y$ have. But these are short: only length 1 for all. Can we do better?
ABACBDCAB and BDCABA, so the longest common substring is BDCAB.
(When do we need LCS?)
Rudimentary Sequence Alignment¶
Given two DNA sequences $v$ and $w$:
$v$: ATATATAT
$w$: TATATATA
How would you suggest aligning these sequences to determine their similarity?
Before we try to align them, we need some objective measure of what a "good" alignment is!
Part 3: Distance Metrics¶
Hopefully, everyone has heard of Euclidean distance: this is the usual "distance" formula you use when trying to find out how far apart two points are in 2D space.
How is it computed?
For two points in 2D space, $a$ and $b$, their Euclidean distance $d_e(a, b)$ is defined as:
$d_e(a, b) = \sqrt{(a_x - b_x)^2 + (a_y - b_y)^2)}$
So if $a = (1, 2)$ and $b = (5, 3)$, then:
$d_e(a, b) = \sqrt{(1 - 5)^2 + (2 - 3)^2} = \sqrt{(-4)^2 + (-1)^2} = \sqrt{16 + 1} = 4.1231$
How can we measure distance between two sequences?
There is a metric called Hamming Distance, which counts number of differing corresponding elements in two strings.
We'll represent the Hamming distance between two strings $v$ and $w$ as $d_H(v, w)$.
$v$: ATATATAT
$w$: TATATATA
$d_H(v, w)$ = 8
That seems reasonable. But, given how similar the two sequences are (after all, the LCS of these two is 7 characters), what if we shifted one of the sequences over by one space?
$v$: ATATATAT-
$w$: -TATATATA
Now, what's $d_H(v, w)$?
$d_H(v, w)$ = 2
The only elements of the two strings that don't overlap are the first and last; they match perfectly otherwise!
Edit distance¶
Hamming distance is useful, but it neglects the possibility of insertions and deletions in DNA (what is the only thing it counts?). So we need something more robust.
The edit distance between two strings is the minimum number of elementary operations (insertions, deletions, or substitutions / mutations) required to transform one string into the other.
Hamming distance: $i^{th}$ letter of $v$ with $i^{th}$ letter of $w$ (how hard is this to do?)
Edit distance: $i^{th}$ letter of $v$ with $j^{th}$ letter of $w$ (how hard is this to do?)
Hamming distance is easy, but gives us the wrong answers. Edit distance gives us much better answers, but it's hard to compute: how do we know which $i$ to pair with which $j$?
What's the edit distance for $v$ = TGCATAT and $w$ = ATCCGAT?
One solution:
- TGCATAT (delete last T)
- TGCATA (delete last A)
- ATGCAT (insert A at front)
- ATCCAT (mutate G to C)
- ATCCGAT (insert G before last A)
ATCCGAT == ATCCGAT, done in 5 steps!
Can it be done in 4 steps?
(...mmmmaybe--but that's for next week!)
Part 4: Global vs Local Alignment¶
indel is a portmanteau of "insertion" and "deletion", so we don't need to worry about which strand we're actually referring to.
What is the edit distance here?
Highly conserved subsequences¶
Things get hairier when we consider that two genes in different species may be similar over short, conserved regions and dissimilar over remaining regions.
Homeobox regions have a short region called the homeodomain that is highly conserved among species--responsible for regulation of patterns of anatomical development in animals, fungi, and plants.
- A global alignment would not find the homeodomain because it would try to align the entire sequence.
- Therefore, we search for an alignment which has a low edit score locally, meaning we have to search aligned substrings of the two sequences.
Here's an example global alignment that minimizes edit distance over the entirety of these two sequences:
Here's an example local alignment that may have an overall larger edit distance, but it finds the highly conserved substring:
"BUT!", you protest.
"If the local alignment has a higher edit score, how do we find it at all?"
We've already seen that we need to consider three separate possibilities when aligning sequences:
- Insertions / Deletions (characters added or removed)
- Mutations / Substitutions (characters modified)
- Matches (characters align)
With Hamming distance, two characters were either the same or they weren't (options 1 and 2 above were a single criterion).
With edit distance, we separated #1 and #2 above into their own categories, but they are still weighted the same (1 insertion = 1 mutation = 1 edit)
Are all insertions / deletions created equal? How about all substitutions?
Scoring Matrices¶
Say we want to align the sequences:
$v$ = AGTCA
$w$ = CGTTGG
But instead of using a standard edit distance as before, I give you the following scoring matrix:
This matrix gives the specific edit penalties for particular substitutions / insertions / deletions.
It also allows us to codify our understanding of biology and biochemistry into how we define a "good" alignment. For instance, this penalizes matching A with G more heavily than C matched with T.
Here is a sample alignment using this scoring matrix:
Making a scoring matrix¶
Scoring matrices are created based on biological evidence.
Some mutations, especially in amino acid sequences, may have little (if any!) effect on the protein's function. Using scoring matrices, we can directly quantify that understanding.
- Polar to polar mutations (aspartate -> glutamate)
- Nonpolar to nonpolar mutations (alanine -> valine)
- Similarly behaving residues (leucine -> isoleucine)
Standard scoring matrices¶
For nucleotide sequences, there aren't really "standard" scoring matrices, since DNA is less conserved overall and less effective to compare coding regions.
There are, however, some common amino acid scoring matrices. We'll discuss two:
- PAM (Point Accepted Mutation)
- BLOSUM (Blocks Substitution Matrix)
PAM¶
PAM is a more theoretical model of amino acid substitutions.
It is always associated with a number, e.g. 1 PAM, written as PAM$_1$. This means the given PAM$_1$ scoring matrix is built to reflect a 1% average change in all amino acid positions of the polypeptide, according to evolution.
Some important notes:
- This is an average. Even with PAM$_{100}$, not every residue will have changed (some are more conserved than others)
- Some residues may have mutated several times!
- Some residues may have mutated back to their original state!
- Some residues may not have changed at all
PAM$_{250}$ is a widely used scoring matrix.
Mutating A to A is clearly the most preferable (highest score in that row of 13 points), but after 250 evolutions, a mutation from A to G also seems very favorable (12 points).
BLOSUM¶
Unlike PAM, scores in BLOSUM are derived from direct empirical observations of the frequencies of substitutions in blocks of local alignments in related proteins. They both, however, often obtain identical alignment scores.
Like PAM, BLOSUM also has a number associated with it, this time to represent the observed substitution rate between two proteins sharing some amount of similarity.
BLOSUM$_{62}$ is a common scoring matrix, representing substitution rates in proteins sharing no more than 62% identity.
Next week¶
We'll look at how to use these matrices to determine the best alignments of sequences!
Administrivia¶
- You should everything you need now for Assignment 2. Any questions so far? Due next week!
Additional Resources¶
- Jones, Neil C. and Pevzner, Pavel A. An Introduction to Bioinformatics Algorithms, Chapter 6. 2004. ISBN-13: 978-0262101066
- Based heavily on the modified slides of Dr. Phillip Compeau.