ARTICLE

TOP ABSTRACT ARTICLE REFERENCES

The rapid generation of single-read sequence from the 3' ends and 5' portions of sufficiently expressed mRNAs (popularly referred to as ESTs) (Adams et al. 1991; Okubo et al. 1991; Wilcox et al. 1991) has resulted in the discovery of many genes well before the projected completion of the human genome project and before the completion of sequencing efforts in other organisms (Adams et al. 1992; Matsubara and Okubo 1993; Venter 1993). Because the source information is available for every EST, the intracluster representation of libraries from discrete disease and developmental states can be contrasted and hence large-scale expression studies can be performed (Okubo et al. 1992, 1994; Adams et al. 1995; Vasmatzis et al. 1998). EST sequence has enabled the construction of a physical map of the human genome (Hudson et al. 1995) as well a gene map that localizes many genes with respect to the markers of the physical map (Schuler et al. 1996). The utility of EST data has also been increased greatly by the establishment of centralized databases (Boguski et al. 1993).

The fragmented nature and vast quantity of EST data pose an obstacle to harvesting the full potential from this data source. Hence, several projects are in progress to construct information frameworks, called gene indices, where the fragmented, error-prone EST data and the known gene sequence data can be consolidated and placed in a correct pathologic and mapping context indexed by gene such that all data concerning a single gene is in a single index class and each index class contains the information for only one gene. Algorithmically, these projects all comprise some type of cluster analysis in which sequence similarity and possibly other criteria are used to form the clusters or index classes. Below, we detail some of the clustering methods used in several gene index projects.

