Skip to main content
Loading Explore guides…
Biological Sequence Analysis · Acedra
Home
Explore
Biology
Biological Sequence Analysis
Sign in to save
Chat
If YouTube asks you to sign in, open the video on YouTube.
Watch on YouTube
Reload player
Transcript
143 segments
Copy
Hide video
0:04
Andy Baxevanis: Okay. So, today’s lecture is devoted to biological sequence analysis, and in the same way that Eric Green next week is going to lay the groundwork for a lot of the lectures that follow in this course by providing you an overview of the field of genomics, my job is to provide you a solid foundation for the upcoming lectures devoted to bioinformatics. As I think is obvious to most people at this point, genomics and bioinformatics go hand
0:30
in hand with one another, and this sentiment is best illustrated by this figure from NHGRI’s current strategic plan showing the spectrum of genomic-based research going from understanding the structure of genomes, understanding the biology that is encoded by those genomes, moving onto disease, advancing the science of medicine, and improving the effectiveness of health care, and our ability to work across this entire spectrum really depends in large
0:57
part on our ability to use computational approaches to analyze the data that are being generated in the course of whole genome sequencing efforts, genome-wide association studies, clinical sequencing effort, pharmacogenomics studies, or any of the many other kinds of genomic science that you’ll be hearing about throughout this series. So, given this inextricable relationship between genomics and bioinformatics, my job, again,
1:21
during the two weeks I’ll be lecturing is to give you as gentle of an introduction as I can to some of the major tools that are used in bioinformatics. You’re going to see these tool being applied over and over again in many of the lectures in this series, and the concepts that are underlying these tools are very important to understand, so having that basic understanding will take these techniques, many of them that
1:43
you may have used before without quite understanding how they actually work, and take them away from the realm of being the black box. So, I hope you’ll find this background useful and pick -- that you’ll pick up a few strategies along the way that will hopefully be of practical use in your own studies. So, we’re going to start off with a discussion of sequence analysis methods by reinforcing the importance of alignments as one of the most powerful approaches we have in the study
2:09
of sequence data, specifically with respect to determining similarity and deducing homology. Those are words that carry special meaning and provide the underpinning for much of the fields of comparative and illusionary genomics. So, why construct sequence alignments? And you could ask this question. It applies to whether we’re talking about a pair-wise alignment, comparing two sequences to one another or whether we’re looking at a multiple sequence alignment comparing
2:35
a collection of presumably related sequences. So, these alignments allow us to come up with some measure of relatedness between nucleotide or, you know, acid sequences. It allows us to draw biological inferences regarding structural, functional, or evolutionary relationships, and because of that, it’s really important to use the correct terminology when we talk about these phylogenetic relationships. So, let’s go back to those terms that I just mentioned a moment ago.
3:02
I want you to imagine that you have two sequences, A and B, and that you align them, and the easiest way to convey how related these sequences are is to simply count the number of matches you -- and mismatches you have between those two sequences, and when you do that, you have a quantitative measure called similarity, so similarity is always based on an observable -- again, in this case, how many residues are lining up with each other.
3:26
It’s usually expressed as a percent identity, and what it allows you to do is quantify changes that occur as the two sequences in question diverge, so looking for substitutions, insertions, or deletions. It also helps to very quickly narrow down which of the residues that have been aligned are important for maintaining a protein’s structure or function, so whenever we see a high degree of similarity—again, this observable metric—it will imply either a
3:55
common evolutionary history or possible commonality in biological function. Okay, so regions that are important from a functional standpoint when thinking about things at the protein level, regions that might be important at the nucleotide level, which are involved in things like expression. So, the second term we need to define is the term that is more often than not misused, and that term is homology. So, homology implies an evolutionary relationship.
4:25
It -- so -- and when we have two proteins that are called homologues, that means that their genes have arisen from a common ancestor. Homologues -- two genes or two proteins either are or are not homologous, okay? This is not measured in degrees. This is the conclusion that one draws based on a similarity metric, and there’s a very good quote that I always like to show that drives this point home. It’s -- and this comes from Walter Fitch, one of the grandfathers of the field of bioinformatics.
4:56
It’s worth repeating here that “homology, like pregnancy, is indivisible.” You are either homologous—pregnant—or you are not, and this I -- quote I show because you will not forget this at this point that this is the conclusion. You never say 80 percent homology. You might say 80 percent identity when comparing two sequences to one another, but the conclusion based on some degree of identity is homology.
5:19
Now, just to make things even more complicated, there are several -- the term homologue might actually describe several different types of evolutionary relationships, and we’re going to focus specifically on two of them here. So, the first one is called orthologues. These are genes that have diverged as a result of a speciation event, so put more simply, it’s just -- you can just think of orthologues as the same gene in different species.
5:46
Because of that, they’re direct descendants of the sequence in a common ancestor. They most likely have similar domain and three-dimensional structure, usually retain the same biological function over time, and can be used to predict gene function in novel genomes. The other term, the other is called a paralogue, and paralogues are genes that arose by the duplication of a single gene in a particular lineage, so here you can simply think of paralogues
6:16
as similar genes in the same species that may have perhaps coopted a preexisting gene for new biological functions, and we see that all the time over evolutionary time, thus the whole concept of evolutionary innovation. Okay, so trying to put this all together in a diagram because I realize that these terms are confusion -- confusing. Let’s start with here, gene A, that is in some common ancestor, and over time there
6:44
were a series of speciation events that gave rise to these three genes up here, one, two, and three, so since they’re all direct descendants of this gene Alpha, one, two, and three are orthologues of one another, okay? All right, now if we roll the clock back a little bit and consider what happens when you have a gene duplication event, so we’ll go back in time prior to Alpha, some event happened here, giving rise to both Alpha and Beta.
7:13
Alpha -- and perhaps it’s an exact copy; perhaps it’s slightly different, and in the same way that this lineage on the left has evolved over time, this lineage on the right has evolved over time, somewhat differently in this case. So, the relationships here -- again, one, two, and three are orthologous because they all arose from the common ancestor A. Four, five, and six are orthologous to one another because they all arose from ancestor B, but any pair of the genes that arose from Alpha
7:44
or from Beta are paralogues of one another, so again, these are the genes that are related through a gene duplication event that you can trace back. So, again, this admittedly can be confusing. I remember when I was starting out I had a diagram like this over my desk. I would recommend you do the same, and, of course, to make things even more complicated, there are sub-definitions of each of these, so I refer you to two papers.
8:07
The first one is the Walter Fitch paper where the pregnancy quote came from. The second one is a very nice review by Eugene Koonin on orthologues, paralogues, and evolutionary genomics, so for those of you who wanted to delve a bit deeper into this subject. Okay, so that takes care of the terminology. How do we actually figure out these relationships? How do we identify these orthologues—again, the same gene in different species?
8:31
The method that’s most commonly used and straightforward is something called looking for reciprocal best BLAST hits, okay, and it doesn’t -- you don’t have to use BLAST. We’re talking about BLAST but we just want to make -- use any method that deduces similarity between two sequences. Okay, so let’s say we have a bunch of sequences from organism one and a bunch of sequences from organism two. Okay, and again, what we’re trying to do is find the genes in this genome that are
9:00
most similar to themselves in this genome, so we’ll start with the first gene up here. We do a similarity search, and it identifies this particular gene in organism two. You now do the reverse. You do the search using this particular sequence, and if it hits the same sequence in organism one, you can make the deduction that those things are -- those two genes are orthologous to one another because they both found each other in these reciprocal sequence searches,
9:31
and this is a method that works really well. It has a very low false positive rate, but it doesn’t always work, so let’s take a second case where we’ll take this particular gene, do the similarity search. It finds this gene in organism two, but then when we use this as the search query searching against the genome, it hits one of the other genes, okay, so you need to probe further here. Obviously, these genes are all related to on another.
9:57
They found each other in a sequence similarity search, but that reciprocal best relationship didn’t happen, and that is more than likely because these two genes are paralogous to one another, okay, so that -- so -- because these two have come from a relatively recent gene duplication event, it ends up confounding the analysis. It breaks the reciprocity, so to figure this out a little bit more, you have to use phylogenetic
10:26
or similar techniques to deduce exactly what’s going on in this kind of situation. Okay, so that’s it for the terminologies. Hopefully, so far, so good. All right, so let’s talk about the alignments themselves. There are two basic kinds of alignments. The first one is called global sequence alignments, where you take two sequences and you align them across their entire length. You try to align as many of the amino acid residues or nucleotides as you can across the entire sequences.
10:55
So, this works really well when you have highly similar sequences that are of relatively similar length, but as that metric declines, as they become less and less similar to one another, it really is hard to get a global alignment across two sequences. You end up forcing an alignment, and you end up missing important structural or functional domains by doing that. So, from a biological point of view, this is not going to be our method of choice.
11:24
Instead, we’re going to look at local sequence alignments, and what this does is it tries to find as many regions of similarity between the two sequences as possible, so something that are called paired subsequences, a glorified term for finding two regions within those sequences in question that are -- that can align with one another. So, any regions that are outside the area of the local alignment are excluded.
11:52
Because of the way this method works, you could actually have multiple alignments for the same two species -- two sequences that you’re trying to align. Okay, and that’s actually a good thing because it focuses on regions of high similarity, and it will pick up regions at a -- of evolutionary importance, either from a structural or a functional standpoint, so this is really the method of choice for sequences that share some similarity or for sequences of different lengths.
12:17
Okay, so this is what we’re going to focus on for most of this morning’s lecture, okay? Before we talk about the methods, a little bit about these things called scoring matrices. So, once we have these two sequences aligned, how do we actually measure how good the alignment actually is between sequence A and sequence B? Okay, so the scoring matrices come into play. They give us a metric to assess how good these alignments are.
12:45
They -- in the numbers that are found in these matrices, what they actually represent are actual biological characteristics and physical characteristics. So, physical, chemical properties of nucleotides and amino acids, so these are dependent on the sidechain structure and chemistry of the amino acids, their -- the function, as well, and it takes into account things like when you look at cysteines and prolines, they’re
13:12
obviously important for structure and functions and especially with prolines, their -- when you find them at the end of an alpha helix, they more normally end up breaking the chain. Things like tryptophan you can only put them in so many places because they have a bulky sidechain, so as far as stearate considerations go. Lysines and arginines, they both have a positive charge, overrepresented in proteins that bind
13:34
to DNA, so they can bind the negatively charged background -- backbone, but, again, taking account in -- you know, charge, shape, all of those kinds of things, so in putting these matrices together, there’s two major things we’re going to think about: conservation. What residues can substitute for another residue and not adversely affect the function of the protein okay? So, we know just empirically things like isoleucine and valine, they’re both small, they’re both hydrophilic.
14:04
They can often substitute for one another. Same thing for serines and threonines. They’re both polar, so the bottom line here, the game is to just make sure that you’re conserving charge size, hydrophobicity, and other physical chemical factors. The other thing that comes into play is frequency. When we look at the entire constellation of proteins and the amino acids in composition of those proteins, how often do we see a particular amino acid occur across all of these proteins?
14:33
So, taking into account rare versus common. Now, okay, so you know say, “All right, Andy, why should I care about any of this?” Anytime you do any sort of sequence comparisons, when you do a BLAST search, when you do a multiple sequence alignment, any kind of comparison like this, in the background, these scoring matrices are being used, and the choice of scoring matrix becomes very, very important because they implicitly represent particular evolutionary patterns, so depending on which
15:08
one you choose or not choose because you’ve decided to use the default values, it can strongly influence the outcomes of your analyses, so the default options may not always be the best choice, so that’s why we’re talking about how these are constructed, and I’m going to give you some guidelines as we go through. So, let’s take the easy case first. Let’s just look at a nucleotide scoring matrix, and, here, a simple match/mismatch scheme is employed.
15:39
So, if you look at this matrix, the four nucleotides are across the top and down the side. You’ll notice on the diagonal, okay, you see a bunch of twos, so let’s say you had an alignment and you have two Cs aligned with each other. Here’s the C, here’s the C. You see where they intersect, and for that particular position in the alignment, you’re going to assess a positive score of two points. Every place else off the diagonal, you see a whole bunch of minus threes, so anytime
16:05
there is a mismatch, you deduct three points. This assumes that each nucleotide occurs 25 percent of the time, which is actually not true, but it’s close enough, and that’s what’s usually used, and -- but this kind of scheme actually works okay because in general, the nucleotides are not all that chemically dissimilar, but when you compare them to amino acids, the situation -- the landscape becomes much more complex. So, this is an example of a protein scoring matrix.
16:38
It’s something called BLOSUM62. This is the default matrix that is used for the BLAST searches at NCBI. You’ll notice across the top, we have the amino -- all of the amino acids across the top and going down the side, so let’s take a couple of cases so you can see how this is different from what happens in the nucleotide world. Let’s say, again, you have sequence A and sequence B. You’ve aligned them, and there’s
17:02
a particular position where you have two tryptophans lining up with each other, so we’re going to look for the W up here, the W over here, and see where they intersect, and in this case, for that particular position because it matched, the value that we assess is 11 points. Now, what I want you to do, and you can probably see this better on the screen. Look at the diagonal, and you’ll see that unlike the nucleotide case where a match was always assessed the same value, here it’s not, okay?
17:31
So, here we’re taking, again, into account how frequently these amino acids occur in the entire constellation of proteins. Okay, so with tryptophan being the rarest of the residues, it had the highest score on the diagonal. Now, let’s say we -- instead of that tryptophan lining up with another tryptophan, the tryptophan lines up with a tyrosine, so a W with the Y, so now, if we look at the intersection of those two, we see a value of plus two, okay, so this takes into account the fact
18:05
that these two can and often are observed to substitute with one -- for one another in proteins, but it’s not an exact match, so we can’t give it as many points as we would with -- for the exact match. Instead, we give it a lower score, but it’s still a positive score. We take the third case where it’s just a total mismatch, so something like a cysteine for a glutamine, in this case, you see that there is a negative value here of minus three,
18:29
okay, so these are two residues that do not normally substitute for one another based on observation, so we subtract points. So, based -- the basic rules here are anytime you have an exact match, you get the most points, rare residues earning more points than an exact match of common residues. You also earn points for conservative substitutions, just not as many as you would for an exact match, and finally, you lose points for mismatches, okay?
18:59
So, that’s the way the rules work no matter what kind of -- this type of matrix that we’re looking at. So, let’s talk a little bit more about the BLOSUM matrices. The good news about these matrices is you never have to construct them yourselves. They’re all constructed for you, but again, it’s knowing which ones you have to pick. So, how was that matrix just constructed, the one we just looked at? Again, it’s looking for differences in conserved, un-gapped regions of a protein family, so
19:27
it’s taking multiple sequence alignments of known protein families and looking for often these residues can substitute for one another and how often they occur, so directly calculated based on local alignments, coming up with the substitution probabilities in the overall frequency of amino acids. Because of how they’ve been constructed, they’re really sensitive to detecting even subtle changes in structural -- structure or function, so because of that they are the
19:57
best performing of the series of matrices available to you. Historically, in the very early days of bioinformatics there was another series called the PAM matrices, and these were derived at a time when very little sequence information was available. They were mostly based on globular proteins with a lot of extrapolations going on. The only reason I point that out to you is because, especially if you’re using phylogenetic
20:20
programs, you’ll often see them as an option, but I would instead recommend that you pick one of the BLOSUM matrices for your studies. So, bottom line, you can use these matrices to find both closely related and distantly related sequences, so here’s where evolution now comes into play. You’ll notice in that matrix you looked at with all of the numbers on it -- it was called BLOSUM62, so what does that number stand for? It always is BLOSUM and a number.
20:50
What the number represents to you is in this first bullet point. It’s that these -- the -- that particular matrix was built using sequences that shared no more than n percent identity. In the BLOSUM62 case, 62 percent identity, and so in order to remove bias, what ends up happening is we go through a process that looks like this. So, let’s say we have a multiple sequence alignment here. There are a couple of conserved positions, which are noted by the asterisks, and let’s
21:24
say we wanted to construct a BLOSUM80 matrix, so we’re going to look for all of the pair Y’s combinations of these things and cluster together any of the sequences that are 80 percent or more similar to one another, okay? In the next step, we actually cluster those together, so this cluster is represented by this sequence, these three by a single sequence, and these two by a single sequence, and at that point we can calculate the matrix. Why do we do this?
21:53
The reason is really more of the history of sequencing than anything else. There are more and more whole genome sequences available to us, but they’re not very perfectly spread across all of the taxa, as you undoubtedly know in doing your own searches at -- over at NCBI, that you will see the same gene has been sequenced many, many times, so when you look at them, which one is the right one? By -- if you took all of that information together, you would end up biasing the analysis,
22:25
so it’s really just a matter of trying to take this vast amount of sequenced data that we have available to us but narrow it down to the -- in a format that will end up being most biologically informative in order to do these comparisons. Okay so implications. Again, when we do this clustering, we reduce the contribution of closely related sequences, less bias. Bottom line: when we reduce n, you get more distantly related sequences.
22:54
When we increase n, you get more closely related sequences, so this is your first cheat sheet of the morning. Put a star on this slide, please, and basically gives you some recommendations of where you should be looking for your -- for selections. If you’re looking for shorter alignments that are highly similar, you’re going to pick a BLOSUM matrix with a higher number, and you can see here the range of similarity that it’s best at detecting.
23:21
At the bottom, if you used a BLOSUM30 matrix, you would get longer, weaker local alignments, but you would use these to detect more distantly related proteins, things that have much lower sequence similarity. The one in yellow, the BLOSUM62, is the default, and just, again, from testing, this is the most effective in finding all potential similarities. This is the default value when you go to NCBI’s BLAST website.
23:49
Okay, so, you can pick one of these, but I would actually defer to the guru in this field, Steve Altschul, who is the driving force behind the development of BLAST many, many years ago, and what he normally recommends when you don’t really know what you’re looking for is to use a triple BLAST strategy. Pick the default, pick a higher number, pick a lower number, and run your query sequence three times against those three -- using those three different matrices, and that’s what
24:22
I would recommend you do as you get more and more familiar with these methodologies because you’ll start to find that by tweaking these values, you’re going to start to deduce relationships that you didn’t know that existed or didn’t even have any inkling that they existed when you started, so start off with doing these triple strategies until you get more familiar again or unless you really, really know what you’re -- what you’re looking for.
24:45
The bottom line is probably as apparent at this point is that no single matrix is the complete answer for all of these sequence comparisons, and -- so, again, I’ve given you a strategy. There is more information in this current protocols and bioinformatics piece by David Wheeler [spelled phonetically] that talks at much more length about the selection of the matrices, but you have enough to get yourselves started, now. Okay, so that’s the matrices.
25:09
Now, many of you have already done pair-wide sequence alignments, multiple sequence alignments, and a lot of times you see gaps that are inserted, and the gaps are put there to improve the alignment, but the thing you have to keep in mind is that in putting those gaps in to improve the alignment, that you have basically implicitly inserted a biological event, okay? It could be a deletion. It could be an insertion, okay?
25:39
So, it’s not just, “Oh, I’m going to make it pretty by putting some gaps in.” There’s actually something that has to happen in the background to make that actually true. So, because of that, you have to keep the gaps to a reasonable number to not reflect a biologically implausible situation, and the rule of thumb that most people use is about one per twenty. Of course, it depends on what case you’re looking at, but going back the scoring matrices,
26:04
you can now ask yourself, “Well, how do you score these things?” because we’ve been talking about aligning residues, coming up with a value for each particular position. How do you score a nucleotide or an amino acid to a gap, which is essentially nothing? Okay, so, it’s not a match. It’s not a mismatch. What do we do with the numbers? So, there is something called the affine gap penalty, and this is because we are again this is an implicit biological event.
26:31
We have to actually deduct points anytime we put a gap into an alignment. The formula is here, so G plus LN. G is a gap-opening penalty, so just when you put that first gapped position in, it takes a big penalty usually, so on the nucleotide side, you’re deducting five points. On the protein side, you’re deducting 11 points right off the bat. Every time you extend that gap by one more position, you assess this extension penalty—two in the case of nucleotides; one in the case of proteins.
27:04
Okay, so, again, keep in mind, it’s not just about making your alignment pretty. It’s about really making it biologically plausible, okay, and you can change -- these values can all be changed in the BLAST algorithm to make the gaps more or less permissive. Okay, so, with all of that by way of background, we now can finally talk about BLAST, the basic local alignment search tool. This is really the cornerstone, the most important technique we use in bioinformatics because,
27:33
again, it allows us to make these deductions about orthology, paralogy, what have you, really trying to figure out how two proteins are available to other -- and what we can learn by aligning them to one another. So, we’re going to hit this pretty hard today, and it -- because it’s important that you know how to use it properly. Unfortunately, there are many cases in the literature because of the improper use of
27:58
BLAST that have led to really wrong conclusions, and at the end of the day it’s because people either didn’t know how to use it or they didn’t interpret the results properly. The most important thing to take away from this lecture or any of the bioinformatics lectures along the way is that just because you see something in a list of results does not automatically imply biological significance. We’re going to see probabilities. We’re going to see scores.
28:26
We’re going to see all of these things, but a lot of those things that appear in the output list actually have biological relevance, and this is where we’re going to help you figure out which ones do and which ones don’t. So, it’s an election year, so let’s go ahead and take a poll. How many of you have used BLAST before? Okay. How many of you know how it actually works? Okay, a little fewer. Good, good for self-awareness, [laughs] and for those few of you who’ve done it before
28:54
-- I know you, John. How many of you have ever changed the default values? Oh, good, a few of you. Thank you. Okay, now, for the rest of you, nothing to be ashamed of. All will be revealed today, and you’ll be able to consider yourselves expert by the -- experts by the end of today’s lecture. So, okay, back to BLAST, the basic local alignment search tool. What BLAST seeks to do is find things called high-scoring segment pairs, HSPs, and these
29:21
quite simply are pairs of sequences that can be aligned with one another. When they are aligned, based on our scoring matrices, they have the largest maximum score possible that you can’t make that score better by making the alignment longer or shorter. They have to break something called a score threshold, and they can be either gapped or un-gapped. By the -- so since this method is looking for these little pair subsequences, again,
29:49
a local alignment method, the results you get back are not limited to just the best alignment between two sequences. They’ll give you all of the plausible alignments between two sequences. Okay. There are five different programs in the BLAST family, and you can just differentiate them by the last letter, so BLASTN starts with a nucleotide sequence as its query, and it does its searches against a nucleotide database.
30:17
BLASTP, the P stands for proteins, so we’re taking a protein sequence and searching it against a database of protein sequences. BLASTX does something a little bit funky. You start off with a nucleotide sequence, a six-train -- a six-frame translation is performed, and based on those translations the search is done against a protein database, so the comparison is actually proteins versus proteins, even though you started with a nucleotide sequence.
30:45
The others are on the slide for sake of completion, but in your own trials, these are the three that you will probably use more often than not. Okay, so how does this actually work? So, let’s say we have this query sequence of a particular length, and what BLAST does is it doesn’t try to take this whole thing at once and do an alignment. Instead, it starts with just a few -- in this case, amino acids, the default being three,
31:17
so the length of what we call the query word is three, so what happens here is in essence, we’re going to work across different windows in this sequence, so the search is going to start with the GSQ, and then the SQS and the QSL and move across. Just for sake of symmetry on this slide, we’re going to start in the middle here with the P, the Q, and the G. Okay, so we want to now look for the occurrence of that PQG sequence
31:41
in all of the other proteins in the database that we’re searching, but think back to the matrices, right? We want to look -- find the exact matches, but we also want to find places where there are conservative substitutions because those are also going to be biologically informative, so we’re not just looking for the PQG but other things that are close to the PQG, again, using our scores as the metric. So, we have here the neighborhood, the exact match, and a bunch of conservative substitutions,
32:14
in this case just focusing on the middle letter. There’s a score next to each one of these, so in this case the PQG, if it matches exactly, earns 18 points, and so if you go back to the -- to the slide with the matrix and just do a P for P, Q for Q, G for G, that you’d end up with seven plus five plus six, and that’s 18 points. Again, every time we have not an exact match but a conservative substitution, we still
32:40
assess a positive value for that position, but just not as many points, so you’ll see that the scores start to go down. At some point, BLAST will automatically draw a line defining what this neighborhood score threshold is, so where is the fence that we’re putting around the neighborhood. In this case, the value is 13, so anything from 13 and up will move to the next step. Everything from 12 down will be rejected
33:04
Okay, so, again, exact matches and conservative substitutions, and that t value changes from search to search automatically. You can set it, though. All right, so now we’ve got the alignment at the bottom, here, the -- here’s the sequence that we started with. We have it now aligned with a different sequence here that it -- starting again with this PQG, it found here a PMG. The PMG is in our neighborhood. You’ll notice that the middle position actually is not assessed as either an exact match or
33:35
a conservative substitution, that’s what the plus signs stand for, but again, it broke that neighborhood score threshold so it carries through. So just from this you can do a qualitative assessment, looking at this middle line to see how these -- how these look to you, how well they align with one another, but now here’s the time to bring the numbers back into play. All right, so here’s the alignment again. Here we have a graph that’s showing the cumulative score on the y-axis and the length
34:03
of the alignment on the x-axis. So here we have our alignment of three letters, that’s this very first point on the graph here. It’s above our neighborhood score threshold and now we’re just going to keep going out in both directions as far as we can go. So in this particular case, we keep adding -- going in both directions and let’s just assume that we have more matches than we do mismatches, so -- or more conservative substitutions than we do non-conservative substitutions.
34:32
So we keep going up as we extend and here’s this neighborhood score threshold, “S,” okay? Now herein lies pitfall number one, okay? Everything above this score value will be reported to you in the BLAST results. Again, just because it breaks the score threshold it does not imply biological significance, okay? I’ll show you what does in a minute. All right, we’re going to keep going up and up and up, going as far as we can in both
35:02
directions, but now at some point I want you to imagine that we’re going to start mismatching more than we’re matching. So the curve is going to turn downwards because we keep subtracting points all over those negative values in the BLOSUM matrices until we have gone a little bit too far. So BLAST will tell us, ‘All right, you’ve now gone a little bit too far, too many mismatches, too many gap penalties. It’s time to go back to this peak,’ and that peak is going to define the maximum length
35:32
of the alignment between these two sequences, okay? So again, we’re not going to pay attention to this. We’re going to have to come up with a different value and let’s explore a little bit more why that is. Okay. Scores and alignments don’t tell the whole story. All right, so here’s an alignment; this is now a contrived alignment -- we’re going to take a little side road here -- here’s a contrived alignment between two sequences,
35:58
20 amino acids long, and what I did here is this is just a bunch of amino acids that more frequently occur amongst the entire constellation of proteins. In this alignment, also just contrived, same length but using the most rare of the residues, okay? So they both are of length 20, so you might say, ‘Well, they’re both of the same length and so they’re probably both, you know -- you know, just as good as the other one,’ but
36:27
then if you look at the scores -- okay, the top one, because it’s made up of common residues, scores 91. The bottom one, because it’s made up of rare residues, scores 138. So you now have a situation where you can’t just look at the length of the alignment and you can’t just look at the score. You need some metric that will normalize all of this, that you can look at with confidence regardless of what search you’re doing, every single time, to let you know which ones
36:55
of these are true positives and which ones of these are false positives. Okay, so how do we do that? There is a formula. I’m sorry to show you a formula, but this [laughs] formula is called the Karlin-Altschul equation and it uses the score, a normalized version of the score, to calculate a probability value called “E.” Okay, and all you need to know here is that E represents the number of these high-scoring segment pairs, these local alignments that were found purely by chance.
37:27
So this is a measure of the number of false positives and because of that, you want this number to be as small as possible. So lower values of E signify higher similarity, okay? Rule of thumb that I’m going to give you, one of the first rules of thumb as far as BLAST goes, a good starting point for determining what is statistically significant: if you’re doing a comparison between two nucleotide sequences, you want to look for probability values 10 to the minus 6th or smaller.
37:59
If you’re doing a protein-based search -- so this is blastp and blastx; because remember blastx, you’re doing that translation before the comparison is done -- 10 to the minus 3rd or better. Okay, so by using these values you don’t have to worry about what’s being printed out on my screen; you can actually start to deduce which ones are actually true. Okay, so with all of that as far as background goes, we’re going to pretend we’re at the computer.
38:28
The screen captures are all in your handout or on your tablets or on your laptops; I see a nice glow coming from the audience from down here and what I’d like you to do -- we’re going to go through these together but I’d like you, when you get back to your labs, to actually start to practice using these examples as well. Okay, so here we are at the NCBI website; lots of things going on on this page, there’s a lot of good stuff here that takes you to PubMed, a bunch of the databases, but we’re
38:55
going to focus today on BLAST. So the link is right there and if we click on that, that takes us to the BLAST homepage. Let me orient you very quickly here: at the top of the page you can go directly to a set of assembled genomes, at the bottom there are some specialized BLAST methodologies available to you. We’re going to talk about a few of those a little bit later and for right now we’re going to focus on the stuff in the middle here, just basic BLAST, and the second link
39:22
is the protein BLAST. So we’re going to pretend that we’re going to do a query. Again, for your practice, on this webpage -- here is the URL -- are all of the sequences we’re going to go through, both this week and when I see you again in three weeks’ time, and again I would encourage you to just put fingers to keyboards and try all of this out. Okay, so this what the query page looks like. So at the very top, of course you have to have a query sequence that you want to start with.
39:51
This is the sequence from the page I just showed you, so you just paste it into the box. If you don’t want to search the entire sequence you have the ability to just say, ‘All right, I want to start this comparison just using residues 50 through 150,’ if you have a particular region you know that you’re interested in or you could just trim the sequence back. There are a whole bunch of databases that are available to you, so N.R. is NCBI’s non-redundant database.
40:21
The one I want to focus on just for discussion this morning is something called RefSeq and so again let’s take a little bit of an aside. What RefSeq attempts to do is to provide a single reference sequence for each molecule in the central dogma. So a single DNA sequence, a single mRNA sequence, and a single protein sequence and again this goes back to something I already mentioned earlier in the construction of a BLOSUM matrices.
40:46
Let’s say you’re interested in the deleted for colorectal cancer gene in human. If you do that search, you’re going to get back somewhere between 20 and 50 different sequences and you kind of look at it all, you know, sort of befuddled and say, ‘Well, which one is the right one, or at least the canonical one?’ and this is what this is intended to help you to do, is select what is considered to be the canonical sequence.
41:08
Obviously, because each one of these molecules is represented only once, it is non-redundant. There are folks at NCBI who update these records to reflect the current knowledge of sequence data and biology, so that’s very good when you’re looking at the actual GenBank entries. It gives you some information about the biological attributes of the gene, the gene transcript or the protein, wide taxonomic range, and again the ongoing curision [spelled phonetically],
41:34
which I think is the most valuable part of this database. Each entry in GenBank as you undoubtedly know carries something called an accession number, so that in essence is that sequence’s social security number; you can just put that in and you know you’re going to get the same sequence back every time. They carry codes and the way to tell the RefSeq entries from the rest of the entries is by looking for these prefixes.
42:04
So if you see the N series here, “NT” means it’s a genomic contig [spelled phonetically]. If you see the “NM,” it’s an mRNA sequence. The “P” again stands for protein and an “NR” means that it’s a non-coding transcript. What you’ll also see are some entries that arise from genome annotation. So taking these genomic contigs, the NT entries, and using computational methods, you can predict the model mRNAs that come from this sequence and in turn, the model proteins that come
42:37
from these model mRNAs. So this is probably the most important thing to take away from this slide is when you’re looking at these that start with an N, these have been directly sequenced or directly determined, experimentally determined. The ones with the Xs arise from predictions, okay? So if you get back a set of results, if the only thing available to you are the Xs, that’s fine; keep it in the back of your head though that they haven’t necessarily been verified.
43:05
Always go for the Ns when you can. The other nice thing about using RefSeq is, because of the non-redundant nature of these databases, you get just a tighter, neater set of results back. Okay, but for purposes of the example I’m just going to leave this at NR. You can limit your searches by organism or by taxonomic group, so you can say “human,” “bilateria,” whatever you want to put in there or you can just leave it blank and it will search everything.
43:32
The program selection is blastp, right down here, and this is where people usually stop, they’ll just click the button and go, but down here tucked away are these -- is a little popup that says “algorithm parameters,” so it will just expand the page now and here’s where a lot of the numbers we just talked about come into play. The first one just gives you the maximum number of target sequences; what this means is how
43:58
many sequences will BLAST return to you based on the search that you’re issuing, based on your search query. The default is 100; all right, in this case I’m changing it to 250. I like to set it on the high end because I don’t want to miss anything, right? [laughs] So I would rather put it -- put it on the higher end and then determine which ones I want to look at. There’s an expect threshold; this is the E-value threshold, so again, those probability values that we talked about today.
44:27
BLAST sets this at 10. Again, I’m going to leave it at 10 for the example, just for illustrative purposes; you could change it to the values that I gave you, but I’ll caution you as to why you may not want to do that later. Again, going back to manual inspection. The word size is three; you can change the word size if you want. Selection of matrices, so you can see everything that’s available to you here. Underneath you can change the gap costs as well, so if you want to make it less gap-prone,
45:02
you can increase the penalties; if you want to be more permissive, you can decrease the penalties. Filters so we can filter out things called low-complexity regions. So what are those things? These are regions of biased composition. BLAST depends on the fact that the sequences that it is analyzing is of relatively uniform composition, but occasionally you’ll see proteins where -- or even nucleotide sequences where there are runs of a single amino acid or runs of a particular nucleotide and that
45:35
ends up messing BLAST up; it often leads to false positives. So things like homopolymeric runs, short period repeats, and just the subtle overrepresentation of several residues. So I would advise that you filter them out, but it’s not enabled by default, so we’re going to just put a checkbox over there. Finally at the bottom here there’s a little checkbox that says “show results in a new window” and I like to do that so I don’t have to keep moving back and forth between
46:02
the pages; it just makes navigation a little bit easier. So with all of those things, we’re going to finally now go ahead and click the BLAST button and see what we get back. It’s not going to be quite that fast, you might have to get coffee certain mornings, but this morning we’re good. Okay, so let me orient you to this screen as well. The first thing I want to point out is right here above this very multi-colored table is says “distribution of 200 BLAST hits on the query sequence.”
46:32
You’ll recall that the default to return is 100 sequences, so if we had not changed our value to 250 we would have missed 100 of the results that came out of this analysis, so again, it’s good that we changed to a higher number. There is a section up here where your sequence is represented by this black bar and what’s right below it: it has now identified a protein superfamily that this belongs to and we’re
46:59
going to talk more about that domain analysis when I see you next and we talk about domain architectures. Let’s pretend to scroll down a little bit and let me decipher what this graph actually means. There’s a color key at the top showing you the alignment scores, not the probabilities, the scores. So each one of these lines, if -- let’s say all of the red ones, it means that they scored more than 200. Again, not really relevant to our analysis and being able to make conclusions.
47:29
They’re all in descending score order. If you mouse over any of them, this box up here will show you what the found sequence actually is, the identity of that sequence. Anytime you see the little gaps here, it -- so here we have an insertion, here we have a deletion, just that little hash mark that goes in -- it’s just a gap in the alignment with what you started with, with the subject -- I’m sorry, with the sequence it found, the subject.
47:59
So it mightn’t be a case where you’ve got more than one high-scoring segment pair. Here we’ve got a sequence where there’s one, two, three, and four different high-scoring segment pairs or it could be a massed region. Again, looking for those regions of biased composition and we’ll see that in a little bit more detail in a couple of slides. All right, again imagining we’re scrolling down and this is our hit list, this is what actually came back.
48:25
You have the identities, the definition line for each one of the sequences, the score values are here, but again we’re going to ignore those and instead we’re going to look at the column that says “E-value,” okay? Now you’ll see that there’s a bunch of zeroes here, okay? Now strictly speaking, in computer speak that means that this is 10 to the minus 1,000 or less. Basically it’s just vanishingly small, it’s just too -- the values get so small at some
48:56
point that they’re not even worth showing on here, plus there’s only so many -- so much room in the field. Let’s take a look further down, there’s one down there that’s [unintelligible] minus 179 and that just means 8 times 10 to the minus 179, okay? So really, really small numbers at the top of our list and what we can -- and there’s percent identity values here as well. As you can see, especially with the ones at the top, they’re all very close to each other, 100 percent, 98 percent and so on.
49:25
Let’s scroll to the bottom of the list, okay, and here’s where we’re going to start applying your -- some rules to help you deduce which ones of these are actually the ones that you want to consider, that might actually have some biological relevance. So we’re going to draw a line down here and we’re going to reject everything that is below the desired threshold, again, from the rules I gave you before, 10 to the minus 3rd or higher in this case.
49:53
So now in this case when you look at these they all meet the criteria, so we’re going to go and accept them but you’ll see here I say “for now,” because there’s another rule coming that we’re going to use to let biology be your guide. We’re going to take a look at one of these alignments, these high-scoring segment pairs in a little bit more detail here. At the very top it gives you the identity of what this found protein was and links back to its corresponding GenBank entry.
50:24
In this case the sequence we started with, the query, it started -- this alignment starts at position 17. The subject, the one that it found, also starts at position 17 and each line is -- shows you what the position numbers are, so this alignment goes from 17 to 704. Okay. At the top you have some information on identities and positives, okay? The identities are how many exact matches are there between these two sequences and
50:58
so 688 of the 688 positions lined up with each other, so this is 100 percent. The positives are the same as well. The difference between these two: “identities” mean exactly that, the exact matches; the “positives” include the conservative substitutions, okay? So those numbers will not always match and in this case we didn’t introduce any gaps. So here’s where the next rule comes into play. So I’ve given you a probability rule, now I’m going to give you an identity rule.
51:25
You’re going to be looking at this value -- if you’re doing a protein-protein comparison, again blastp or blastx, you want that value to be at least 25 percent. If you’re doing a nucleotide search, you want that value to be at least 70 percent for starters, okay? All right, you see a bunch of these that are in lowercase here and these are representative of the low-complexity regions that I talked about before, these regions of biased composition,
51:52
so it just points out to you where they actually are. Okay. In the second example -- this is still the same entry, this is still the same hit, here’s the 704 that you saw in the previous slide -- here, above this little statistics row -- in the last one, we’re going to just go back for a minute -- you notice that there was -- the identity of the protein was given to you. Here you don’t see that, so it’s another high-scoring segment pair for the same two proteins.
52:19
The one you started with aligned to the one that it found here, so this is just another part of the molecule that aligned; in this case the alignment was from 777 to 1,374. So here it says “range two,” that just tells you very quickly that the second high-scoring segment pair was identified. Okay, in this particular case we have a gap that’s introduced here as well, so you can see where that gap is and you can see also that that’s reflected in the numbers up
52:48
here, that it’s four percent of the total alignment. Okay, so let’s bring the example all the way around; so here’s our overview of our results again. The two lines that we saw, that gave us our statistics -- so for the first high-scoring segment pair and the second high-scoring segment pair, it beat our E-value that we were looking at okay? Again, percent identity: these are proteins, so 25 or better, so we’ve got 195, we’ll
53:18
accept those as well and just -- you can see if you were to do this on the screen and mouse over it, it would show you that this was the particular hit that gave rise to these two high-scoring segment pairs, okay? So it’s pretty much as easy as that once you understand what the rules are, right? So the manipulations of cutting and pasting and clicking buttons, that’s not the hard part; it’s really understanding what you should be doing with respect to all of these default values.
53:45
The reason they’re the default values is in most cases they’ll work, but now hopefully you understand a little bit better what can possibly go wrong if you don’t pick the right values, especially based on the kinds of biological questions you’re trying to answer. So the scoring matrices come into play if you’re looking, again, for things that are more related or less related, so closely or distantly. The distantly related ones are going to give you some great phylogenetic information about
54:08
the origin over different genes and proteins. The closely related ones are going to give you better information about the particular properties of a protein family and what makes -- what is distinctive about that particular protein family and these are the rules that are going to help you find all of that out. So again, second cheat sheet of the day, put a star on this one as well: here are the rules for the E-value and the sequence identity again for both the nucleotides and the proteins.
54:34
Now these are guidelines, okay? These are not holy rules, these are not -- you know, these are really just a starting point, so I don’t want you to use these cutoffs blindly. This is just, again, someplace to start until you get more comfortable with using BLAST. You should also pay attention to alignments on either side of the dividing line. So you’ll remember we drew a red line, right, but take a look -- if there are things below
54:57
that line that might actually make biological sense to you -- and it’ll happen more often than not that it might just be under the cutoff, but you have biological knowledge that tells you that the things under that line are related. Biology always trumps these results, okay, because remember these are predictions, okay? So when you have biological knowledge, you use that biological knowledge; so never ignore the biology that you know, okay? All right, so that’s it for BLAST proper.
55:27
There were some other things on the BLAST page that I want to point out to you. So let’s say instead of just having a sequence of interest and trying to deduce new relationships by comparing it to a database of sequences, instead let’s say you have two sequences, A and B, that you just want to align with each other. You’re doing some study and you just want to see how well they align, so we’re going to use something called “BLAST two sequences” to do that.
55:53
We’re going to come back to our BLAST homepage and this particular program is found in the lower part here: bl2seq: align two or more sequences using BLAST. This page looks very similar to the one that we saw before. The only difference is that when you check the little box that says “align two or more sequences,” this part of the screen will then expand, where you can put the second sequence in. So we’re going to take a sequence here called “SOX10,” another one called “the sex-determining
56:24
region Y.” These are related to one another, we’re just going to see how related they actually are and where. So again, sequence A, sequence B. As before we could define query subranges on these things. We’re going to expand the algorithm parameter box one more time to see what’s down here, so again, maximum target sequence is 100; in this case we don’t care because we know that we’re only going to get one back because we’re only doing a pair-wise alignment.
56:51
Changing the E-value here is also not that important because the sequences were preselected, so you were already starting from some place of biological knowledge, but you’ll still get back a measure of whether or not the alignments are statistically similar. So word sizes, all of those things remain the same. Below, again you can change the matrices, you can filter for the low-complexity regions again, but we’re just going to go ahead, check the “show results in new window”
57:22
box off, and go ahead and do the BLAST search and what you end up getting back looks very similar to what we saw before. Here is the overview. Again, you started -- you’re only comparing against one sequence, so you only see one line here. One is much shorter than the other. If we look at the results table here it gives us an E-value, one times 10 to the minus 26, so that meets our E-value criterion. The percent identity of the match is 46 so that matches our identity criterion as well
57:55
and here are the two alignments. So the first scoring -- high-scoring segment pair, the second high-scoring segment pair and you can just see for yourselves how well these align to one another. This one is pretty short, I’m not sure I would buy that one, but here -- you know, this is the region of interest here that’s biologically important between these two sequences. Okay, so pretty easy way to do it. Now we focused on proteins for most of the morning because I think it’s easier to understand
58:22
how BLAST works because of how the matrices work as well, but of course we need to be able to do nucleotide searches as well. We’re not going to spend much time on this this morning other than for me to orient you to what’s where. So we come back to this page and again, the nucleotide BLASTs are in this second section within the basic BLAST section and when you move down to the program selection you can pick one of three different nucleotide algorithms to do this comparison.
58:56
One is called “MegaBLAST,” one is called “discontiguous MegaBLAST,” and the other one is called “blastn.” So what’s the difference between the three -- and they’re summarized on this slide for you. So MegaBLAST, which is the default, is really optimized for very long or highly similar sequences, so when you want to do a search against complete genomes or against chromosome length contigs. So this really allows you to map transcripts that you might have, mRNA sequences, things
59:27
like that, onto a sequence in seconds. You’ll see the word length here is 28, so you get longer and highly similar alignments and that’s what you want in this particular case. You’re looking for almost exact matches or just exact matches and it runs a lot faster because that word length is longer. The middle block is probably what you’re going to be doing more often when it comes to searches at the nucleotide level.
59:55
The results between these two methods will actually be very similar, but MegaBLAST is going to run faster in relative terms. Pay attention here to the kinds of gap penalties here. Discontiguous MegaBLAST is using the affine gap penalty that we talked about before. In this case MegaBLAST uses a linear BLAST penalty, so there’s no penalty for opening the gap, it’s just a defined number of points for each position in the gap, okay? Finally the old stalwart, blastn.
1:00:25
It’s best for finding short, nearly exact matches and you can see the values that are being applied here as well. Okay. So that’s really all we’re going to talk about on the nucleotide side today, because it works exactly the same as what we talked about on the protein side and again, you have different programs available to you, but it’s going to come down between BLAST and discontiguous MegaBLAST as far as your choices go on this side of the house.
1:00:52
The last thing I want to talk to you about this morning is something called Blat. So we are now going to leave the Lister Hill Center on the other side of campus and move ourselves over to the sunny shores of California and a tool that’s made available through the genome center at the University of California at Santa Cruz called Blat, and what Blat just stands for is the BLAST-like Alignment Tool and the method is somewhat similar to MegaBLAST
1:01:19
in that it really is optimized to align longer nucleotide sequences, things that are more than 40 in length, having 95 percent or more sequence similarity and because of the way these settings are set, you can find exact matches reliably down to a length of 33, so really the method of choice when you’re looking for exact matches in nucleotide databases. It’s also much faster, 500 times faster, than BLAST for doing either mRNA or DNA searches.
1:01:51
It may however miss divergent or shorter sequence alignments, so again here’s a case where you want to use BLAST and you want to use Blat depending what you’re looking for. You can use it on protein sequences; I never have, I feel that blastp is more efficient for doing those kinds of searches. Okay, so when do you use this? You’re looking to characterize some unknown gene or sequence fragment, so where does it
1:02:17
live on a particular chromosome, finding its genomic coordinates, determining gene structure, so the presence or absence of exons, what other information might be available in intergenic regions, and looking for markers of interest in the vicinity of a sequence. Obviously it helps you find highly similar or identical sequences based on what I’ve already told you, so you can align mRNA sequences onto a genome assembly, find other gene family
1:02:44
members, or do some cross species alignment to identify putative homologs, going back to the evolutionary question. It also allows you to display a specific sequence as a separate track within the browser, so we now start to introduce the concept of a track; we’re going to see an example of that shortly. Dr. Wolfsberg [spelled phonetically] is going to be spending quite a bit of time discussing the use of this browser when she presents her lecture in two weeks’ time, but let
1:03:11
me give you a first gentle introduction to it. So here we are at the UCSC genome bioinformatics website. You can get to the Blat program by just clicking on the Blat link which is in the left sidebar; takes you to a very simple search page, much simpler than what we saw with BLAST. So I’ve gone ahead here and just pasted in a sequence of a cDNA clone that comes from the cancer genome anatomy project. This is from Macaca mulatta, so we’re going to change here the genome to rhesus.
1:03:43
The assembly automatically will change for you, it’s as a DNA sequence and you have two choices here, “submit” or “I’m feeling lucky.” So if you’re feeling lucky, what this is going to do is just return to you the highest-scoring alignment between this sequence and the ones from -- the ones against the rhesus genome. Instead what we’re going to do is click on the submit button. Okay. All right, so here is what you get back.
1:04:10
You’ll notice that for each one of these -- here is what we started with, you have a score for each one of them, it gives you a start and end position, okay, then puts it -- tells you which chromosome in the rhesus genome it actually -- the match is actually on and how long the alignment is. So since this is the longest alignment, that’s the one we’re going to consider for this example. In order to visualize this within the UCSC genome browser we’re just going to click
1:04:36
on the browser link all the way over there on the left and this is what you get. If we had clicked the I’m feeling lucky button you would have come directly to this screen. So right away we can see the chromosomal position of our particular sequence; so this is on chromosome six, this is the sequence we -- where is -- this is the sequence we started with, your sequence from Blat search, and it in turn now has found a whole bunch of non-rhesus
1:05:03
RefSeq genes -- so again, that non-redundant database -- from a wide variety of organisms. So very quickly you can now start to identify orthologous genes to the one that we started with and actually the identity of that -- those orthologous genes. So if we started with an unknown, we now can actually assign it to a gene family. We’re not going to talk about this much more than that; let me just show you one other view of this.
1:05:29
If we click on the “details” button we can look at a more traditional view of the alignment. So here is what we started with. You’ll see a whole bunch of colors going on and the key is deciphered for you at the top; so matching bases in cDNA and genomic sequences are colored blue and capitalized, light blue bases mark the boundaries of gaps in either sequences, often splice sites, and if we just keep scrolling down we now just see what the pair-wise alignment looks more
1:05:55
traditionally -- more reminiscent of the BLAST view, okay? So that is it for today as far as what we’re going to cover. We’re going to build upon these concepts when I’m back with you in week four, in discussing how we can go beyond just looking at these letter-by-letter alignments to find the relationships between gene and protein sequences. Next week Eric Green will be behind the podium; he is going to give his normal tour de force
1:06:22
and give you a really good, comprehensive view of the genomic landscape as of this year. It is -- it has always one of the most wonderful lectures in this series because it really gives you a great overview of everything going on in genomics currently, but also giving you a lot of the history of how we’ve gotten to the position where we can do a lot of the great work that we currently do. So with that, I’m happy to take questions at the podium.
1:06:48
Thanks for joining us today and we’ll see you again next week. [applause] [end of transcript]
Here to help you learn
Sign in to keep a private conversation with your course tutor.
Ask me anything about the material…
Sign in to chat