This chapter describes genome sequence data and explains the relevance of the statistics, computation and algebra that we have discussed in Chapters 1–3 to understanding the function of genomes and their evolution. It sets the stage for the studies in biological sequence analysis in some of the later chapters.

Given that quantitative methods play an increasingly important role in many different aspects of biology, the question arises: why the emphasis on genome sequences? The most significant answer is that genomes are fundamental objects that carry instructions for the self-assembly of living organisms. Ultimately, our understanding of human biology will be based on an understanding of the organization and function of our genome. Another reason to focus on genomes is the abundance of high fidelity data. Current finished genome sequences have less than one error in 10,000 bases. Statistical methods can therefore be directly applied to modeling the random evolution of genomes and to making inferences about the structure and organization of functional elements; there is no need to worry about extracting signal from noisy data. Furthermore, it is possible to validate findings with laboratory experiments.

The rate of accumulation of genome sequence data has been extraordinary, far outpacing Moore's law for the increasing density of transistors on circuit chips. This is due to breakthroughs in sequencing technologies and radical advances in automation. Since the first completion of the genome of a free living organism in 1995 (Haemophilus Influenza [Fleischmann et al., 1995]), biologists have completely sequenced over 200 microbial genomes, and dozens of complete invertebrate and vertebrate genomes.

Some of the statistical models introduced in Chapter 1 have the feature that, aside from the observed data, there is hidden information that cannot be determined from an observation. In this chapter we consider graphical models with hidden variables, such as the hidden Markov model and the hidden tree model. A natural problem in such models is to determine, given a particular observation, what is the most likely hidden data (which is called the explanation) for that observation. This problem is called MAP inference (Remark 4.13). Any fixed values of the parameters determine a way to assign an explanation to each possible observation. A map obtained in this way is called an inference function.

