EMBL-Outstation - The European Bioinformatics Institute (EBI), Wellcome Trust Genome Campus, Hinxton, Cambridge CB10 1SD, UK;

Corresponding Author:
Jeremy D. Parsons
Boinformatics Consultancy
37 Walpole Road,
Cambridge,
CB1 3TH
England

Tel. +44 7986 534110

http://www.littlest.co.uk/

Abstract

Motivation:

Results:

Availability:

Contact:

Introduction

CORBA and biological data

Sequence databases (for example the EMBL NSDB Stoesser et al, 1998) are currently the most familiar biological database and a typical source for EST sequences. Traditionally these databases are presented and distributed using "flat file" views of the sequence originally submitted, but these static raw views are being superseded by novel interfaces better suited to the needs of academic and industrial researchers. The new browsing interfaces are needed because current research may need to investigate nucleotide or protein sequences grouped by their function, phylogeny, map position, cellular location, expression regulation, disease association, annotation similarity or metabolic pathway etc. To support these new investigative cross-database clients, it is useful to have standardized API level access to both data and services - a requirement where CORBA is starkly appropriate. In summary, ESTs are  inherently difficult data to handle (yet extremely useful) so they are a natural target for new efforts to improve data quality, cross database browsing and novel visualization tools.

EST databases and clustering extensions

ESTs have many problems that stem from their means of production. The cDNA library from which any ESTs are drawn will sample the levels of expression in a particular tissue at a particular time: rare transcripts will be missed and highly expressed genes will be overly abundant. This latter problem, redundancy, is both wasteful and difficult to handle due to the extreme volume of error prone data. Redundancy can be reduced in the laboratory using normalization techniques but these are imperfect though additional subtractive techniques are being investigated (Bonaldo et al, 1996).

To address the redundancy problem using computer techniques, a variety of EST clustering programs have been implemented for example Parsons et al (1992), Boguski and Schuler (1995), Parsons (1995), Eckman et al. (1998), and Yee and Conklin (1998). Once clustered, gene fragments can be assembled using sequence assembly programs, or both steps can be combined in large-scale EST assembly and gene indexing protocols. An assembly stage adds two main benefits: first, it produces contigs and consensus sequences which can completely hide EST redundancy and, second, it should also improve the length and quality of the gene reconstructions beyond that available from any one EST (usually single pass reads). For some applications, the clustering is all that is required, for example, the early stages of STS mapping but for others, sophisticated cluster partitioning (Burke et al, 1998) or complex assemblies are more suitable.

Whilst the redundancy problem can be mitigated, other deficiencies of EST databases are computationally insoluble and most EST based projects have centered on pragmatic scientific goals cognizant that perfect results are impossible. Early projects that used EST clustering include TIGR's THC database (Adams et al, 1995), Geneexpress (Houlgatte et al, 1995), and UniGene (Boguski and Schuler, 1995) but new genomes (e.g. Arabidopsis, Rounsley et al., 1996), new methods (e.g. SANBI STACK, http://www.sanbi.ac.za/Dbases.html), and new public data (Aaronson, 1996) have kept the field busy with important discoveries and enhanced views (Eckman et al, 1998).

Assembly and Gene Indexing

EST clustering on typical cDNA/EST collections will produce clusters which may vary in size from one to approximately 20,000 which is a size range commonly encountered with larger shotgun assembly projects. Some DNA shotgun assembly programs (necessary for all large-scale genomic sequencing) have been applied to EST assembly, however there is a significant danger of producing artificial sequences by incorrectly joining, or homogenizing transcripts from different genes. In summary, postulated transcription products will ultimately need experimental biological confirmation and the totality of products from a single gene may not fit well into a FASTA format file: instead, a description is needed with the potential to record all of the reconstructed biological complexity.

EST clusters can be combined with other sources of gene sequence information, (for example full-length mRNA, and genomic DNA) to better sample an organisms entire genome, (especially for rare transcripts) to create "gene indexes", and consensus sequence databases. As gene mapping information accumulates, so the gene indexes lead to the creation of "gene maps" which link genes and their functions to locations on a physical map. Gene maps supported by polymorphism information make an extremely valuable resource for positional candidate disease gene hunting.

Introduction to EST Clustering Algorithms

Computational Complexity

Experimental and Biological Complexity

Comparisons of Existing Programs

