This appendix describes a simple yet very efficient Perl solution to a problem known as the disjoint sets problem, the dynamic equivalence relation problem, or the unionfind problem. This problem appears in applications with the following scenario.

1. Each one of a finite set of keys is assigned to exactly one of a number of classes. These classes are the “disjoint sets” or the partitions of an equivalence relation. Often, the set of keys is known in advance, but this is not necessary to use our Perl package.

2. Initially, each key is in a class by itself.

3. As the application progresses, classes are joined together to form larger classes; classes are never divided into smaller classes. (The operation of joining classes together is called union or merge.)

4. At any moment, it must be possible to determine whether two keys are in the same class or in different classes.

To solve this problem, we create a package named UnionFind. Objects of type UnionFind represent an entire collection of keys and classes. The three methods of this package are:

• $uf = UnionFind– > new(), which creates a new collection of disjoint sets, each of which has only one element;

• $uf– > inSameSet($key1,$key2), which returns true if its two arguments are elements of the same disjoint set or false if not;

• $uf– > union($key1,$key2), which combines the sets to which its two arguments belong into a single set.

This book is designed to be a concrete, digestible introduction to the area that has come to be known as “bioinformatics” or “computational molecular biology”. My own teaching in this area has been directed toward a mixture of graduate and advanced undergraduate students in computer science and graduate students from the biological sciences, including biomathematics, genetics, forestry, and entomology. Although a number of books on this subject have appeared in the recent past – and one or two are quite well written – I have found none to be especially suitable for the widely varying backgrounds of this audience.

My experience with this audience has led me to conclude that its needs can be met effectively by a book with the following features.

• To meet the needs of computer scientists, the book must teach basic aspects of the structure of DNA, RNA, and proteins, and it must also explain the salient features of the laboratory procedures that give rise to the sorts of data processed by the algorithms selected for the book.

• To meet the needs of biologists, the book must (to some degree) teach programming and include working programs rather than abstract, high-level descriptions of algorithms – yet computer scientists must not become bored with material more appropriate for a basic course in computer programming.

• Justice to the field demands that its statistical aspects be addressed, but the background of the audience demands that these aspects be addressed in a concrete and relatively elementary fashion.

Once DNA fragments have been sequenced and assembled, the results must be properly identified and labeled for storage so that their origins will not be the subject of confusion later. As the sequences are studied further, annotations of various sorts will be added. Eventually, it will be appropriate to make the sequences available to a larger audience. Generally speaking, a sequence will begin in a database available only to workers in a particular lab, then move into a database used primarily by workers on a particular organism, then finally arrive in a large publicly accessible database. By far, the most important public database of biological sequences is one maintained jointly by three organizations:

• the National Center for Biotechnology Information (NCBI), a constituent of the U.S. National Institutes of Health;

• the European Molecular Biology Laboratory (EMBL); and

• the DNA Databank of Japan (DDBJ).

At present, the three organizations distribute the same information; many other organizations maintain copies at internal sites, either for faster search or for avoidance of legal “disclosure” of potentially patent-worthy sequence data by submission for search on publicly accessible servers. Although their overall formats are not identical, these organizations have collaborated since 1987 on a common annotation format. We will focus on NCBI's GenBank database as a representative of the group.

Once a new segment of DNA is sequenced and assembled, researchers are usually most interested in knowing what proteins, if any, it encodes. We have learned that much DNA does not encode proteins: some encodes catalytic RNAs, some regulates the rate of production of proteins by varying the ease with which transcriptases or ribosomes bind to coding sequences, and much has no known function. If study of proteins is the goal, how can their sequences be extracted from the DNA? This question is the main focus of gene finding or gene prediction.

One approach is to look for open reading frames (ORFs). An open reading frame is simply a sequence of codons beginning with ATG for methionine and ending with one of the stop codons TAA, TGA, or TAG. To gain confidence that an ORF really encodes a gene, we can translate it and search for homologous proteins in a protein database. However, there are several difficulties with this method.

• It is ineffective in eukaryotic DNA, in which coding sequences for a single gene are interrupted by introns.

• It is ineffective when the coding sequence extends beyond either end of the available sequence.

• Random DNA contains many short ORFs that don't code for proteins. This is because one of every 64 random codons codes for M and three of every 64 are stop codons.

• The proteins it detects will probably not be that interesting since they will be very similar to proteins with known functions.

The human genome comprises approximately three billion (3 ×109) base pairs distributed among 23 pairs of chromosomes. The Human Genome Project commenced in the 1990s with the primary goal of determining the sequence of this DNA. The task was completed ahead of schedule a few months before this book was completed.

Scientists have been anxious to remind the public that the completion of the Human Genome Project's massive DNA sequencing effort will hardly mark the end of the Project itself; we can expect analysis, interpretation, and medical application of the human DNA sequence data to provide opportunities for human intellectual endeavor for the foreseeable future.