Examples of inference functions include gene-finding functions which were discussed in [Pachter and Sturmfels, 2005, Section 5]. These inference functions of a hidden Markov model are used to identify gene structures in DNA sequences (see Section 4.4. An observation in such a model is a sequence over the alphabet Σ′ = {A, C, G, T}.

After a short introduction to inference functions, we present the main result of this chapter in Section 9.2. We call it the Few Inference Functions Theorem, and it states that in any graphical model the number of inference functions grows polynomially if the number of parameters is fixed. This theorem shows that most functions from the set of observations to possible values of the hidden data cannot be inference functions for any choice of the model parameters.

We present a new, statistically consistent algorithm for phylogenetic tree construction that uses the algebraic theory of statistical models (as developed in Chapters 1 and 3). Our basic tool is Singular Value Decomposition (SVD) from numerical linear algebra.

Starting with an alignment of n DNA sequences, we show that SVD allows us to quickly decide whether a split of the taxa occurs in their phylogenetic tree, assuming only that evolution follows a tree Markov model. Using this fact, we have developed an algorithm to construct a phylogenetic tree by computing only O(n2) SVDs.

We have implemented this algorithm using the SVDLIBC library (available at http://tedlab.mit.edu/~dr/SVDLIBC/) and have done extensive testing with simulated and real data. The algorithm is fast in practice on trees with 20–30 taxa.

We begin by describing the general Markov model and then show how to flatten the joint probability distribution along a partition of the leaves in the tree. We give rank conditions for the resulting matrix; most notably, we give a set of new rank conditions that are satisfied by non-splits in the tree. Armed with these rank conditions, we present the tree-building algorithm, using SVD to calculate how close a matrix is to a certain rank. Finally, we give experimental results on the behavior of the algorithm with both simulated and real-life (ENCODE) data.

The general Markov model

We assume that evolution follows a tree Markov model, as introduced in Section 1.4, with evolution acting independently at different sites of the genome.

Graphical models are powerful statistical tools that have been applied to a wide variety of problems in computational biology: sequence alignment, ancestral genome reconstruction, etc. A graphical model consists of a graph whose vertices have associated random variables representing biological objects, such as entries in a DNA sequence, and whose edges have associated parameters that model transition or dependence relations between the random variables at the nodes. In many cases we will know the contents of only a subset of the model vertices, the observed random variables, and nothing about the contents of the remaining ones, the hidden random variables. A common example is a phylogenetic tree on a set of current species with given DNA sequences, but with no information about the DNA of their extinct ancestors. The task of finding the most likely set of values of the hidden random variables (also known as the explanation) given the set of observed random variables and the model parameters, is known as inference in graphical models.

Clearly, inference drawn about the hidden data is highly dependent on the topology and parameters (transition probabilities) of the graphical model. The topology of the model will be determined by the biological process being modeled, while the assumptions one can make about the nature of evolution, site mutation and other biological phenomena, allow us to restrict the space of possible transition probabilities to certain parameterized families. This raises several questions.

Using homologous sequences from eight vertebrates, we present a concrete example of the estimation of mutation rates in the models of evolution introduced in Chapter 4. We detail the process of data selection from a multiple alignment of the ENCODE regions, and compare rate estimates for each of the models in the Felsenstein hierarchy of Figure 4.7. We also address a standing problem in vertebrate evolution, namely the resolution of the phylogeny of the Eutherian orders, and discuss several challenges of molecular sequence analysis in inferring the phylogeny of this subclass. In particular, we consider the question of the position of the rodents relative to the primates, carnivores and artiodactyls; we affectionately dub this question the rodent problem.

Estimating mutation rates

Given an alignment of sequence homologs from various taxa, and an evolutionary model from Section 4.5, we are naturally led to ask the question, “what tree (with what branch lengths) and what values of the parameters in the rate matrix for that model are suggested by the alignment?” One answer to this question, the so-called maximum-likelihood solution, is, “the tree and rate parameters which maximize the probability that the given alignment would be generated by the given model.” (See also Sections 1.3 and 3.3.)

There are a number of available software packages which attempt to find, to varying degrees, this maximum-likelihood solution. For example, for a few of the most restrictive models in the Felsenstein hierarchy, the package PHYLIP [Felsenstein, 2004] will very efficiently search the tree space for the maximum-likelihood tree and rate parameters.

INTRODUCTION

An oligonucleotide probe is a short piece of single-stranded DNA complementary to the target gene whose expression is measured on the microarray by that probe. In most microarray applications, oligonucleotide probes are between 20 and 60 bases long. The probes are either spotted onto the array or synthesised in situ, depending on the microarray platform (Chapter 1).

Usually, oligonucleotide probes for microarrays are designed within several hundred bases of the 3′ end of the target gene sequence. So for a fixed oligonucleotide length, there are several hundred potential oligonucleotides, one for each possible starting base. Some of these oligonucleotides work better than others as probes on a microarray. This chapter describes methods for the computer selection of good oligonucleotide probes.

What Makes a Good Oligonucleotide Probe?

Good oligonucleotide probes have three properties: they are sensitive, specific and isothermal.

Asensitive probe is one that returns a strong signal when the complementary target is present in the sample. There are two factors that determine the sensitivity of a probe:

1. ▪ The probe does not have internal secondary structure or bind to other identical probes on the array.

2. ▪ The probe is able to access its complementary sequence in the target, which could potentially be unavailable as a result of secondary structure in the target.

A specific probe is one that returns a weak signal when the complementary target is absent from the sample; i.e., it does not cross-hybridise. There are two factors that determine the specificity of a probe:

1. ▪ Cross-hybridization to other targets as a result of Watson–Crick base-pairing

2. ▪ Non-specific binding to the probe; e.g., as a result of G-quartets

DNA array technology is almost fifteen years old, and still rapidly evolving. It is one of very few platforms capable of matching the scale of sequence data produced by genome sequencing. Applications range fromanalysing single base changes, SNPs, to detecting deletion or amplification of large segments of the genome, CGH. At present, its most widespread use is in the analysis of gene expression levels. When carried out globally on all the genes of an organism, this analysis exposes its molecular anatomy with unprecedented clarity. In basic research, it reveals gene activities associated with biological processes and groups genes into networks of interconnected activities. There have been practical outcomes, too. Most notably, large-scale expression analysis has revealed genes associated with disease states, such as cancer, informed the design of new methods of diagnosis, and provided molecular targets for drug development.

At face value, the method is appealingly simple. An array is no more than a set of DNA reagents for measuring the amount of sequence counterparts among them RNAs of a sample. However, the quality of the result is affected by several factors, including the quality of the array and the sample, the uniformity of hybridisation process, and the method of reading signals. Errors, inevitable at each stage, must be taken into account in the design of the experiment and in the interpretation of results. It is here that the scientist needs the help of advanced statistical tools.

Dr. Stekel is a mathematician with several years of experience in the microarray field. He has used his expertise in a company setting, developing advanced methods for probe design and for the analysis of large, complex data sets.

INTRODUCTION

The image of the microarray generated by the scanner (Section 1.3) is the raw data of your experiment. Computer algorithms, known as feature extraction software, convert the image into the numerical information that quantifies gene expression; this is the first step of data analysis. The image processing involved in feature extraction has a major impact on the quality of your data and the interpretation you can place on it.

In Chapter 1 we discussed three technologies by which microarrays are manufactured: in-situ synthesis with the Affymetrix platform, inkjet in-situ synthesised arrays (Rosetta, Agilent and Oxford Gene Technology) and pin-spotted microarrays. This chapter focusses on pin-spotted arrays. Affymetrix has integrated its image processing algorithms into the Gene Chip experimental process and there are no decisions for the end-user to make. Inkjet arrays are of much higher quality than pin-spotted arrays and do not suffer from many of the image-processing difficulties of spotted arrays; also, Agilent provides image-processing software tailor-made for their platform, so there are no decisions for the end-user either.

Pin-spotted arrays, on the other hand, provide the user with a wide range of choices of how to process the image. These choices have an impact on the data, and so this chapter describes the fundamentals of these computational methods to give a better understanding as to how they impact the data.

FEATURE EXTRACTION

The first step in the computational analysis of microarray data is to convert the digital TIFF images of hybridisation intensity generated by the scanner into numerical measures of the hybridisation intensity of each channel on each feature. This process is known as feature extraction.

INTRODUCTION

Chapter 1 introduced microarray technologies and discussed the use of microarrays in the laboratory. The remainder of the book is dedicated to microarray bioinformatics. This chapter, together with the next chapter, discusses the bioinformatics required to design a DNA microarray. In this chapter, we look at the sequence databases that are used to select and annotate the genes that the microarray detects and, thus, the sequences that will appear on the array. Chapter 3 looks at the computer design of oligonucleotide probes for oligonucleotide arrays.

There are two broad questions and one more specific consideration that this chapter seeks to address:

1. What resources could I use to design my own custom array?

If you are designing a custom microarray to study a particular disease, tissue or organism, you will need to identify the genes that might be expressed in your samples and identify the sequences of those genes. One of the aims of this chapter is to give an understanding of which databases you could use to select such genes.

How can I find more information about the sequences of the genes on my array?

DNA microarrays contain sequences that will have derived from DNA sequence databases. The output file containing the numerical results of the microarray experiment that you will analyse also contains a number of fields that relate these sequences to the databases from which they derive. This chapter describes the meanings of these fields and the nature of the databases.

INTRODUCTION

One of the most exciting areas of microarray research is the use of microarrays to find groups of genes that can be used diagnostically to determine the disease that an individual is suffering from, or prognostically to predict the success of a course of therapy or results of an experiment.

In these studies, samples are taken from several groups of individuals with known pathologies, outcomes or phenotypes and hybridised to microarrays. The aim is to find a small number of genes that can predict to which group each individual belongs. These genes can then be used in the future as part of a molecular test on further individuals, either using a focussed microarray, or a simpler method such as quantitative polymerase chain reaction (PCR).

EXAMPLE 9.1 DATA SET 9A

Bone marrow samples are taken from 27 patients suffering from acute lymphoblastic leukemia (ALL) and 11 patients suffering from acute myeloid leukaemia (AML) and hybridised to Affymetrix arrays. We want to be able to diagnose the leukemia in future patients using either Affymetrix technology or usingmore focussed arrays with a small number of genes. How do we choose a set of rules to classify these samples?

The development of such predictive models depends on statistical and computational techniques, many of which are still the subject of active research. There are essentially three parts to developing a predictive model, and so the chapter is arranged into three further sections:

1. Section 9.2: Methods of Classification, looks at a number of commonly used methods for distinguishing between groups of individuals based on a given set of measurements. There are several well-established methods for doing this, many of which have been shown to work well with microarray data.

2. […]