The following content is provided
under a Creative Commons license.
Your support will help MIT OpenCourseWare continue
to offer high quality educational resources for free.
To make a donation or view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at ocw.mit.edu.
PROFESSOR: All right.
We should probably get started.
So RNA plays important regulatory and catalytic roles
in biology, and so it's important to understand
its function.
And so that's going to be the main theme of today's lecture.
But before we get to that, I wanted
to briefly review what we went over last time.
So we talked about hidden Markov models, some
of the terminology, thinking of them as generative models,
terminology of the different types of parameters,
the initiation probabilities and transition probabilities
and so forth.
And Viterbi algorithm, just sort of the core algorithm
used whenever you apply HMMs.
Essentially, you always use the Viterbi algorithm.
And then we gave as an example the CpG Island HMM,
which is admittedly a bit of a toy example.
It's not really used in practice,
that illustrates the principles.
And then today we're going to talk about a couple
of real world HMMs.
But before we get to that, I just
wanted to-- sort of toward the end,
we talked about the computational complexity
of the algorithm, and concluded that if you have a case state
HMM run on a sequence of length L, it's order k squared L.
And this diagram is helpful to many people in sort
of thinking about that.
So you can have transitions from any state--
for example, from this state-- to any of the other five
states, and there's five-state HMM.
And when you're doing the Viterbi,
you have to maximize over the five possible input transitions
into each state.
And so the full set of computations
that you have to do from going from position i to i plus 1
is k squared.
Does that make sense?
And then there's L different transitions you have to do,
so it's k squared L.
Any questions about that?
OK.
All right and, so the example that we gave is shown here.
And what we did was to take an example sort of where you could
sort of see the answer-- not immediately see it,
but if we're thinking about it a little, figure out the answer.
And then we talked about how the Viterbi algorithm actually
works, and why it makes the transitions at the right place.
It seems to intuitively like it would make a transition later,
but actually transitions at the right place.
And one way to think about that is
that these are not hard and fast decisions because you're
optimizing two different paths.
At every state, you're considering two possibilities.
And so you explore the possibility of-- the first time
you hit a c, you explore the possibility of transitioning
from genome to island, but you're not
confirming whether you're going to do that yet until you get
to the end and see whether that path ends up having a higher
probability at the end of the sequence than the alternative.
So that's sort of one way of thinking about that.
Any questions about this sort of thing,
how to understand when a transition will be made?
And I want to emphasize, for this simple HMM,
we talked about you can kind of see
what the answer's going to be.
But if you have any HMM, any sort of interesting real world
HMM with multiple states, there's
no way you're going to be able to see it.
Maybe you could guess what the answer might be,
but you're not going to be able to be confident of what that
is, which is why you have to actually implement it.
All right, good.
Let's talk about a couple of real world HMMs.
So I mentioned gene finding.
That's been a popular application of HMMs,
both in prokaryotes and eukaryotes.
There's some examples discussed in the text.
Another very popular application are so-called profile HMMs.
And so this is a hidden Markov model
that's made based on a multiple alignment of proteins which
have a related function or share a common domain.
For example, there's a database called
Pfam, which includes profile HMMs for hundreds
of different types of protein domains.
And so once you have many dozens or hundreds or thousands
of examples of a protein domain, you
can learn lots of things about it--
not just what the frequencies of each residue
are in each position, but how likely
you are to have an insertion at each position.
And if you do have an insertion, what
types of amino acid residues are likely to be inserted
in that position, and how often you
are likely to have a deletion at each position
in the multiple alignment.
And so the challenge then is to take a query protein
and to thread it through all of these profile HMMs and ask,
does it have a significant match to any of them?
And so that's basically how Pfam works.
And the nice thing about HMMs is that they allow you
to-- if you want to have the same probability
of an insertion at each position in your multiple alignment,
you can do that.
But if you have enough data to observe that there's
a five-fold higher likelihood of having an insertion at position
three in a multiple alignment than there is at position two,
you can put that in.
You just change those probabilities.
So in this HMM, each of the hidden states
is either an M state, which is a match state, or an I state,
or an insert state.
And so those will emit actual amino acid residues.
Or it could be a delete state, which
is thought of as emitting a dash, a placeholder
in the multiple alignment.
So these are also widely used.
And then one of my favorite examples-- it's fairly simple,
but it turns out to be quite useful--
is the so-called TMHMM for prediction
of transmembrane helices in protein.
So we know that many, especially eukaryotic proteins,
are embedded in membranes.
And there's one famous family of seven transmembrane helix
proteins, and there are others that
have one or a few transmembrane helices.
And knowing that a protein has at least one transmembrane
helix is very useful in terms of predicting its function.
You predict it's localization.
And knowing that it's a seven transmembrane helix protein
is also useful.
And so you want to predict whether the protein has
transmembrane helices and what their orientation is.
That is, proteins can have their end terminus
either inside the cell or outside the cell.
And then of course, where exactly those helices are.
And this program has about a 97% accuracy,
according to [? the author. ?] So it works very well.
So what properties do you think--
we said before that you have to have
strongly different emission probabilities
in the different hidden states to have a chance of being
able to predict things accurately.
So what properties do you think are captured
in a model of transmembrane helices?
What types of emission probabilities
would you when you have for the different states in this model?
Anyone?
So for this protein, what kind of residues
would you have in here?
Oops, sorry.
I'm having trouble with this thing.
All right, here in the middle of the membrane,
what kind of residues are you going to see there?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Those are going to be hydrophobic.
Exactly.
And what about right where the helix emerges
from the membrane?
[INAUDIBLE] charge residue's there to kind of anchor
it and prevent it from sliding back into membrane.
And then in general, both on the exterior and interior,
you'll tend to have more hydrophilic residues.
So that's sort of the basis of TMHMM.
So this is the structure.
And you'll notice that these are not exactly the hidden states
that correspond to individual amino acid residues.
These are like meta states, just to illustrate
the overall structure.
I'll show you the actual states on the next slide.
But these were the types of states
that the author, Anders [? Crow ?], decided to model.
So he has sort of a-- focuses here on the helix core.
There's also a cytoplasmic cap and a non-cytoplasmic cap.
Oops, didn't mean that.
And then there's sort of a globular domain on each side--
both on the cytoplasmic side, or you
could have one on the non-cytoplasmic side.
OK, so there's going to be different compositions in each
of these regions.
Now one of the things we talked about with HMMs
is that if you were-- now let's think about the helix core.
The simplest model you might think of
would be to have sort of a helix state,
and then to allow that state to recur to itself.
OK, so this type of thing where you then have some transition
to some sort of cap state after, this would allow you
to model helices of any length.
But now how long are transmembrane helices?
What does that distribution look like?
Anyone have an idea?
There's a certain physical dimension.
[INAUDIBLE]
It takes a certain number residues to get across here,
and then that number is about 20-ish.
So transmembrane helices tend to be sort of
on the order of 20 plus or minus a few.
And so it's totally unrealistic to have a transmembrane helix
that's, like, five residues long.
So if you run this algorithm in generative mode,
what distribution of helix lengths will you produce?
We're running in generative mode where
we're going to let, remember, to generate
a series of hidden states and then
associated amino acid sequences.
It's coming from some, let's say-- I don't know.
What kind of states are there here? [INAUDIBLE] plasmic.
Let's say goes into helix, hangs out here.
I'm sorry, is there an answer to this question?
Anyone?
I don't know how long-- if I let it run,
it'll generate a random number.
It depends on what this probability is here.
Let's call this probability p, and then this
would be 1 minus p.
OK, so obviously if 1 minus p is bigger,
it'll tend to produce longer helices.
But in general, what is the shape
of the distribution there of consecutive helical states
that this model will generate?
AUDIENCE: Binomial.
PROFESSOR: Binomial.
OK, can you explain why?
AUDIENCE: Because the helix would
have to have probable-- the helix of length n
would occur 1 minus p to the n power.
PROFESSOR: OK, so a helix of length 10 with a probability
of then, say, let's call it L, for the length of the helix,
equals n is 1 minus p to the n, right?
Is that binomial?
Someone else?
AUDIENCE: Yeah.
Is it a negative binomial?
PROFESSOR: Negative binomial.
OK.
AUDIENCE: [INAUDIBLE] states and a helix state before moving out
[INAUDIBLE].
PROFESSOR: Yeah.
So the distribution is going to be like that.
You have to stay in here for n and then leave.
So this is the simplest-- you can
have special cases of binomial and negative binomial.
But in general, this distribution
is called the geometric distribution.
Or a continuous version would be the exponential distribution.
So what is the shape of this distribution?
If I were to plot n down here on this axis, and the probability
that L equals n on this axis, what kind of shape--
could someone draw in the air?
So you had up and then down?
OK, so actually, it's going to be just down.
Like that, right?
Because as n increases, this goes down
because 1 minus p is less than 1.
So it just steadily goes down.
And what is the mean of this distribution?
Anyone remember this?
Yeah, so there's sort of two versions of this
that you'll see.
One of them is the 1 minus p n minus 1 p, and one of them
is this.
And so this is the number of failures before a success,
if you will.
Successes lead to the helix.
And this is the number of trials till the first success.
So one of them has a mean that's 1/p,
and the other has a mean that's 1 minus p over p.
So usually, p is small, and so those are about the same.
So 1/p.
You could think that 1/p is roughly right.
And so if we were to model transmembrane helices,
and if transmembrane heresies are about--
I said about 20 residues long-- you
would set p to what value to get the right mean?
AUDIENCE: 0.05.
PROFESSOR: Yeah.
0.05.
1/20, so that 1 over that will be about 20, right?
And then 1 minus p would, of course, be 0.9.
So if I were to do that, I would get a distribution
that looks about like this with a mean of 20.
But if I were to then look at real transmembrane helices
and look at their distribution, I
would see something totally different.
It would probably look like that.
It would have a mean around 20.
But the probability of anything less than 15 would be 0.
That's too short.
It can't go across the membrane.
And then again, you don't have ones that are 40.
They don't kind of wiggle around in there and then come out.
They tend to just go straight across.
So there's a problem here.
You can see that if you want to make a more accurate model,
you want to not only get the right emission probabilities
with the right probabilities of hydrophobics and hydrophilics
and the different states, but you also
want to get the length right.
And so the trick that-- well, actually, yeah.
Can anyone think of tricks to get the right length
distribution here?
How do we do better than this?
Basically, hidden Markov models where
you have a state that will recur to itself,
it will always be a geometric distribution.
The only choice you have is what is that probability.
And so you can get any mean you want,
but you always get this shape.
So if you want a more general shape,
what are some tricks that you could do?
How could you change the model?
any ideas?
Yeah, go ahead.
AUDIENCE: [INAUDIBLE] have multiple helix states.
PROFESSOR: Multiple helix states.
OK.
How many?
AUDIENCE: Proportional to the length we want, [INAUDIBLE].
PROFESSOR: Like one for each possible length.
AUDIENCE: It'd be less than one length.
PROFESSOR: Or less than one.
OK.
So you could have something like-- I mean,
let's say you have like this.
Helix begin-- or, helix 1, helix 2.
You allow each of these to recur to themselves.
What does that get you?
This actually gets you something a little bit better.
It gives you a little bit about of-- it's more like that.
So that's better.
But if I want to get the exact distribution, then actually
one-- so this is the solution that the authors actually used.
They made essentially 25 different helix states,
and then they allowed various different transitions here.
So it's a larger arbitrary here, but they
have this special state three that can kind of take a jump.
So it can just continue on to four,
and that'll make your maximum length helix core.
Or it can skip one, go to five, and that'll
make a helix core that's one residue shorter than that,
or it can skip two, and so forth.
And you can set any probabilities
you want on these transitions.
As so you can fit basically an arbitrary distribution
within a fixed range of lengths that's
determined by how many states you have.
OK, so they really wanted to get the length distribution right,
and that's what they did.
What's the cost of this?
What's the downside?
Simona?
AUDIENCE: I was just going to ask,
it looks like from this your minimum helix
length could be four.
PROFESSOR: Yeah.
That's a good question.
Well, we don't know what the probabilities-- they say
said on that.
Well, did they really mean that?
And also, that's only the core, and maybe these cap things
can be-- yeah, that seems a little short to me.
So yeah, I agree.
I'm not sure.
It could just be for the sake of illustration,
but they don't actually use those.
But anyway, I'll probably have to read the paper.
I haven't read this paper for many years
so I don't remember exactly the answer to that.
But I have a citation.
You can look it up if you're curious.
But the main point I wanted to make with this
is just that by setting an arbitrary number of states
and putting in possible transitions between them,
you can actually construct any length of distribution
you want.
But there is a downside, and what is that downside?
AUDIENCE: Computational cost.
PROFESSOR: Yeah, the computational cost.
Instead of having one helix state,
now we've got 25 or something.
So and the time goes up by the square of the number of states,
so it's going to run slower.
And you also have to estimate all these parameters.
OK, so here's an example of the output of the TMHMM
program for a mouse chloride channel gene, CLC6.
So the program predicts that there
are seven transmembrane helices, as shown
by these little red blocks here.
You can see they're all about the same-- about 20 or so--
and that the program starts outside and ends inside.
So let's say you were going to do some experiments
on this protein to test this prediction.
So one of the types of experiments people do
is they put some sort of modifiable
or modified residue into one of the spaces
between the transmembrane helices.
And then you can test, by modifying this cell
with something that's a non-permeable chemical,
can you modify that protein?
So only if that stretches on the outside of the cell
will you be able to predict it.
So that's a way of testing the topology.
So if you were doing those types of experiments,
you might actually-- like maybe you're not sure
if every transmembrane helix is correct.
There could be some where the boundaries were
a little off, or even a wrong helix.
And so one of the things that you often
want with a prediction is not only
to know what is the optimal or most likely prediction,
but also how confident is the algorithm in each
of the parts of its prediction.
How confident is it in the location of transmembrane helix
three or the probability that actually there
is a transmembrane helix three.
And so the way that this program does that is using something
called the forward-backward algorithm.
So those of you who read the Rabener tutorial,
it's described pretty well there.
The basic idea is that I mentioned
that this Po-- the probability of the observable sequence
summing over all possible HMM structures
or all possible sequences of hidden states--
that is possible to calculate.
And the way that you do it is you
run an algorithm that's similar to the Viterbi,
but instead of taking the maximum entering
each hidden state at intermediate positions,
you sum those inputs.
So you just do the sum at every point.
And it turns out that will calculate the sum of the two
values at the end-- or the k values at the end
will be equal to the sum of the probabilities of generating
the observable sequence over all possible sequences
of hidden states.
OK, so that's useful.
And then you can also run it backwards.
There's no reason it has to be only going in one direction.
And so what you do is you run these sort of summing versions
of the Viterbi in both the forward direction
and also run one in the backward direction.
And then you take a particular position here--
like let's say this is your helix state, for example.
And we're interested in this position
somewhere in the middle of the protein.
Is that a helix or not?
And so basically you take the value
that you get here from the forward
in your forward algorithm and the value
that you get here in the backward algorithm,
and multiply those two together, and divide by this Po.
And that gives you the probability.
So that ends up being a way of calculating
the sum of all the parses that go
through this particular position i in the sequence
in that particular state.
I mean, I realize that may not have been totally clear,
and I don't want to take more time to totally go into it,
but it is pretty well described and Rabener.
And I'll just give you an example.
So if you're motivated, please take a look at that.
And if you have further questions,
I'd be happy to discuss during office hours next week.
And this is what it looks like for this particular protein.
So you get something called the posterior probability, which
is the sum of the probabilities of all the parses.
And they've plotted it for the particular state that
is in the Viterbi path, that is in the optimal parse--
so for example, in blue here.
Well, actually, they've done it for all the different states
here.
So blue is the probability that you're outside.
OK, so it's very, very confident that the end terminus
of the protein is outside the cell.
It's very, very confident in the locations
of transmembrane helices one and two.
It actually more often than not thinks
there's actually a third helix right here,
but that didn't make it in the optional parse.
That actually occurs in the majority of parses,
but not in the optimal.
And it's probably because it would then cause other things
to be flipped later on if you had transmembrane helix there.
It's not sure whether there's a helix there or not,
but then it's confident in this one.
OK, so this gives you an idea.
Now if you wanted to do some sort of test of the prediction,
you want to test probably first the higher confidence
predictions, so you might do something right here.
Or if maybe from experience you know
that when it has a probability that's that high,
it's always right, so there's no point testing it.
So you should test one of these kind of less confident regions.
So this actually makes the prediction much more useful
to have some degree of confidence assigned
to each part of the prediction.
So for the remainder of today, I want
to turn to the topic of RNA secondary structure.
So at the beginning, I will sort of
get through some nomenclature.
And then to motivate the topic, give some biological examples
of RNA structure.
Gives me an excuse to show some pretty pictures of structure.
And then we'll talk about two approaches which
are two of the most widely used approaches toward predicting
structure.
So using evolution to predict structure
by method of co-variations, which works well when you
have many homologous sequences.
And then using sort of first principles
thermodynamics to predict secondary structure
by energy minimization where obviously you
don't need to have a homologous sequence present.
And the nature biotechnology primer
on RNA folding that I recommended
is a good intro to the energy minimization approach.
So what is RNA secondary structure?
So you all know that RNAs, like proteins,
have a three-dimensional tertiary fold structure that,
in many cases, determines their function.
But there's also sort of a simpler representation
of this structure where you just describe which pairs of bases
are hydrogen bonded to one other.
OK, and so for RNA-- so it's a famous example
of an RNA structure, this sort of clover leaf structure
that all tRNAs have.
The secondary structure of the tRNA is the set of base pairs.
So it's this base pair here between the first base
and this one toward the end, and then
base right here, and so forth.
And so if you specify all those base pairs,
then you can then draw a picture like this, which gives you
a good idea of what parts of the RNA molecule are accessible.
So for example, it won't tell you
where the anticodon loop is, which
is sort of the business end of the tRNA.
But it narrows it down to three possibilities.
You might consider that, or that, or down here.
It's unlikely to be something in here
because these bases are already paired.
They can't pair to message.
So it gives you sort of a first approximation toward the 3D
structure, and so it's quite useful.
So how do we represent secondary structure?
So there's a few different common representations
that you'll see.
So one is-- and this is sort of a computer-friendly but not
terribly human-friendly representation,
I would say-- is this sort of dot
in parentheses notation here.
So the dot is an unpaired base and the parenthesis
is a paired base.
And how do you know-- chalk is sort of non-uniformly
distributed here-- so if you have a structure like this
and you have these three parentheses, what
are they paired to?
Well, you don't know yet until you get further down.
And then each left parenthesis has
to have a right parenthesis somewhere.
So now if we see this, then we know
that there are two unpaired bases here,
and then there's going to be three
in a row that are paired-- these guys.
We don't know what they're paired to yet.
Then there's going to be a five base pair loop, maybe
a little pentagon type thing.
Two, three, four-- oops-- four, five.
And this one would be the right parentheses
that pair with the left parentheses over here.
I should probably draw this coming out
to make it clearer that it's not paired.
So this notation you can convert to this.
So after a while, it's relatively easy to do this,
except when they're super long.
So that's what the left part of that would look like.
So what about the right part?
So the right part, we have something like one, two, three,
four, bunch of dots, and then we have two, and then a dot,
and then two.
What does that thing look like?
So that's going to look like four bases here in a stem.
Big loop, and then there's going to be two bases that
are paired, and then a bulge, and then
two more that are paired.
These things happen in real structures.
OK and then the arced notation is
a little more human-friendly.
It actually draws an arc between each pair
of bases that are hydrogen bonded.
So I'm sure you can imagine what those structures would
look like.
And it turns out that the arcs are very important.
Like whether those arcs cross each other
or not is sort of a fundamental classification of RNA
secondary structures, into the ones that are tractable
and the ones that are really difficult.
So pretty pictures of RNA.
So this is a lower resolution cryo-EM structure
of the bacterial ribosomes.
Remember, ribosomes have two sub-units-- a large sub-unit,
50S, and a small sub-unit, 30S.
And if you crack it open-- OK, so you basically split.
You sort of break the ribosome like that, and you look inside,
they're full of tRNAs.
So there are three pockets that are normally
distinguished within ribosomes.
The A site-- this is the site where
the tRNA enters that's going to add
a new amino acid to the growing peptide chain.
The P site, which is this tRNA will have it
[INAUDIBLE] with the actual growing peptide.
And then the exit tunnel where this tRNA will eventually--
the exit, the E site, which is the one that
was added a couple of residues ago.
So people often think of RNA structure
just in terms of these secondary structures
because they're much easier to generate
than tertiary structures, and they give you-- like for tRNA,
it gives you some pretty good information about how it works.
But for a large and complex structure like the ribosome,
it turns out that RNA is actually not
bad at building complex structures.
I would say it's not as good as protein,
but it is capable of constructing something
like a long tube.
And in fact, in the ribosome, you
find such a long tube right here.
That is where the peptide that's been synthesized
exits the ribosome.
And you'll notice it's not a large cavity in which
the protein might start folding.
It's a skinny tube that is thin enough that the polypeptide has
to remain linear, cannot start folding back on itself.
So you sort of extrude the protein
in a linear, unfolded confirmation,
and let it fold outside of the ribosome.
If it could fold inside that, that might clog it up.
That's probably one reason why it's not designed that way.
I'm sure that was tried bye evolution and rejected.
So if you look at the ribosome-- now remember,
the ribosome is composed of both RNA and protein--
you'll see that it's much more of one than the other.
And so it's really much more of the fettuccine, which
is the RNA part, than the linguini of the protein.
And if you also look at the distribution
of the proteins on the ribosome, you'll
see that they're not in the core.
They're kind of decorated around the edges.
It really looks like something that was originally made out
of RNA, and then you sort of added proteins as accessories
later.
And that's probably what happened.
This is based on the structures that
were solved a few years ago.
If you then look at where the nearest proteins are
to the active site-- actual catalytic site-- remember,
the ribosome catalyzes peptide in addition to an amino acid
to a growing peptide, so peptide bond formation--
you'll find that the nearest proteins are around
18 to 20 angstroms away.
And this is too far to do any chemistry,
so the active site residues or molecules
need to be within a few angstroms
to do any useful chemistry.
And so this basically proves that the ribosome.
Is a ribozyme.
That is, it's an RNA enzyme.
RNAs is [INAUDIBLE].
So here is the structure of a ribosome.
It's very kind of beautiful, and it's impressive
that somebody can actually solve the structure of something
this big.
But what is actually the practical use
of this structure?
Turns out there's quite an important practical application
of knowing the structure.
Any ideas?
AUDIENCE: Antibiotics.
PROFESSOR: Antibiotics.
Exactly.
So many antibiotics work by taking advantage of differences
between the prokaryotic ribosome structure
and eukaryotic ribosome structure.
So if you can make a small molecule--
these are some examples-- that will inhibit
prokaryotic ribosomes but hopefully not
inhibit eukaryotic ribosome, then
you can kill bacteria that might be infecting you.
So non-coding RNA.
So there's many different families of non-coding RNAs,
and I'm going to list some in a moment.
And I'm going to actually challenge
you, see if you can come up with any more
families of non-coding RNAs.
But they're receiving increasing interest,
I would say, ever since micro RNA's were discovered.
Sort of a boom in looking at different types
of non-coding RNAs.
Link RNA is also important and interesting, as well as many
of the classical RNA's like tRNAs and rRNAs and snoRNAs.
There may be new aspects of their regulation and function
that will be interesting.
And so when you're studying a non RNA,
it's very, very helpful to know its structure.
If it's going to base pair in trans with some other RNA--
as tRNAs do, as micro RNA's do, for example, or snRNAs
and snoRNAs-- then you want to know
which parts of the molecule are free
and which are internally based paired.
And if you want to predict non RNAs genes in a genome,
you may want to look for regions that
are under selection for conservation of RNA structure,
for conservation of the potential
to base pair at some distance.
If you see that, it's much more likely
that that region of the genome encodes a non-coding RNA
than it codes, for example-- there's a coding axon
or that it's a transcription factor binding
site or something like that that functions at the DNA level.
So having this notion of structure--
even just secondary structure-- is helpful for that application
as well, and predicting functions as well,
as I mentioned.
So co-variation.
So let's take a look at these sequences.
So imagine you've discovered a new class of mini micro RNA's.
They're only eight bases long, and you've sequence five
homologues from your five favorite mammals.
And these are the sequences that you get.
And you know that they're homologous
by [? a centimeter ?], they're in the same place
in the genome, and they seem to have the same function.
What could you say about their secondary structure
based on this multiple alignment?
You have to stare at it a little bit to see the pattern.
There's a pattern here.
Any ideas?
Anyone have a guess about what the structure is?
Yeah, go ahead.
AUDIENCE: There's a two base pair stem, and then
a four base loop.
PROFESSOR: Two base pair stem, four base loop,
and you have of the stem.
So how do you know that?
AUDIENCE: So if you look at the first two
and last two bases of each sequence,
the first and the eighths nucleotide
can pair with each other, and so can the second and the seventh.
PROFESSOR: Yeah.
Everyone see that?
So in the first column you have AUACG,
and that's complementary to UAUGC.
Each base is complementary.
And the second position is CAGGU complementary to GUCUA.
There's one slight exception there.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah.
Well, it turns out that that RNA-- although the Watson Crick
pairs GC and AU are the most stable-- GU pairs
are only a little bit less stable than AU pairs,
and they occur in natural RNA molecules.
So GU is allowed in RNA even though you would never
see that in DNA.
OK, so everyone see that?
So the structure is-- I think I have it here.
This would be co-variation You're changing the bases,
but preserving the ability to pair.
So when one base change-- when the first base changes from A
to U, the last base changes from U to A
in order to preserve that pairing.
You wouldn't know that if you just had two sequences,
but once you get several sequences,
it can be pretty compelling and allow
you to make a pretty strong inference that that
is the structure of that molecule.
So how would you do this?
So imagine you had a more realistic example where
you've got a non-coding RNA that's 100 or a few hundred
bases long, and you might have a multiple alignment of 50
homologous sequences.
You want something, you're not going
to be able to see it by eye.
You need sort of a more objective criterion.
So one method that's commonly used
is this statistic IX mutual information.
So if you look in your multiple alignment--
I'll just draw this here.
You have many sequences.
You consider every pair of columns--
this is a multiple alignment, so this column and this column--
and you calculate what we're going
to call-- what are we going to call it? f ix.
That would be the frequency of a nucleotide x.
You're in column i, so you just count how many A's, C's, G's,
and T's there are.
And similarly, f jy for all the possible values
of x and all the possible values of y.
So these are the base frequencies in each column.
And then you calculate the dinucleotide frequencies xy
at each pair of columns.
So in this colony, you say if there's an A here and a C here,
and then there's another AC down here,
and there's a total of one, two, three, four, five, six,
seven sequences, then f AC ij is 2/7.
So you just calculate the frequency of each dinucleotide.
These are no longer consecutive dinucleotides in a sequence
necessarily there.
They can be in arbitrary spacing.
OK, so you calculate those and then
you throw them into this formula,
and out comes a number.
So what does this formula remind of?
Have you seen a similar formula before?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Someone said [INAUDIBLE] Yeah, go ahead.
AUDIENCE: It reminds me of the Shannon entropy [INAUDIBLE].
PROFESSOR: Yeah, it looks like Shannon entropy,
but there's a log of a ratio in there,
so it's not exactly Shannon entropy.
So what other formula has a log of a ratio in it?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Relative.
So it actually looks like relative entropy.
So relative entropy of what versus what?
Who can sort of say more precisely if it's-- we'll say
it's relative entropy of something versus a p versus q.
And what is p and what is q?
Yeah, in the back.
AUDIENCE: Is it relative entropy of co-occurrence
versus independent occurrence?
PROFESSOR: Good.
Yeah. co-occurence-- everyone get that?
Co-occurrence of a pair of nucleotide xy at positions ij.
Versus q is an independent occurrence.
So if x and y occurred independently,
they would have this frequency.
So if you think about it, you calculate the frequency
of each base at each column in the multiple alignment.
And this is like your null hypothesis.
You're going to assume, what if they're evolving independently?
So if it's not a folded RNA-- or if it's a folded RNA
but those two columns don't happen to interact--
there's no reason to suspect that those bases would
have any relationship to each other.
So this is like your expected value
of the frequency of xy in position ij.
And then this p is your observed value.
So you're taking relative entropy of basically observed
over expected.
And so relative entropy has-- I haven't proved this,
but it's non-negative.
It can be 0, and then it goes up to some maximum,
a positive value, but it's never negative.
And what would it be if, in fact, p were equal to q?
What would this formula give?
This is where we're saying suppose.
Suppose this.
In general, this won't be sure, but suppose
it was equal to that.
We've got mi ij equals summation of what?
That log of this, which is equal to this,
so it's fx i fy j over the same thing--
hope you can see that-- log of-- log of 1 is 0, right?
So it's just 0.
So if the nucleotides of the two columns
occur completely independently, mutual information is 0.
And that's one reason it's called mutual information.
There's no information.
Knowing what's in column i gives you
no information about column j.
So remember, relative entities are measures of information,
not entropy.
And what is the maximum value that the mutual information
could have?
Any ideas on that?
Any guesses?
Joe, yeah.
AUDIENCE: You could have log base 2 log over f sub x,
f sub y.
PROFESSOR: Of 1?
OK, so you're saying if one of the particular dinucleotides
had a frequency of 1?
AUDIENCE: Yeah.
So if they're always the same whenever there's-- like an A,
there's always going to be a T.
PROFESSOR: Right.
So whenever there's an A, there's always a G or a T.
AUDIENCE: So then you'd get a 1 in the numerator,
and they're relative probably in the bottom, which
would be maximized if they were all even.
PROFESSOR: If they were all?
[INTERPOSING VOICES]
PROFESSOR: If they were uniform.
Yeah.
So did everyone get that?
So the maximum occurs if fx i and j-- they're both uniform,
so they're a quarter for every base at both positions.
That's the maximum entropy in the background distribution.
But then if fx y ij equals 1/4, for example, x equals y--
or in our case, we're not interested in that.
We're interested in x equals complement of y.
C of y is going to be the complement of y.
And 0 otherwise for x not equal complement of y.
OK, so for example, if we have only the dinucleotides AT,
CG, GC, and TA occur, and each of them
occurs with a frequency of 1/4, then
you'll have four terms in the sum because, remember,
the 0 log 0 is 0.
So you'll have four terms in the sum, and each of them
will look like 1/4 log 1/4 over a 1/4 times 1/4.
And so this will be 4, so log 2 of 4 4 is 2.
And so you have four terms that are each 1/4 times 2.
And so you'll get 2.
Well, this is not a sum.
These are the four terms.
These are the individual nonzero terms in that sum.
Does that make sense?
Everyone get this?
So that's why this is a useful measure of co-variation.
If what's in one column really strongly
influences what's in the other column,
and there's a lot of variation in the two columns,
and so you can really see that co-variation well,
then mutual information is maximized.
And that's basically what we just said,
is written down here.
So it's maximal.
They don't have to be complementary.
It would achieve this maximum of 2 if they are complementary,
but it would be also if they had some other very specific
relationship between the nucleotides.
So if you're going to use this, the way you would use it
is take your multiple alignment, calculate
the mutual information of each pair of columns--
so you actually have to make a table, i versus j,
all possible pairs of columns-- and then
you're going to look for the really high values.
And then when you find those high values, when
you look at what actual bases are tending to occur together,
you'll want to see that they're bases
that are complementary to one another.
And another thing that you'd want to see
is you'd want to see that consecutive positions in one
part of the alignment are co-varying
with consecutive positions in another part of the alignment
in the right way, in this sort of inverse complementary way
that RNA likes to pair.
Does that make sense?
So in a sort of nested way in your multiple alignment,
if you saw that this one co-varied with that,
and then you also saw that the next base co-varied
with the base right before this one,
and this one co-varies with that one,
that starts to look like a stem.
It's much more likely that you have a three-base stem
than that you just have some isolated base
pair out in the middle of nowhere.
It turns out it takes a few bases to make
a good thermodynamically stable stem,
and so you want to look for blocks of these things.
And so this works pretty well.
Yeah, actually, one point I want to make first
is that mutual information is nice
because it's kind of a useful concept
and it also relates to some of the entropy and relative
entropy that we've been talking about in the course before.
But it's not the only statistic that would work in practice.
You can use any measure of basically non-independence
between distributions.
A chi square statistic would probably
work equally well in practice.
And so here is a multiple alignment
of a bunch of sequences.
And what I've done is put boxes around columns
that have significant
mutual information with other sets of columns.
So for example, this set of columns here at the left-- the
far left-- has significant mutual information
with the ones at the far right.
And these ones, these four positions
co-vary with these four, and so forth.
So can you tell, based on looking
at this pattern of co-variation, what
the structure is going to be?
OK, let's say we start up here.
The first is going to pair with the last,
with something at the end.
Then we're going to have something
here in the middle that pairs with something else nearby.
Then we have something here that pairs
with something else nearby, then we have another like that.
Does that make sense?
So that there's these three pairs of columns
in the middle-- these two, these two, and these two-- and then
they're surrounded by this thing,
the first pairing with the last.
And so it's a clover leaf, so that's tRNA.
Yeah?
AUDIENCE: So with that previous slide, this table here,
you could create a co-variation matrix.
How would that-- or, and it could be--
PROFESSOR: How does that co-variations matrix--
how do you convert it to this representations?
AUDIENCE: I'm just wondering how this would go up.
Like let's say you took the co-variation matrix--
PROFESSOR: Oh, what would it look like?
AUDIENCE: --and visualized it as a heat map--
PROFESSOR: In the co-variation matrix.
AUDIENCE: Yeah.
What would it look like in this particular example?
PROFESSOR: Yeah, that's a good question.
OK, let's do that.
I haven't thought about that before,
so you'll have to help me on this.
So here's the beginning.
We're going to write the sequence from 1
to n in both dimensions.
And so here's the beginning, and it co-varies with the end.
So this first would have a co-variation with the last,
and then the second would co-vary with the second
to last, and so forth.
So you get a little diagonal down here.
That's this top stem here.
And then what about the second stem?
So then you have something down here
that's going to co-vary with something kind of near by it.
So block two is going to co-vary with block three.
And again, it's going to be this inverse complementary kind
of thing like that.
It's symmetrical, so you get this with that.
But you only have to do one half,
so you can just do this upper half here.
So you get that.
So it would look something like that.
AUDIENCE: So with the diagonal line orthogonal
to the diagonal of the matrix--
PROFESSOR: Yeah, that's because they're inverse complementary.
AUDIENCE: OK.
PROFESSOR: That make sense?
Good question.
But we'll see an example like that later
actually, as it turns out.
All right, so here's my question for you.
You're setting this non-coding RNA.
It has some length.
You have some number of sequences.
They might have some structure.
Is this method going to work for you, or is it not?
What is required for it to work?
For example, would I want to isolate
this gene-- this non-coding RNA gene--
just from primates, from like human, gorilla,
chimp, orangutan, and do that alignment?
Or would I want to go further?
Would I want to go back to the rodents and dog, horse--
how far do you want to go?
Yeah, question.
AUDIENCE: I think we a need a very strong sequence alignment
for this, so we cannot go very far,
because if you don't have a high percentage homology,
then you will see all sorts of false positives.
PROFESSOR: Absolutely.
So if you go too far, your alignment will suffer,
and you need an alignment in order
to identify the corresponding columns.
So that puts an upper limit on how far you can go.
But excellent point.
Is there a lower limit?
Do you want to go as close as possible,
like this example I gave with human, chimp, orangutan?
Or is that too close?
Why is too close bad?
Tim?
AUDIENCE: Maybe if you're too close,
then the sequence is having to [INAUDIBLE]
to give you enough information [INAUDIBLE].
PROFESSOR: Yeah, exactly.
They're all the same.
Actually, you'll get 1 times 1 over 1
in that mutual information statistic, which log of that
is going to be 0.
There's zero mutual information if they're all the same.
So there has to be some variation,
and the structure has to be conserved.
That's key.
You have to assume that the structure is well conserved
and you have to have a good alignment
and there has to be some variation,
a certain amount of variation.
Those are basically the three keys.
Secondary structure has a more highly conserved sequence.
Sufficient divergence so that you have these variations,
and sufficient number of homologues you have to get good
statistics, and not so far they your alignment is bad.
Sorry about that.
Sally?
AUDIENCE: It seems like another thing
that we assume here is that you can project it onto a plane
and it will lie flat.
So if you have some very important, weird folding
that allows you to, say, crisscross the rainbow thing.
PROFESSOR: Yeah, crisscross the rainbow.
Yeah, very good question.
So in the example of tRNA, if you
were to do that arc diagram for tRNA,
it would look like another big arc--
that's the first and last-- and then you
have these three nested arcs.
Nothing crisscrossing.
What if I saw-- [INAUDIBLE]-- two blocks of sequence that
have a relationship like that?
Is that OK?
With this method, the co-variation, that's OK.
There's no problem there.
What does this structure look like?
So [INAUDIBLE] you have a stem, then you have a loop,
and then a stem.
So this is 1 pairs with 3.
That's 1.
That's 3.
Then you've got 2 up here, but 2 pairs with 4.
So here's 4 over here, so 4 is going
to have to come back up here and pair with 2.
This is 2 over here.
So that is called a pseudoknot.
It's not really a knot because this thing doesn't
go through the loop, but it kind of
behaves like a knot in some ways.
And so do these actually occur in natural RNAs?
Yes, Tim is nodding.
And are they important?
Can you give me an example where they are important
biologically?
AUDIENCE: [INAUDIBLE]
[INTERPOSING VOICES]
PROFESSOR: Riboswitches.
We're going to come to what riboswitches
are in a moment for those not familiar.
And I think I have an example later
of a pseudoknot that's important.
So that's a good question.
I think I should have added to this list the point
that you made in the back that they
have to be close enough that you can get a good alignment.
I should add that to this last.
Thanks.
It's a good point.
All right, so classes of non-coding RNAs.
As promised, my favorites listed here.
Everyone knows tRNAs, rRNAs.
You can think of UTRs as being non RNAs.
They often have structure that can
be involved in regulating the message.
snRNAs involved splicing.
snoRNAs-- small nucleolar RNAs-- are
involved in directing modification
of other RNAs, such as ribosomal RNAs and snRNAs, for example.
Terminators of transcription in prokaryotes
are like little stem loop structures.
RNaseP is an important enzyme.
SRP is involved in targeting proteins with signal peptides
to the export machinery.
We won't go into tmRNA.
micro RNAs and link RNAs, you probably
know, and riboswitches.
So Tim, can you tell us what a riboswitch is?
AUDIENCE: A riboswitch is any RNA structure
that changes confirmation according
to some stimulus [INAUDIBLE] or something in the cell.
It could be an ion, critical changes in the structure.
[INAUDIBLE]
PROFESSOR: Yeah, that was great.
So just for those that may not have heard,
I'll just say it again.
So a riboswitch is any RNA that can
have multiple confirmations, and changes confirmation
in response to some stimulus-- temperature, binding
of some ligand, small molecules, something like that, et cetera.
And often, one of those structures
will block a particular regulatory element.
I'll show an example in a moment.
And so when it's in one confirmation,
the gene will be repressed.
And when it's in the other, it'll
be on. so it's a way of using RNA's secondary structure
to sense what's going on in the cell
and to appropriately regulate gene expression.
All right, so now we're going to talk about a second approach.
So this would be the approach.
You've got some RNA.
It may not do something, and maybe you can't find any
homologues.
It might be some newly evolved species-specific RNA,
or your studying some obscure species
where you don't have a lot of genomic sequence around.
So you want to use the first principles, approach,
the energy minimization approach.
Or maybe you have the homologues,
but you don't trust your alignment.
You want a second opinion on what
the structure is going to be.
So just in the way that protein folding--
you could think of an equilibrium model
where it's determined by folding free energy,
and enthalpy will favor base pairing.
You get gain some enthalpy when you form a hydrogen bond,
and entropy will tend to favor unfolding.
So an RNA molecule that's linear has
all this confirmational flexibility,
and lose some of that when you form a stem.
It forms a helix.
Those things don't have as much flexibility.
And even the nucleotides in the loop are a little bit
confirmationally-- they're not as flexible
as they were when it was linear.
So that means that at high temperatures,
it'll favor unfolding.
So the earliest approaches were approaches
that sought to maximize the number of base pairs.
So they basically ignore entropy and focus on the enthalpy
that you gain from forming base pairs.
And so Ruth Nussinov described the first algorithm
to figure out what is the maximum number of base pairs
that you can form in an RNA.
And so a way to think about this is
imagine you've got this sequence.
What is the largest number of base pairs
I can form with this sequence?
I could just draw all possible base pairs.
That A can pair with that T. This A can pair with that T.
They can't both pair simultaneously, right?
And this C can pair with that G. So if we don't allow crossing,
which-- coming back to Sally's point--
this would cross this, right?
So we're not going to allow that.
So the best you could do be to have this A pair with this C
and this C pair with this G and form this little structure.
This is not realistic because RNA loops can't be one base.
They minimum is about three.
But just for the sake of argument,
you can list all these out, but imagine now
you've got 100 bases here.
Every base will on average potentially
be able to pair with 24 or 25 other bases.
So you're just going to have just an incredible mishmash
of possible lines all crisscrossing.
So how do you figure out how to maximize that pairing?
Any ideas?
Don, yeah?
AUDIENCE: You look for sections of homology.
PROFESSOR: We're not using homology.
We're doing [INAUDIBLE]
AUDIENCE: I'm sorry, not homology, but sections where--
PROFESSOR: Complementary?
AUDIENCE: Complementary.
Yeah, that's the word I was thinking.
PROFESSOR: The blocks are complementary.
AUDIENCE: And then so--
PROFESSOR: You could blast the sequence against inverse
complements itself and look for little blocks.
You could do that.
That's not what people generally do,
mostly because the blocks of complementarity in real RNA
structures are really short.
They can be two, three, four, bases.
Sally, yeah?
AUDIENCE: Could you use [INAUDIBLE] approach
where you just start with a very small case and build up?
PROFESSOR: So we've seen that work for protein sequence
alignment.
We've seen it work for the Viterbi algorithm.
So that is sort of the go-to approach in bioinfomatics,
is to use some sort of dynamic programming.
Now this one for RNA secondary structure
that Nussinov came up with is a little bit
different than the others.
So you'll see it has a kind of different flavor.
It turns out to be actually it's a little hard
to get your head around at the beginning,
but it's actually easier to do by hand.
So let's take a look at that.
OK, so recursive maximization of base pairing.
Now the thing about base pairing that's
different from these other problems
is that the first base in the sequence
can base pair with the last.
How do you chop up a sequence?
Remember with Needleman-Wunsch and with Viterbi
we go from the beginning to the end,
and that's a logical order.
But with base pairing, that's actually not a logical order.
You can't really do it that way.
So instead, you go from the inside out.
You start in the middle of a sequence
and work your way outwards in both directions.
Or another way to think about it is you start with you write
the sequence from 1 to n on both axes,
and then actually we'll see that we initiate the diagonal all
to 0's.
And then we think about these positions here next.
So 1 versus 2.
Could 1 pair with 2?
And could 2 pair with 3?
Those are like little bits of possible RNA
secondary structure.
Again, we're ignoring this fact that loops
have to be certain minimum.
This is sort of a simplified case.
And then you build outwards.
So you conclude that base 4 here could pair with base 5,
so we're going to put a 1 there.
And then we're going to build outward
from that toward the beginning of the sequence
and toward the end, adding additional base pairs
when we can.
That's basically the way the [INAUDIBLE] works.
And so that's one key idea, that we
go from sort of close sequences, work
outward, to faraway sequences.
And the second key idea is that the relationship
that, as you add more bases on the outside of what you've
already got, that the optimal structure in that larger
portion of sequence space is related
to the optimal structures of smaller portions of it
in one of four different ways.
And these are the four ways.
So let's look at these.
So the first one is probably the simplest
where if you're doing this, you're here somewhere,
meaning you've compared sequences from position,
let's say, i minus 1 to j minus 1 here.
And then we're going to consider adding-- actually,
it depends how you number your sequence.
Let me see how this is done.
Sorry. i plus 1.
i plus 1 to j minus 1.
We figured out what the optimal structure is in here,
let's suppose.
And now we're going to consider adding one more
base on either end.
We're going to add j down here, and we're
going to ask if it pairs with i.
And if so, we're going to take whatever the optimal structure
was in here and we're going to add one base pair,
and we're going to add plus 1 because now
it's got one additional.
We're counting base pairs.
So that's that first case there.
And then the second case is you could also consider just
adding one unpaired base onto whatever structure you had,
and then you don't add one.
And you could go in either direction.
You can go sort of toward of the beginning of the sequence
or toward the end of the sequence.
And then the third one is the tricky one,
is what's called a bifurcation.
You could consider that actually i and j
are both paired, but not with each other.
That i pairs with something that was inside here
and j pairs with something that was inside here.
So your optimal parse from i to j, if you will,
is not going to come from the optimal parse from i plus 1
to j minus 1.
It's going to come from rethinking this and doing
the optimal parse from here to here and from here to here,
and combining those two.
So you're probably confused by now,
so let me try to do an example.
And then I have an analogy that will confuse you further.
So ask me for that one.
This was the simplest one I could come up
with that has this property.
OK, so we said before that if you were doing the optimal
from 1 to 5, that it would be the AC pairing with the GT.
We do that one.
And now if you notice, this guy is kind of a similar sequence.
I just added a T at the beginning and an A at the end.
And so you can probably imagine that the best structure of this
is here, those three.
You've got three pairs of this sub-sequence here.
That's as good as you can do with seven bases.
You can only get three pairs.
And this is as good as you can do with five,
so these are clearly optimal.
So the issue comes that if you're starting from somewhere
in the middle here-- let's say you are-- let's see,
so how would you be doing this?
You start here.
Let's suppose the first two you consider are these two.
You consider pairing that T with that A.
You can see this is not going to go well.
You might end up with that as your optimal substructure
of this region.
Remember, you're working from the inside out,
so you're going from here to here, and you end up with that.
And what do you do here?
You don't have a G to pair the C to,
so you add another unpaired base.
Now you've got this optimal substructure
of a sequence that's almost the whole sequence.
It's just missing the first and last bases,
but it only has three base pairs.
So when you go to add this, you can say,
oh, I can't add any more base pairs, so I've only got three.
But you should consider that we've already
solved the optimal structure of that,
and we had two nice pairs here.
We had that pair and that pair, and we already
solved the substructure of the optimal structure
of this portion here, and you had those three pairs.
And so you can combine those two and all of a sudden
you can do much better.
So that's what that bifurcation thing is about.
So this is the recursion working out,
and you can see that's the base pairing one.
You can add one, or you can just add an unpaired base
and you don't add anything.
Or you consider all the possible locations
of bifurcations in-between the two positions you're adding,
i and j, and you consider all the possible pairs.
And you just sum up each pair and go-- I'm sorry,
you don't sum them up.
You consider them all, and then you take the maximum.
All right, so the algorithm is to take an n by n matrix,
initialize the diagonal to 0, and initialize the sub-diagonal
to 0 also.
Just don't think too much about that.
Just do it.
And then fill in this matrix recursively
from the diagonal up and to the right.
And it actually doesn't matter what order you fill it in
as long as you're kind of working your way up
into the right.
You have to have the thing to the left and the thing
below already filled in if you're going to fill in a box.
And then you keep track of the optimal score, which
is going to be the sum of base pairs.
And then you also keep track of how you got there.
What base pair did you add so that you can trace back?
And then when you get up to the upper right
corner of this matrix, you then trace back.
So here is a partially filled in this matrix.
This is from that the Nature Biotechnology Review.
And the 0's are filled in.
So here's what I want you to do at home,
is print out, photocopy or whatever-- make this matrix,
or make a bigger version of it perhaps--
and look at the sequence and fill in this matrix,
and fill in the little arrows every time you add a base pair.
It's actually not that hard.
There are no bifurcations in this, so that's the tricky one.
Ignore that one.
You'll just be adding base pairs.
It'll be pretty easy.
And then you can reconstruct the sequence.
So here it is filled in.
And the answer is given, so you can check yourself.
But do it without looking at the answer.
And then you go to the upper right corner.
That means that the optimal structure
from the beginning of the sequence to the end-- which,
of course, was our goal all along.
And then you trace back and you can
see whenever you're moving diagonally here,
you're adding a base pair.
Remember, you add one on each end,
and so you're moving diagonally and adding the base pair,
and you get this little structure here.
So computational complexity of the algorithm.
You could think about this but I'll just tell you.
It's memory n squared because you've
got to fill in this matrix, so square
of the length of the sequence.
Time n cubed.
This is bad now.
And why is it n cubed?
It's n cubed because you have to fill in a matrix that's n by n.
And then when you do that maximization step,
that check for bifurcations, that's sort of of order n,
as well.
So n cubed-- so this means that RNA folding is slow.
And in fact, some of the servers won't
allow you to fold anything more than a thousand
bases because they'll take forever or something like that.
And it cannot handle pseudoknots.
If you think through the recursion,
pseudoknots will be a problem.
I'm going to just show you-- yeah,
I'll get to this-- that these are from the viruses.
Real viruses, some of them have pseudoknots
like these ones shown here, and some even
have these kissing loops, which is another type where
the two stem loops, the loops interact.
And the pseudoknots in particular
are important in the viral life cycle.
They can actually cause programmed ribosomal frame
shifting.
When the ribosomes hits one of the things,
normally it just denatures RNA secondary structure.
When it hits a pseudoknot, it'll actually
get knocked back by one and will start
translating in a different frame.
And that's actually useful to the virus
to do that under certain circumstances.
That's how HIV makes the replicated polymerase,
is by doing a frame shift on the ribosome using a pseudoknot.
So these things are important.
And there's fancier methods that use
more sophisticated thermodynamic models where
GC counts more than AU.
And I won't go into the details, but I just
wanted to show you some pretty pictures here
that the Zuker algorithm-- this is
a real world RNA folding algorithm-- calculates not only
the minimum energy fold, but also sub-optimal folds,
and the probabilities of particular base pairs,
summing over all the possible structures that RNA could form,
weighted by their free energy.
So it's the full partition function.
It's not perfectly accurate.
It gets about 70% of base pairs correct,
which means it usually gets things right,
but occasionally totally wrong.
And there's a website for the Mfold server, which is actually
one of the most beautiful websites in bioinfomatics,
I would say.
And also if you want to run it locally,
you should download the Vienna RNAfold package,
which has a very similar algorithm.
And I just wanted to show you one or two examples.
So this is the U5 snRNA.
This is the output of Mfold.
It predicts this structure.
And then this what's called the energy dot plot, which
shows the bases in the optimal structure down below here
and then sort of these suboptimal structures here.
And you can see there's no ambiguity.
It's totally confident in this structure.
Then I ran the lysine riboswitch through this program,
and I got this.
I got the minimum for energy structure
down in the lower left.
And then you see there's a lot of other colored dots.
Those are from the suboptimal structures.
So it looks like this thing has multiple structures, which
of course it does.
So the way that this one works is, in the absence of lysine,
it forms this structure where the ribosome binding
sequences-- this is prokaryotic-- is exposed.
And so the ribosome can enter and translate
these lysine biosynthetic enzymes.
But then when lysine accumulates to a certain level,
it can interact with the RNA and shift it's structure
so that you now form this stem, which
sequesters the ribosome binding sequence
and blocks lysine biosynthesis.
So a very clever system.
And it turns out that there's dozens
of these things in bacterial genomes,
and they control a lot of metabolism.
So they're very important.
And there may be some in eukaryotes, too,
and that would be good.
If anyone's looking for a product, not happy
with their current project, you might
think about looking for more riboswitches.
So I'm going to have to end there.
And thank you guys for your attention,
and good luck on the midterm.