Published references for EST comparison programs show a variety of specified goals and assumptions. The simplest use scripting languages (e.g. Perl) to run and parse sequence database searches by BLASTN (Altschul et al, 1997, http://www3.ncbi.nlm.nih.gov/BLAST/) or FASTA (Pearson et al., 1997, http://www.cs.virginia.edu/brochure/profs/pearson.html). EST sequences may be clipped to various levels of quality, vector trimmed, and low complexity masked before being sent for BLASTN/FASTA comparisons. BLASTN and FASTA hash the query sequence before scanning a sequence database. JESAM and ACEmbly hash the whole database (indexing) and find hash hits for the query against all database sequences simultaneously before checking for alignments. Construction of Compugen's EST clusters (LEADS, http://www.cgen.com/solutions/index.html) uses a proprietary non-alignment based algorithm. Other (published) non-alignment algorithms include EMBLSCAN (Bishop and Thompson, 1984) D-squared Cluster,  and RAPID (Miller et al, 1999) try to estimate the significance of discovering oligos (DNA words) in common between two sequences without the cost of a full alignment.

The simple ideal that all cognate ESTs should cluster together could be achieved easily if all 3-prime ESTs were anchored to the same region of the cDNA and gene. In this situation (ignoring sequencing errors) any EST from the same gene would have the same basecalled sequence . The original ICAtools (Parsons et al., 1992) were the first tools built to automate this simple analysis. Unfortunately, as described by Gautheret et al. (1998), alternate polyadenylation is common (in humans) so ESTs may not be well anchored and discovering cognate ESTs reverts to solving a peculiar massive shotgun assembly. The anchored EST assumption will, however, still give useful results for strictly relative comparisons of the complexity within individual cDNA libraries, especially if transcript heterogeneity is temporally and anatomically restricted (Gautheret et al, Burke et al).

BLAST and FASTA based scripts

Glaxo's "Dynamic" assembler (Gill et al, 1997) uses interactive rounds of BLASTN searches of  EST databases to identify possible overlapping sequences which are fed to Gap (Bonfield et al, 1995) for Greedy assembly. The approach taken by Glaxo was quick to implement, and effective for simple assemblies but suffers pauses between each round of assembly, the Greedy assembly is prone to errors with repetitive sequence regions, the static HTML output limits interactivity/editing and alternative BLASTN alignments displayed simultaneously in one alignment may cause confusion for the user.

The original UniGene package (Boguski and Schuler, 1995, Schuler et al, 1996) relied on a BLASTN-like hash based search to identify pairs of sequences with at least two 13-mer words separated by no more than two bases. These sequence pairs were then aligned more accurately within a 10-base window either side of the diagonal identified by the oligo match. The alignment was scored +1 for a nucleotide identity, -2 for a mis-match, -1 for a gap, and 0 for an N ambiguous base. An alignment quality was calculated by dividing the score by the alignment length and had to exceed 0.91 to be accepted. An extra constraint was that the aligned region had to extend within 35 bases of the end of either sequence. UniGene itself only accepts data of known good quality, (high quality base count, presence of poly-A tail etc.) so not all ESTs that might pass the above alignment criteria are ever checked for matches. Explicit lists of excluded sequences (>200,000) are not published nor assemblies, alignments or consensus sequences. UniGene is available online at http://www.ncbi.nlm.nih.gov/UniGene/index.html with extra descriptions at  http://www.ncbi.nlm.nih.gov/Schuler/Papers/ESTtransmap/

Specialized Clustering Programs

CLEANUP (Grillo et al, 1996) is similar to the ICAtools but performs a unique hashed-query pairwise sequence comparison which makes it slower than ICAass. Query sequences are encoded as hashed oligos along with one base mutated versions of the oligos to enhance query sensitivity. The largest published data set: 2400 Drosophila sequences, was self-compared in 160 seconds. All code is written in C, and publicly available.

Small-Scale Assembly Programs

Gap4 ( http://www.mrc-lmb.cam.ac.uk/pubseq/overview.html ) implements a simple algorithm (relatively slow) where new shotgun reads are compared against contig consensus sequences looking for hash hits that extend the contigs in a Greedy fashion. The program is memory efficient and scales to megabase-sized YAC assembly projects with more than 10,000 reads. Gap4 can read and extend contigs produced by other programs e.g.. Phrap, and CAP2.

Phrap (unpublished, http://bozeman.genome.washington.edu/index.html) has incorporated error probability information deeper into assembly and editing than any other assembly package. Alignments are stored in detail, as a directed graph in RAM, leading to large memory requirements (1 GB for a YAC assembly). A consensus sequence is extracted by using a mosaic approach where the best read in any multiple alignment is directly copied to the result. Phrap has no specific features to handle EST idiosyncrasies like alternative splicing etc. Assemblies are fast (6 hours for 49,000 sequences on a single 195 MHz R10000) however a single sequence addition requires all previous calculations to be repeated. An optimized and parallelised version of Phrap is commercially available (http://www.spsoft.com/).

CAP3 (unpublished, http://genome.cs.mtu.edu/cap/cap3.html) is an enhancement of CAP2 (Huang, 1996) - both are relatively slow. CAP3 can utilize ancillary information to guide sequence assembly: mate-pairs, clone size, orientation/location, and basecall error estimates. CAP3 is used at MIPS to assemble UniGene clusters (http://www.mips.biochem.mpg.de/proj/human/)

TIGR Assembler, (Sutton et al, 1995) is a relatively fast sequence assembly program that has been proved assembling a diverse collection of repetitive microbial genomes. Memory usage for the query hashing stage is 32 * 4ⁿ bytes where n is the word size (oligo size) e.g. 512 MB for 12-mer oligos. Four criteria control the match decision: minimum length of overlap, maximum number of local errors, percentage of best possible score, and maximum length of overhang. TIGR use a combination of BLASTN and TIGR Assembler to produce their own consensus sequences. This combination is relatively slow: multiple-months to produce THC (Tentative Human Consensus) sequences for all the human ESTs. The THC database is biased towards splitting EST clusters due to this stringent assembly stage (Burke et al.,1998). TIGR Assembler on its own has been used to assemble 186,000 ESTs on a single SPARC processor, with 512 MB RAM, in seven days.

ACEmbly (http://alpha.crbm.cnrs-mop.fr/acembly/) is a fast hashed-database program where oligos are sampled from all sequences and used to prime alignments and then assemblies. ACEmbly fragment ordering is initially based on a distance metric calculated from the number of oligos two sequences have in common divided by their total number of sampled oligos. This "zipping" is fast (< 10 seconds for 600 reads) but prone to errors with repeated sequences. JESAM, the subject of this paper is another hashed-database sequence comparison program.

Word-Based Clustering

The published work on D-Squared Cluster (Hide et al, 1994, http://www.bchs.uh.edu/stemple/D2/d2.html) was preliminary and experimental with only limited exploration of the kinds of optimizations available that might make improvements to the programs specificity. D-Squared cluster uses windowed (i.e. local) n-mer oligo frequency comparisons between a query sequence and a compressed database. Using a Cray Y/MP-48 supercomputer, 884 bases could be compared against 19 million bases in less than one minute.

RAPID (Rapid analysis of Pre-Indexed Data structures) Miller et al., 1999 and http://www.bioinf.man.ac.uk/rapid/) uses oligo frequencies to alter the significance of oligo matches before scoring particular pairwise matches. No clustering is offered but the program has been tested on the analogous problem of removing redundancy and contamination from DNA sequence databases. RAPID uses Receiver Operator Characteristic (ROC) curves to tune parameters such as oligo size to get database search results with maximal specificity to the query. Rapid is memory intensive but relatively fast.

Systems and Methods

Algorithm

Sequence overlap detection

To ensure that the published alignments were both biologically useful and mathematically optimal, it was thought necessary to use a dynamic programming algorithm with a sophisticated gap penalty scheme (Altschul and Erickson, 1986). However, computation time makes it impractical to compare all sequences against each other for large EST datasets with workstation  implementations of derivatives of the Smith-Waterman algorithm. Specialized hardware was an unacceptable solution due to perceived problems of cost, availability, portability, and ease of algorithm development. The JESAM alignment algorithm therefore uses dynamic programming only for the final alignments relying on the gross overall overlap being easy to find because the goal was only to discover potentially overlapping subsequences, not distant homologues mutated apart though millennia

Alignment Algorithm Stages

For all the alignments described in this paper, an affine gap scoring scheme (+5 match, -15 mismatch, -15 gap-start, -15 per-gap-base) was used with end gaps only penalized if not-consistent with contiguity and an overall score of 200 needed for an alignment to be stored. To ensure sensitive alignments are calculated and stored, no predefined repeats are masked. Furthermore, repeat discovery is easier to do accurately after all alignments are calculated.

Supercluster creation

Implementation

Alignment calculation and storage.

The distributed architecture of JESAM provides a flexible framework to allow extension to more computationally intensive sequence comparisons (for example base call error probabilities) and could also support other costly calculations like assembly contig building via directed graph explorations. The PVM slaves running sequence comparisons are typically CPU-limited so the performance of the PVM implementation on networks of cheap uniprocessor PCs is equivalent to that of more expensive multi-processor servers.

Updates do not require recomputation of known alignments because all sequence overlaps are automatically stored by the Berkeley DB routines. If a new alignment crosses a user-determined threshold score, then it is saved complete with the position of all indels, starts and ends. The alignments are offered by a CORBA interface provided by a C++ wrapper around the same ANSI C code used to create the alignments.

A Java client recursively steps through the matches to make superclusters which are stored in a serialized Java Hashtable object which loads and saves quickly. The Supercluster server took only one afternoon to write thanks to the alignment data CORBA interface, simple algorithm, and the innate productivity of Java. Both the alignment and supercluster data can be viewed interactively using web-accessible Java applets or larger datasets can be saved to client file systems for use in other file based scripts etc. The JESAM clusters have been integrated with other databases at the European Bioinformatics Institute, see URL http://corba.ebi.ac.uk/EST for browser demonstrations.

Figure 1

Figure 2

Results and Discussion

JESAM Performance

Table 1

Total time taken to align all 123459 Rat ESTs ( in EMBL and EMBLnew on 15th September 1999, totaling 48461025 bases) against each other and then store the alignments when using different numbers of slaves. The master processes were run on a two processor (400 MHz) UltraSPARC server and the slaves were run on separate (248 MHz) UltraSPARC servers all running Solaris 2.5 or 2.6.

Table 2

Peak Core Memory Usage (KBytes) for all processes creating the alignments for Table 1. All processes were sampled at 10-minute intervals. Figures shown are for the largest memory usage by any instance of each process at any sample time so actual memory usage is less than shown.

Table 3

Peak Virtual Memory Usage (KBytes) calculated as per Table 2
 

Clustering Comparisons

The UniGene investigation detailed in Fig. 3 and Fig. 4 was small but it sampled examples of many of the known or expected properties of the UniGene databases, and some of the advantages and disadvantages of JESAM's simplistic superclustering client. UniGene Rat excluded >20,000 sequences included in the corresponding JESAM alignments and clusters, and UniGene mouse similarly excluded >150,000 sequences. The two rat databases were generally more similar than the two mouse databases (data not shown). Perhaps most significantly, JESAM both clustered more sequences and created more clusters for both data sets despite the use of genomic clones in UniGene (excluded from JESAM) that helped clustering in two of the ten mouse clusters examined (but no rat cluster). However, one rat UniGene cluster (Rn. 34368) did probably correctly combine two large divergent cDNAs which may encode alternative transcripts that JESAM rejected as being clearly non-overlapping sequences. The worst feature of the current published JESAM superclusters was demonstrated with the 12 sequences in UniGene Mm. 2477 which were combined into one supercluster of size 116977 with many other sequences presumably joined on repetitive DNA. This supercluster would not form with a more selective scoring scheme but sensitivity would then be lost. Biasing towards sensitivity minimizes information loss yet still significantly reduces the problem scale for later cluster partitioning software components. With all the alignments available on-line, it should be relatively easy to develop and test some quite sophisticated partitioning or assembly algorithms so these superclusters should not be seen as final answers in contrast to the other EST cluster databases which cannot be programmatically interrogated.

Figure 3

Figure 4

Figure 5 optional suggestions:

Rat Cluster Size Histogram

Extensions

The mouse EST and mRNA clusters (downloaded directly from the CORBA servers) have been successfully assembled using CAP3 by Juha Muilu (personal communication) as a useful extension of the JESAM superclusters. CAP3 however ignores the pairwise alignment data so an opportunity still exists for a more integrated, higher performance assembly tool tailored to EST data. A suitable post-assembly goal would be something like the Pangea CRAW tools (Burke et al., 1998 and Chou and Burke, 1999) which demonstrate how complicated patterns of expression variations can be summarized and presented in a Java applet alongside supporting analyses.

Linking the EST clusters with other databases would make them more valuable. The EMBL outstation European Bioinformatics Institute (EBI) amongst others (including the OMG's Life Science Research (LSR) Domain Task Force (DTF) (URL: http://lsr.ebi.ac.uk/) are working towards standards for object interfaces and components for many kinds of biological data. Some of the EBI's effort's include physical and genetic maps (Jungfer and Rodriguez-Tomé, 1998), annotated DNA and protein sequences, and sequence traces (Parsons et al., 1999).

JAVA RMI could have offered a CORBA alternative but was not investigated due to the lack of relevant biological standards efforts, frameworks, language independence, services and local support. AcePerl (Stein and Mieg, 1998), is another potential alternative.

Conclusion

The use of  CORBA interfaces in the JESAM package requires more software libraries (and greater expertise to install and maintain them) than would a normal file-based Input/Output system but they also offer greater flexibility. Furthermore, by maintaining all information accessible online, it is possible to give users an evidence trail that they can follow all the way from cluster, to alignment, to sequence and even back to the original sequence trace.

There are currently a few problems deploying CORBA based applications over the Internet. These problems include: firewalls blocking the IIOP protocol, the need to download ORB classes to clients due to the lack of a guaranteed local ORB, and a lack of support for multiple applet signing which would allow applets to follow object references to objects on computers other than the original applet's host. Some of the Java and CORBA problems should decrease once  JAVA 2 from Javasoft (URL:http://www.javasoft.com/) becomes a standard component of every client workstation. Amongst the many improvements over Java 1.1 is a built-in CORBA ORB allowing distributed computing to extend out of the Web browser and directly into a user's general molecular biology application set.

References

Acknowledgements