The Institute for Genome Research (TIGR) Gene Index (TGI) (http://www.tigr.org/tdb/hgi/hgi.html; Adams et al. 1995; Sutton et al. 1995; White and Kerlavage 1996) is constructed by assembling full-length sequence [from the Expressed Gene and Anatomy Database (EGAD); White and Kervalage 1996)] and ESTs to form tentative human consensi (THCs). The THC_BUILD program (G. Sutton, pers. comm.) constructs index groups according to the following schedule: (1) BLAST and FASTA (Pearson 1990) are used to identify all sequence overlaps, (2) all detected overlaps are stored in a relational database, (3) the CLUB program forms transitive closure groups from the overlap database, and (4) groups are subjected to assembly using TIGR assembler (Sutton et al. 1995; http://www.tigr.org/hgi/hgi_info.html). The assembler also imposes matching constraints on the ends of sequences and a minimum sequence identity within an index group. Most sequence assembly programs share similar properties with TIGR assembler; however, it is worth noting that some assemblers, such as the PHRAP package, incorporate sequence quality data derived from sequence traces into the assembly process (P. Green, pers. comm.) allowing for the incorporation of higher error data.

In UniGene (Boguski and Schuler 1995; Schuler at al. 1996), genes are indexed by forming initial groups with full-length sequences (mRNA and genomic). Then groups are formed within the EST set and between the ESTs and the initial groups. EST matches are not allowed to join distinct initial groups and clusters without a polyadenylation signal or unless at least two 3' ESTs are discarded. Finally, clone information was used to further join clusters and singleton clusters, and unmatched ESTs were added to nonsingleton clusters at lower stringency levels (http://ncbi.nlm.nih.gov/UniGene/TXT/build.html). To enhance the speed of the clustering, a two-phase searching process was used in which two sequences were compared with a constrained Smith-Waterman algorithm only if they shared two common words of length 13 separated by no more than 2 bases. The STACK gene indexing system that uses d2_cluster is covered in the discussion section of this paper and in a companion paper.

In this article we describe d2_cluster, an algorithm for clustering sequences into index classes. First, we describe the basic method (which we call D20 to distinguish it from future variants of the algorithm). We then demonstrate the utility of d2_cluster by performing a data analysis of the results of clustering a moderate sized data set (~43,000 sequences). We characterize the relative performance of d2_cluster and UniGene clustering showing that although the behavior of the two algorithms is very similar, there are measurable differences in the rate of merging sequences into clusters. Finally, we derive estimates of the absolute type I and type II error rates (the probability of under or over clustering) for d2_cluster.

Description of the d2_cluster Method (D20)

d2_cluster is an agglomerative clustering method [every sequence begins in its own cluster, and the final clustering is constructed through a series of mergers (Johnson and Witchern 1994)]. d2_cluster can be described in terms of minimal linkage clustering (sometimes called single linkage or transitive closure in the sequence analysis literature). The term transitive closure refers to the property that any two sequences with a given level of similarity will be in the same cluster; hence, A and B are in the same cluster even if they share no similarity but there exists a sequence C with enough similarity to both A and B. The criterion for joining clusters is the detection of two sequences that share a window of (Window_Size) bases that is (Stringency) percent or more identical. The only criterion for clustering is sequence overlap and source or annotation information is not used. To detect the overlap criterion, we use the d2 algorithm and set parameters and threshold values as described in previous work (Torney et al. 1990; Hide et al. 1994; Wu et al. 1997). The initial and final state of the algorithm is a partition of the input sequences in which each sequence is in a cluster and no sequence appears in more than one cluster.

For ease of notation, let the following conventions hold:

1.	We signify the d2 distance between two sequences, say A and B, as d2(A,B).
2.	Given two clusters, e.g., clusters i and j, the operation MERGE(cluster i, cluster j), also denoted MERGE(cluster i  cluster j), means that all sequences in cluster j are assigned to cluster i.
3.	The database to be clustered contains N sequences that are numbered 0 through (N  1). Let sequence (i) be denoted Si or S(i).
4.	The membership of sequence Si is denoted Ci.

The notation d2(A,B) is conveniently used, but, of course, d2(,) is not a function of only A and B but also of various parameters (specified in Torney et al. 1990; Hide et al. 1994; Wu et al. 1997). The MERGE operation can be expressed in terms of convention 4 above: For all sequences, Sr, such that Cr = j, Cr is reset to be Cr = i.

We describe the progression of d2_cluster inductively in that we first detail what happens in the first two iterations (I0 and I1) and then describe how one performs iteration (i + 1) given that iteration (i) has been completed. Technically speaking, it is sufficient to state only the first step and then to give the step (i) to step (i + 1) instructions, but we detail the first two steps for clarity.

Initial state: Each sequence is in its own cluster. (i.e., Si is in cluster i or Ci = i). First iteration I1: The first sequence in the database, S0, is selected as a query. For each sequence in Si (1  i < N), MERGE(cluster C0)  (cluster Ci) if d2(S0,Si) < THRESHHOLD. Second iteration I2: The second sequence in the database (S1) is now selected as a query. Note that C1 = 1 unless sequence 1 was merged into cluster 0 during step I1. For all sequences, Si (2  i < N), MERGE(cluster C1  cluster Ci) if d2(S1,Si) < THRESHHOLD. (k)th iteration I(k): Suppose we have completed (k  1) iterations. We select sequence Sk as a query. For all seqs, Si (k + 1  i < N), MERGE(cluster Ck  cluster Ci) if d2(Sk,Si) < THRESHHOLD.

The clustering is finished when N iterations are completed. Transitive closure is obtained because clusters are joined if they contain any sequences with sufficient identity. d2_cluster, as described above, can be mapped to the minimal linkage algorithm commonly seen in the statistics and engineering texts, and details of this are given in the online supplement to this paper (http://www.genome.org and www.pangeasystems.com).

Data and Error Analysis

Relative Performance of d2_cluster and UniGene

View larger version (29K): [in this window] [in a new window]   Figure 1   Subsetting comparison of UniGene and d2_cluster. Cluster equivalence means that all elements in one cluster are also present in the other cluster, and vice versa. Out of 14,989 (= 15,226  237) original UniGene clusters and 13,755 d2 clusters, 12,389 (or 83% of UniGene clusters and 90% of d2 clusters) are equivalent. Two hundred thirty-seven UniGene clusters were not considered in the analysis because they were composed of sequences that were screened out in our vector and repetitive elements screening stage.

View larger version (12K): [in this window] [in a new window]   Figure 2   (A) CRAW report (Burke et al. 1998) for a cluster formed by d2_cluster that contains two UniGene clusters: Rn.8 and Rn.3110. (B) (Available as on online supplement to this paper at www.genome.org and at the authors web site at www.pangeasystems.com) Interleaved sequence alignment shows a >300-bp region of near perfect match.

View larger version (63K): [in this window] [in a new window]   Figure 3   (A) CRAW report for a d2 cluster containing isozymes of mouse cytochrome P-450. Seven UniGene clusters (Rn.10843, Rn.3586, Rn.18603, Rn.10842, Rn.9104, Rn.11043, and Rn.15544) are merged. (B) (Online supplement available at www.genome.org and www.pangeasystems.com) Interleaved multiple alignment showing a region of 240 bases with high identity alignment between all four cluster assemblies. d2_cluster has put all of these sequences together because of regions of high identity (as seen in Fig. 3B). UniGene has separated isozymes into distinct clusters, although UniGene clusters Rn.18603, Rn.10842, and Rn.9104 should probably form a single cluster according to reasonable clustering rules due to their perfect assembly into subgroup 1 and high overlap.

Estimating the Absolute Error Rate

DISCUSSION

We have characterized d2_cluster and described the algorithm in terms that should be familiar to statisticians, computer scientists, and biologists alike. It has been shown that d2_cluster performs quite favorably to, and is consistent with, current EST clustering methods. In an empirical study based on >40,000 available rat EST and mRNA sequences, d2_cluster produced results that were between 83% and 90% identical with UniGene although d2_cluster created 10% fewer clusters and 20% fewer singleton clusters than UniGene. Three different measures of joining strength are used to compare the overall sensitivity of d2_cluster and UniGene, and these numbers are averaged to provide an estimate that d2_cluster joins sequences at a rate ~13% higher than UniGene. It remains undetermined, however, if this higher join rate results in more accurate index classes or is simply due to the joining of paralogous genes or other phenomena. Additionally, the absolute correctness of groups formed by d2_cluster has been quantified, and the sensitivity and selectivity are shown to be >99% (i.e., type I and type II error rates are bounded above by 1%).

d2_cluster has found application in the STACK project (Hide et al. 1997) in which ESTs are hierarchically clustered within tissue and arbitrary source categories. d2_cluster is set to join ESTs that are >96% identical over a window of 150 bases. More detail on STACK is given elsewhere (R.T. Miller, A.G. Christoffels, C. Gopalakrishnan, J. Burke, A.A. Ptitsyn, T.R. Broveak, and W.A. Hide, in prep.).

Unfortunately, space limitations prevent elaboration of the many other tools for sequence clustering that have been developed to cluster DNA sequence or to remove redundancies from sequence sets (Houlgatte et al. 1995; Parsons 1995; Grillo et al. 1996; Gill et al. 1997; Eckman et al. 1998; Yee and Conklin 1998; Pietu et al. 1999). Significant research has also been put into the grouping of protein sequence and the determination of domain boundaries (Sonnhammer and Kahn 1994; Worley et al. 1995; Sonnhammer et al. 1997; Gracy and Argos 1998a,b).

Sequence clustering is of little consequence in and of itself, and the prime motivation is to obtain biological knowledge. Effective sequence clustering is an organizing principle that serves as a starting point for discovery. For example, with all sequence information corresponding to a single gene in a cluster, features such as alternative exons and aberrant splicing, among others, can be modeled with greater ease. It is difficult to tune the parameters of primary sequence clustering such that features like splicing differences and even artifacts, such as chimerism, are accounted for while simultaneously generating proper index classes. Instead, decoupling this feature detection and artifact correction from the clustering step allows these problems to be handled in a more robust fashion. Hence, postprocessing steps, such as CRAW, are used to contrast gene variants and correct for artifact. Figure 4 shows how ESTs and mRNAs from a cluster of human sequences corresponding to human RAD1/REC1 cell-cycle control checkpoint protein are placed in the same cluster, whereas CRAW is used to separate distinct splice variants into separate subclusters. In a similar fashion, CRAW is also used to isolate and correct for chimeric sequence and other artifacts. Full details and additional examples are found in previous work (Burke et al. 1998; Chow and Burke 1999). Work has also been done to associate discovered multiple gene forms with sequence source information to infer state specificity or associate novel exon/UTR usage with disease (Burke et al. 1998; Gautheret et al. 1998). The d20 algorithm and others specified here are available at no cost for university researchers, and commercial licenses are available (details available from authors upon request). The commercial versions of CRAW and CRAWview (www.pangeasystems.com) were used in the preparation of this manuscript.

View larger version (73K): [in this window] [in a new window]   Figure 4   Alternative splice forms of the RAD1/REC1 gene are placed in the same cluster by d2_cluster, and the splice variants are separated into distinct subclusters by CRAW.