Many interesting RNAs conserve a secondary structure of base-pairing interactions more than they conserve their sequence. This makes RNA sequence analysis more complicated and difficult than protein or DNA sequence analysis. RNA secondary structure problems are a natural application for probabilistic models based on the stochastic context-free grammars introduced in Chapter 9. In this chapter, we will examine two RNA analysis problems of biological interest.

The first problem is RNA secondary structure prediction for a single sequence. We will outline two well-known dynamic programming algorithms for RNA secondary structure prediction, the Nussinov and the Zuker algorithms. Then we will use RNA secondary structure prediction as an introductory example for the use of SCFGs for RNA analysis, by developing a small SCFG that implements a probabilistic version of the Nussinov algorithm.

The second is a related set of problems, having to do with the analysis of multiple alignments of families of related RNAs. Like Chapter 5, where profile HMMs were used for both multiple alignment and for database searching, we develop RNA structure profiles called ‘covariance models’ (CMs) for dealing with RNA multiple alignments with secondary structure constraints included. Covariance models are used for both RNA multiple alignment and database searches. Consensus structure prediction from RNA multiple alignments, a process called comparative RNA sequence analysis, is also somewhat automated by RNA covariance model training algorithms.

As you read this chapter, bear in mind that SCFG-based RNA analysis methods are not widely known or used.

So far we have concentrated on the intrinsic properties of single sequences, such as CpG islands in DNA, or on pairwise alignment of sequences. However, functional biological sequences typically come in families, and many of the most powerful sequence analysis methods are based on identifying the relationship of an individual sequence to a sequence family. Sequences in a family will have diverged from each other in their primary sequence during evolution, having separated either by a duplication in the genome, or by speciation giving rise to corresponding sequences in related organisms. In either case they normally maintain the same or a related function. Therefore, identifying that a sequence belongs to a family, and aligning it to the other members, often allows inferences about its function.

If you already have a set of sequences belonging to a family, you can perform a database search for more members using pairwise alignment with one of the known family members as the query sequence. To be more thorough, you could even search with all the known members one by one. However, pairwise searching with any one of the members may not find sequences distantly related to the ones you have already. An alternative approach is to use statistical features of the whole set of sequences in the search. Similarly, even when family membership is clear, accurate alignment can be often be improved significantly by concentrating on features that are conserved in the whole family.

Until now, we have treated biological sequences as one-dimensional strings of independent, uncorrelated symbols. This assumption is computationally convenient but not structurally realistic. The three-dimensional folding of proteins and nucleic acids involves extensive physical interactions between residues that are not adjacent in primary sequence. Can probabilistic models of proteins and nucleic acid sequences be developed that allow for longer range interactions? Can we compute efficiently with such models? In this chapter, we will step back from models of particular sequence problems and address these more theoretical issues. We will see how many of the methods described in previous chapters fit into a more general view of modelling sequences

A general theory for modelling strings of symbols has been developed by computational linguists [Chomsky 1956; 1959]. This theory is known as the Chomsky hierarchy of transformational grammars. In the Chomsky hierarchy, most of the models we have used so far in this book are the lowest of four types of model of increasing complexity and descriptive power. Transformational grammars were developed in an attempt to understand the structure of natural languages. They became important in theoretical computer science [Hopcroft & Ullman 1979; Gersting 1993] because computer languages, unlike natural languages, can be precisely specified as formal grammars. Recently, transformational grammars have been applied to sequence analysis problems in molecular biology [Searls 1992; Dong & Searls 1994; Rosenblueth et al. 1996].

Astronomy began when the Babylonians mapped the heavens. Our descendants will certainly not say that biology began with today's genome projects, but they may well recognise that a great acceleration in the accumulation of biological knowledge began in our era. To make sense of this knowledge is a challenge, and will require increased understanding of the biology of cells and organisms. But part of the challenge is simply to organise, classify and parse the immense richness of sequence data. This is more than an abstract task of string parsing, for behind the string of bases or amino acids is the whole complexity of molecular biology. This book is about methods which are in principle capable of capturing some of this complexity, by integrating diverse sources of biological information into clean, general, and tractable probabilistic models for sequence analysis.

Though this book is about computational biology, let us be clear about one thing from the start: the most reliable way to determine a biological molecule's structure or function is by direct experimentation. However, it is far easier to obtain the DNA sequence of the gene corresponding to an RNA or protein than it is to experimentally determine its function or its structure. This provides strong motivation for developing computational methods that can infer biological information from sequence alone. Computational methods have become especially important since the advent of genome projects. The Human Genome Project alone will give us the raw sequences of an estimated 70,000 to 100,000 human genes, only a small fraction of which have been studied experimentally.

Having introduced some methods for pairwise alignment in Chapter 2, the emphasis will switch in this chapter to questions about a single sequence. The main aim of the chapter is to develop the theory for a very general form of probabilistic model for sequences of symbols, called a hidden Markov model (abbreviated HMM). The types of question we can use HMMs and their simpler cousins, Markov models, to consider are: ‘Does this sequence belong to a particular family?’ or ‘Assuming the sequence does come from some family, what can we say about its internal structure?’ An example of the second type of problem would be to try to identify alpha helix or beta sheet regions in a protein sequence.