It is equally true, though less well understood, that this milestone in the Human Genome Project will not mark the end of massive sequencing. Homo sapiens is just one species of interest to Homo economicus, and major advances in agricultural productivity will result from ongoing and new sequencing projects for rice, corn, pine, and other crops and their pests. Because of the large differences in size, shape, and personality of different breeds, the Dog Genome Project promises to bring many insights into the relative influences of nature and nurture. Even recreational sequencing by amateur plant breeders may not lie far in the future.

Suppose we are given a strand of DNA and asked to determine whether it comes from corn (Zea mays) or from fruit flies (Drosophila melanogaster). One very simple way to attack this problem is to analyze the relative frequencies of the nucleotides in the strand. Even before the double-helix structure of DNA was determined, researchers had observed that, while the numbers of Gs and Cs in a DNA were roughly equal (and likewise for As and Ts), the relative numbers of G + C and A + T differed from species to species. This relationship is usually expressed as percent GC, and species are said to be GC-rich or GC-poor. Corn is slightly GC-poor, with 49% GC. Fruit fly is GC-rich, with 55% GC.

We examine the first ten bases of our DNA and see: GATGTCGTAT. Is this DNA from corn or fruit fly?

First of all, it should be clear that we cannot get a definitive answer to the question by observing bases, especially just a few. Corn's protein-coding sequences are distinctly GC-rich, while its noncoding DNA is GC-poor. Such variations in GC content within a single genome are sometimes exploited to find the starting point of genes in the genome (near so-called CpG islands). In the absence of additional information, the best we can hope for is to learn whether it's “more likely” that we have corn or fly DNA, and how much more likely.

The following books are especially recommended for further reading. The citations in each category are arranged roughly from most to least accessible.

Algorithms and Programming in General Bentley (1999); Garey and Johnson (1979); Cormen et al. (1990)

Computational Molecular Biology Mount (2001); Setubal and Meidanis (1997); Salzberg, Searls, and Kasif (1998b); Lander andWaterman (1995); Pevzner (2000)

General Genetics and Molecular Biology Brown (1999); Nicholl (1994); Lewin (1999); Snustad and Simmons (1999); Russell (1992)

Perl Wall et al. (1997); Vromans (2000); Srinivasan (1997)

Public-Domain Databases and Software Gibas and Jambeck (2001); Baxevanis and Ouellette (2001)

Protein Structure Goodsell (1996); Brandon and Tooze (1999)

Statistical Modeling of Sequences Durbin et al. (1998); Baldi and Brunak (1998)

In Chapter 3, we introduced the sequence assembly problem and discussed how to determine whether two fragments of a DNA sequence “overlap” when each may contain errors. In this chapter, we will consider how to use overlap information for a set of fragments to reassemble the original sequence. Of course, there is an infinite number of sequences that could give rise to a particular set of fragments. Our goal is to produce a reconstruction of the unknown source sequence that is reasonably likely to approximate the actual source well.

In Section 13.1 we will consider a model known as the shortest common superstring problem. This simple model assumes that our fragments contain no errors and that a short source sequence is more likely than a longer one. If the input set reflects the source exactly, then only exact overlaps between fragments are relevant to reconstruction. This approach's focus on exact matching allows us to use once more that extremely versatile data structure from Chapter 12, the suffix tree.

Beginning in Section 13.2, we consider another solution, a simplified version of the commonly used PHRAP program. PHRAP assumes that errors are present in the sequence files and furthermore that each base of each sequence is annotated with a quality value describing the level of confidence the sequencing machinery had in its determination of the base. It will not be possible to present all of the issues and options dealt with by PHRAP, but in Sections 13.3–13.7 we will present several important aspects of its approach in a Perl “mini-PHRAP”.

In the previous chapter we saw that finding the very best multiple alignment of a large number of sequences is difficult for two different reasons. First, although the direct method is general, it seems to require a different program text for each different number of sequences. Second, if there are K sequences of roughly equal length then the number of entries in the dynamic programming table increases as fast as the K th power of the length, and the time required to fill each entry increases as the K th power of 2.

In this chapter, we will address both of these difficulties to some degree. It is, in fact, possible to create a single program text that works for any number of strings. And, although some inputs seem to require that nearly all of the table entries be filled in, careful analysis of the inputs supplied on a particular run will often help us to avoid filling large sections of the table that have no influence on the final alignment. Part of this analysis can be performed by the heuristic methods described in the previous chapter, since better approximate alignments help our method search more quickly for the best alignment.

Pushing through the Matrix by Layers

As a comfortable context for learning two of the techniques used in our program, we will begin by modifying Chapter 3's subroutine similarity for computing the similarity of two sequences.