As well as giving examples from the biological sequence world, we also give the mathematics and algorithms for many of the operations on HMMs in a more general form. These methods, or close analogues of them, are applied in many other sections of the book. This chapter therefore contains a fairly large amount of mathematically technical material. We have tried to organise it so that the first half, approximately, leads the reader through the essential algorithms using a single biological example. In the later sections we introduce a variety of other examples to illustrate more complex extensions of the basic approaches.

In the next chapter, we will see how HMMs can also be applied to the types of alignment problem discussed in Chapter 2, in Chapter 5 they are applied to searching databases for protein families, and in Chapter 6 to alignment of several sequences simultaneously.

Introduction

Our goal in this chapter is to formulate probabilistic models for phylogeny and show how trees can be inferred from sets of sequences, either by maximum likelihood or by sampling methods. We also review the phylogenetic methods of the previous chapter, and show that they often have probabilistic interpretations, though they are not usually presented this way.

Overview of the probabilistic approach to phylogeny

The basic aim of probability-based phylogeny is to rank trees either according to their likelihood P(data∣tree), or, if we are taking a more Bayesian view, according to their posterior probability P (tree∣data). There may be subsidiary aims, such as finding the likelihood or posterior probability of some particular taxonomic feature, such as a grouping of a set of organisms on a single branch. To achieve any of these aims, we must be able to define and compute P(x•∣T, t•), the probability of a set of data given a tree. Here the data are a set of n sequences xj for j = 1…n, which we write compactly as x•. T is a tree with n leaves with sequence j at leaf j, and the t• are the edge lengths of the tree. To define P(x•∣T, t•) we need a model of evolution, i.e. of the mutation and selection events that change sequences along the edges of a tree.

At a Snowbird conference on neural nets in 1992, David Haussler and his colleagues at UC Santa Cruz (including one of us, AK) described preliminary results on modelling protein sequence multiple alignments with probabilistic models called ‘hidden Markov models’ (HMMs). Copies of their technical report were widely circulated. Some of them found their way to the MRC Laboratory of Molecular Biology in Cambridge, where RD and GJM were just switching research interests from neural modelling to computational genome sequence analysis, and where SRE had arrived as a new postdoctoral student with a background in experimental molecular genetics and an interest in computational analysis. AK later also came to Cambridge for a year.

All of us quickly adopted the ideas of probabilistic modelling. We were persuaded that hidden Markov models and their stochastic grammar analogues are beautiful mathematical objects, well fitted to capturing the information buried in biological sequences. The Santa Cruz group and the Cambridge group independently developed two freely available HMM software packages for sequence analysis, and independently extended HMM methods to stochastic context-free grammar analysis of RNA secondary structures. Another group led by Pierre Baldi at JPL/Caltech was also inspired by the work presented at the Snowbird conference to work on HMM-based approaches at about the same time.

By late 1995, we thought that we had acquired a reasonable amount of experience in probabilistic modelling techniques. On the other hand, we also felt that relatively little of the work had been communicated effectively to the community.

In the previous chapter, we considered the problem of multiple alignment of sets of sequences. One can argue [Sankoff, Morel & Cedergren 1973] that alignment of sequences should take account of their evolutionary relationship. For example, an alignment that implies many substitutions between closely related sequences is less plausible than one that makes most of its changes over large evolutionary distances.

Some multiple alignment algorithms use a tree; for instance, we have seen that several progressive alignment algorithms use a ‘guide tree’. As the name suggests, this tree is meant to guide the clustering process rather than satisfy a taxonomist. In this chapter we shift emphasis, and begin to take a serious interest in building trees. However, we do not lose sight of alignment: the last section describes methods for simultaneous alignment and tree building.

We concentrate here on two general approaches to tree building: distance methods and parsimony; the next chapter formulates phylogeny probabilistically.

The tree of life

The similarity of molecular mechanisms of the organisms that have been studied strongly suggests that all organisms on Earth had a common ancestor. Thus any set of species is related, and this relationship is called a phylogeny. Usually the relationship can be represented by a phylogenetic tree. The task of phylogenetics is to infer this tree from observations upon the existing organisms.

Traditionally, morphological characters (both from living and fossilised organisms) have been used for inferring phylogenies.

In Chapter 5, we assumed that a reasonable multiple sequence alignment was already known and provided the starting point for constructing a profile HMM. We now look at what a ‘reasonable’ multiple alignment is, and at ways to construct one automatically from unaligned sequences.

Multiple alignments must usually be inferred from primary sequences alone. Biologists produce high quality multiple sequence alignments by hand using expert knowledge of protein sequence evolution. This knowledge comes from experience. Important factors include: specific sorts of columns in alignments, such as highly conserved residues or buried hydrophobic residues; the influence of secondary and tertiary structure, such as the alternation of hydrophobic and hydrophilic columns in exposed beta sheet; and expected patterns of insertions and deletions, that tend to alternate with blocks of conserved sequence. Furthermore, the phylogenetic relationships between sequences dictate constraints on the changes that occur in columns and in the patterns of gaps. RNA alignments involve similar knowledge but additionally they are often strongly constrained by a secondary structure model that in many cases has also been inferred from primary sequence data (Chapter 10).

Manual multiple alignment is tedious. Automatic multiple sequence alignment methods are a topic of extensive research in computational biology. In general, an automatic method must have a way to assign a score so that better multiple alignments get better scores. We should carefully distinguish the problem of scoring a multiple alignment from the problem of searching over possible multiple alignments to find the best one.