2710–2725 Nucleic Acids Research, 2002, Vol. 30, No. 12© 2002 Oxford University Press

Dictionary-driven prokaryotic gene finding
Tetsuo Shibuya and Isidore Rigoutsos1,*
Exploratory Technology, IBM Tokyo Research Laboratory, 1623-14 Shimotsuruma, Yamato-shi, Kanagawa
242-8502, Japan and 1Bioinformatics and Pattern Discovery Group, Computational Biology Center,
IBM Thomas J. Watson Research Center, PO Box 218, Yorktown Heights, NY 10598, USA
Received December 18, 2001; Revised and Accepted April 2, 2002

ABSTRACT
Gene identification, also known as gene finding or
gene recognition, is among the important problems
of molecular biology that have been receiving
increasing attention with the advent of large scalesequencing projects. Previous strategies for solving
this problem can be categorized into essentially two
schools of thought: one school employs sequence
composition statistics, whereas the other relies on
database similarity searches. In this paper, we
propose a new gene identification scheme thatcombines the best characteristics from each of these
two schools. In particular, our method determines
gene candidates among the ORFs that can be
identified in a given DNA strand through the use of
the Bio-Dictionary, a database of patterns that covers
essentially all of the currently available sample of thenatural protein sequence space. Our approach relies
entirely on the use of redundant patterns as theagents on which the presence or absence of genes ispredicated and does not employ any additional
evidence, e.g. ribosome-binding site signals. The
Bio-Dictionary Gene Finder (BDGF), the algorithm’s
implementation, is a single computational engine
able to handle the gene identification task acrossdistinct archaeal and bacterial genomes. The engine
exhibits performance that is characterized bysimultaneous very high values of sensitivity andspecificity, and a high percentage of correctlypredicted start sites. Using a collection of patterns
derived from an old (June 2000) release of the
Swiss-Prot/TrEMBL database that contained 451 602
proteins and fragments, we demonstrate ourmethod’s generality and capabilities through an
extensive analysis of 17 complete archaeal andbacterial genomes. Examples of previously
unreported genes are also shown and discussed indetail.
INTRODUCTION
As a testimony to the accelerated pace of genome sequencing
projects, almost 80 complete genomes have been deposited in

the public databases to date, whereas many more genomes are
currently at various stages of sequencing. Consequently, theautomated identification of the protein coding regions in a
newly sequenced genome is attracting increasing attention.
Accurate gene prediction is of relevance to many biologicalapplications. For example, the predicted coding regions can beused to generate probes for a DNA microarray, they can form
the basis for knockout experiments, the candidate proteins
corresponding to these predicted genes might be used as new
drug targets, etc. In this paper, we focus on the prokaryotic
gene identification problem. Gene identification is also known
as ‘gene discovery’, ‘gene recognition’ or ‘gene finding’—thelatter is the term we use in this discussion.
With the exception of a handful of reported instances inarchaeal organisms, splicing does not generally occur inprokaryotes. Thus, the problem of gene identification in theseorganisms is generally considered to be simpler than its
eukaryotic counterpart. The schemes which have been
proposed over the years have permitted great advances in thein silico prediction of genes in prokaryotic genomes but,arguably, have shortcomings. As such, the demand for increasingly accurate prediction schemes continues.
Over the years, a large number of methods have been
proposed that address the gene finding problem. Thesemethods can be largely divided into two categories. The
methods in the first category make use of the statistics of DNA
sequences to determine the location of genes. In fact, it wasobserved very early that the statistical properties of nucleotide
usage differ inside DNA regions which code for genes andoutside those regions: the concept of the CpG island (1) is ademonstration of such a difference in statistical behavior.
Among the gene identification methods that make use of thisobservation, those which are based on Markov models are themost popular to date (2–5).
The second category comprises methods that are based on
similarity search in large databases of genomic information(6–11). These methods carry out database searches in an effort
to determine DNA regions (respectively amino acidsequences) that share similarities with the DNA regions(respectively amino acid translations) of ORFs from thegenome under consideration. For details on such techniques,
the reader is referred to one of the many review papers on the
topic: several of these papers also address the gene identification
problem in eukaryotic organisms (12–17).
Despite the notable success of the methods that have beendeveloped over the years, each of these two basic strategies has
its own shortcomings. Statistical methods such as those based*To whom correspondence should be addressed. Tel: +1 914 945 1384; Fax: +1 914 945 4104; Email: rigoutso@us.ibm.com

Nucleic Acids Research, 2002, Vol. 30, No. 12 2711

on Markov models can identify coding regions whose statisticalbehavior is similar to that of the used training set. If no appropriate
training set is available, one resorts to using sets derived fromdatabase searches, or simply assumes that very long ORFs do
code for genes. The statistics of coding regions often differ
from organism to organism and, ideally, if one wishes toachieve high prediction ratios using such approaches then one
ought to employ models with organism-dependent parameters.
Essentially, a different Markov model must be built for eachtargeted genome. Moreover, short genes (e.g. <60–80 aminoacids) cannot be predicted reliably using statistical methods.
Finally, genes that are statistically distinct from other genes ofthe same organism, e.g. genes that are the result of horizontaltransfer (18,19), typically represent challenges for methods
based on statistical schemes.
Unlike statistical methods, similarity-based approaches are
more effective in finding short genes or genes that are statistically distinct from the majority of the genes in the organism
being studied. The implicit assumption here is that similar
genes or similar proteins are already present in the databasesthat are searched. Clearly, there is no dependence on training
sets since no such sets are needed. Problems can arise if theshared similarity between a candidate gene and its databasecounterpart is very low and not detectable. In general,similarity-based methods have an improved ability to determine the correct location of genes over statistical methods, avery desirable property for gene finding tools.
Given that the best characteristics of these two categories
complement one another, genome sequencing projects typically use representative methods from both categories togenerate results (20).
In what follows, we present a new method which borrows
the best attributes from similarity searches, while at the sametime relying on implicit sequence statistics as in the case of
Markov models. Our method builds on a new paradigm that
describes amino acid sequences with the help of patterns thatare present in these sequences. The patterns are derived byprocessing very large public databases of amino acidsequences with the help of an unsupervised discovery algorithm. It is worth noting that our method does not make any useof additional evidence, e.g. ribosome-binding sites, in order todecide the presence, absence and location of genes. Clearly,
incorporating such information provides additional constraints
that can further increase the quality and accuracy of the results
generated by a gene finding algorithm. Nonetheless, in thisfirst installment of our approach, we have decided against
incorporating such modules. We made this decision because ofour desire to present a clear assessment of our method’s potential as a generic, alternative scheme to gene finding. We plan toincorporate such additional modules in follow-up work.
In Methods and Algorithms we describe our approach indetail. The Experimental Details and Results contains details
on the experimental setup as well as a presentation and analysisof the results obtained from applying Bio-Dictionary Gene Finder(BDGF), our algorithm’s implementation, to 17 archaeal andbacterial genomes. The presentation concludes with a discussion.
METHODS AND ALGORITHMS
In this section, we present methodological details and giveinformation on the algorithms that we have employed.

Notation and definitions
Let Σ denote the alphabet of all 20 amino acids. When
processing an input dataset containing a collection of strings
from Σ+ with the Teiresias algorithm (21,22), we cansuccinctly capture the patterns that can be discovered with theregular expression Λ(ΛU{‘.’})*Λ where Λ = (ΣU[ΣΣ*Σ]), and‘.’ is a ‘don’t care’ character which stands for any character inΣ. In other words, the generated patterns can either be a single
alphabet symbol or strings that begin and end with a symbol ora bracket with two or more characters, and contain an arbitrary
combination of zero or more residues, brackets with at least
two alphabet characters, and don’t care characters. A bracket
denotes a ‘one of’ choice, i.e. [CPM] denotes exactly one of C,
P or M. Also, a bracket can have a minimum of two alphabet
characters but obviously not more than |Σ| –1.
A pattern t is called an <L,W> pattern with (L ≤ W) if every
substring of t of length W which begins and ends with a literalcomprises L or more positions that are occupied by literals. The
smallest length of an <L,W> pattern is obviously equal to L
whereas its maximum length is unbounded. Any given choicefor the parameters L and W has direct bearing on the degree ofremaining similarity among the instances of the sequence fragments that the pattern captures: the smaller the value of theratio L/W, the lower the degree of local similarity. Also associated with each pattern t is its support which is equal to thenumber of t’s instances in the processed input database.
Finally, K denotes the minimum required support and represents the minimum number of instances that a pattern t musthave before it can be reported.
The Bio-Dictionary
In previous work, we introduced and discussed in detail theconcept of the Bio-Dictionary (23–25). The Bio-Dictionary is
a collection of patterns that we refer to as seqlets (for ‘small
sequences’) and which completely describes and accounts forthe sequence space of natural proteins at the amino acid level.
The seqlets are derived by processing a large public databaseof proteins and fragments using the Teiresias algorithm (21,22)
and for an appropriate choice of L, W and K. In Rigoutsos et al.
(23), we described our computation of the Bio-Dictionary fromthe GenPept release from February 10, 1999 using L = 6,
W = 15 and K = 2. The processed input contained ∼387 000
sequences amounting to a grand total of ∼120 000 000 aminoacids. The computation resulted in a Bio-Dictionary that at thetime comprised ~26 000 000 seqlets and accounted for (i.e.
covered) 98.12% of the amino acid positions in the processed
input. The reader is referred to that publication for details
regarding the computation, example seqlets with discussion,and an extensive description of possible applications. Indeed,
the availability of such a complete collection of seqlets permits
one to effectively and successfully tackle a number of tasks
that include similarity searching (26), functional annotation(25), phylogenetic domain analysis (24), gene identificationand others. Below, we describe gene identification in detail.
The key idea behind dictionary-driven gene finding
As mentioned already, the Bio-Dictionary concept seeks tosubstitute a given database of proteins and fragments such as
GenPept or SwissProt/TrEMBL (27) by an equivalent collectionof regular expressions (= seqlets) each of which represents

2712 Nucleic Acids Research, 2002, Vol. 30, No. 12

Figure 1. The basic idea behind our dictionary-driven gene finding method.

combinations of amino acids that appear two or more times in
the processed input. To the extent that the input sequences insuch a public database correspond to a representative samplingof the sequence space of natural proteins, the seqlets of the
Bio-Dictionary represent an exhaustive collection of intra- andinter-family signals that are discovered in an unsupervised andexhaustive manner. Our computation of intra- and inter-family
signals becomes possible due to the fact that during processing
we consider all publicly available sequences without any of thefiltering that one encounters in databases such as PROSITE(28), Pfam (29) or PRINTS (30) and which is based on thesequences’ known functional behavior. The properties of theTeiresias algorithm guarantee that all <L,W> patterns will be
discovered and that they have the maximum possible extent
and specificity.
Two requirements need to be fulfilled for the Bio-Dictionary
approach to be successful in tackling the kinds of problems inwhich we are interested. First and foremost, the input to beprocessed should be a large and diverse collection of proteinsand fragments [see Rigoutsos et al. (24) for details], a condition
that can be satisfied given the large number of completed andongoing genome sequencing projects which contribute to the
public databases. Second, the pattern discovery process shouldbe able to generate patterns that are specific enough, notaccidental, and account for as much of the processed input aspossible. As we have demonstrated previously (23) this isindeed possible, thus the second requirement is satisfied as well.
Given this description, our strategy for gene finding should
now be evident. Figure 1 graphically depicts the basic idea.
First, we compute all possible ORFs in each of the three
reading frames and for both the forward and reverse strands ofthe given DNA sequence. Clearly, the number of true coding
regions will be a proper subset of this collection of ORFs.

Then, for each ORF we generate its amino acid translation: if
the ORF under consideration is indeed a coding one, then weshould be able to locate instances of many of the Bio-Dictionary’s
seqlets across the span of the ORF’s translation, and vice versa.
If the number of seqlets that we can locate exceeds a predetermined threshold, we report the ORF as a putative gene. We
discuss the details of thresholding below; as a rule of thumb,the higher the number of Bio-Dictionary seqlets that can befound in a given ORF, the more likely it is that the ORF iscoding for a gene.
Incorporating a weighting scheme
The basic strategy for gene identification which we justdescribed is straightforward. It is easy to see that in addition to
the number of seqlets that can be located within the translationof an ORF, the very composition of these seqlets can have an
impact when deciding whether the ORF codes for a gene. In
general, any two Bio-Dictionary seqlets that match an aminoacid translation affect this final decision differently. One canthink of each seqlet being associated with a specific score: bysumming up the scores of the individual seqlets that match an
ORF, we can compute a quality measure that will allow us todetermine whether to report the ORF as a putative gene or todiscard it. Next, we examine how to appropriately weigh eachof the seqlets.
Let T = {t1,t2,...,tn} be the complete collection of seqlets ina given Bio-Dictionary. Let us consider the amino acidtranslation s of an ORF from a given DNA sequence and let lbe the length of s. We say that a seqlet matches at position j of
the amino acid sequence s if an instance of the seqlet can befound beginning at the jth location of s. For example,
G..G.GK[ST]TL matches the sequence MTHVLIKGAGGSGKSTLAFW beginning at position 8 of the sequence. Let

Nucleic Acids Research, 2002, Vol. 30, No. 12 2713

Tsj denote the set of those seqlets that match beginning at position
j of s, and let T′sj be the set T\ Tsj . Also let Ts = { t v1, t v2, ...,t vm}
denote the concatenated list of Tsj values for all j(1 ≤ j ≤ l).
Similarly, let T ′s = { t u1, t u2, ...,t u m′ } denote the concatenatedlist of T′sj values for all j(1 ≤ j ≤ l). Note that Tsj values fordifferent js can contain the same seqlet, thus Ts is in general amultiset. Similarly, T ′s can also be a multiset.
Let pi be the probability that seqlet ti matches a databaseprotein or fragment at a fixed location, and qi be the probability
that ti matches the amino acid translation of a non-coding ORF
at a fixed location. If all the seqlets in the Bio-Dictionary are
assumed to be statistically independent then the probability
that a given ORF corresponds to an actual gene will be equal torPs /[rPs + (1 – r)Qs] where Ps = Πi m= 1 pvi × Πim=′ 1 ( 1 – pui ) ,Qs = Πi m= 1 qvi × Πim=′ 1 ( 1 – qui ), and r is the ratio of ORFs within
the given ORF set that correspond to actual genes.
Let us examine this independence assumption a little further
by considering two specific seqlets that contain don’t carecharacters and have several instances in the database fromwhich they were derived. It is easy to see that because of the
don’t care characters the seqlets can be overlapping andmatching the sequence under consideration (the ‘query’) eventhough they are derived from two distinct groups of unrelatedsequences—in such a case, the seqlets cannot really be considered dependent. Of course, it can also happen that the seqlets
are overlapping and matching the query and they are also
derived from related groups of sequences in the processed
database—in this case, the two seqlets would be dependent.
Frequently though, two seqlets would be dependent only as aresult of a small subset of sequences that is shared by the twogroups of input sequences that gave rise to the two seqlets inthe first place. It is conceivable that one could keep track ofsuch situations for all groups of seqlets whose instances
overlap within some sequences of the database. Something likethis would of course require exceptionally demanding bookkeeping (in terms of required space and time) and would makethe application of our approach prohibitive for problems ofpractical size. It should be clear that the combinations ofseqlets that ought to be considered dependent due to overlapsof some of their database instances are substantially fewer thanthe combinations of seqlets that are independent. Moreover, agiven amino acid position of the query is typically ‘covered’ by
5–10 seqlets, a small number compared with the total numberof seqlets that will match somewhere within it. Thus, any
seqlet dependence that manifests itself in the outlined gene
finding process will involve only a small number of seqlets
each time and only because of a small fraction of the total
number of their instances in the database from which they wereoriginally derived. Consequently, the assumption of independence is a reasonable simplifying approximation that also
allows for speedy computations and, as evidenced by the
results we present in the next section, does not have any noticeable adverse impact on our results.
We can use the ratio Rs = Ps /Qs to determine the likelihoodthat a candidate ORF corresponds to a gene. Note that theprobability will be >0.5 if Rs > 1 and r = 0.5. Let N = Π in= 1 (1– pi) and N′ = Π in= 1 (1 – qi). N denotes the probability that noseqlet matches at a fixed position of an actual protein, and N′
the probability that no seqlet matches at a fixed position of theamino acid translation derived from a non-coding ORF. Recall
that the Bio-Dictionary is derived from a database of proteinsand fragments, thus N is smaller than N′ in most cases. Considering that the number of seqlets that match at a given position
is much smaller than the number of seqlets not matching at thesame position, and that the p is and qis are very small
numbers, we can approximate the terms Π im=′ 1 ( 1 – p u i ) andΠ im′= 1 ( 1 – q ui ) by Nl and N′l respectively. Thus we can use Rs′= (Ps′/Qs′)(N/N′ )l instead of Rs, where Ps′ = Π im= 1 p v i and Qs′= Π im= 1 q vi . Note that if N ∼ N′, we can use R″s = P′s/Q′s insteadof R′s.
By definition, ORFs do not contain any stop codons internally, and consequently long ORF-like stretches are not likely to
be random. Let L denote the probability that a stop codon isobserved at a fixed position of a random DNA sequence. Since
there exist three stop codons among the 64 possible codons andassuming that all four possible nucleotides appear with equal
probability in the random sequence, then L is equal to 3/64.
With this in mind, we can use Rs′′′ = R′/(1 – L)l = (Ps′/Qs′)Ml
where M = N/N′(1 – L) instead of Rs′. If M ∼ 1, we can use
Rs′′ = Ps′/Qs′ instead of Rs″′.
Let wi = log pi – log qi be the weight associated with seqlet ti.
Let us also consider the sum of weights of the seqlets matching
anywhere in the translation s of an ORF as the measure Ws thatis characteristic of the coding quality of the ORF under consideration. It is easy to see that we can write the following
equation for the coding quality measure of an ORF:
W s = Σ mi = 1 w vi = Σ mi = 1 ( log p vi – log qvi ) = log Rs′′. If M cannot beignored, we can instead use the following expression: W′s = Ws
+ l × log M = log Rs′′′ to define the coding quality of an ORF.
In actuality, the value of l × log M is far smaller than the valueof Ws, and we can safely ignore the term during the actual
computations.
At times, we are faced with a situation where we havemultiple start codons matching the same stop codon and must
decide which start/stop pair to report. Our solution amounts topicking the start codon that will result in the highest value forthe coding quality measure. Due to the fact that seqlets can also
have negative associated weights, and even if we ignore the log
M term, it should be evident that selecting the start codon insuch a way will not necessarily result in the reporting of thelongest ORF as coding.
On a related note, ATG is the most frequently used startcodon but it is not the only one. Consequently, it is inappropriate to treat the different start codons in a uniform manner.
Let {c1,c2,...,ck} denote the set of possible start codons. Let fi be
the probability that ci is the start codon of a randomly chosen
coding region, f ′i be the probability that ci is observed in noncoding regions, and gi be log fi – log f ′i. We can then use Ws +gi instead of Ws as the measure of coding quality for the aminoacid translation s of an ORF that is initiated by the start codon ci.
In order to compute the coding quality measure, we need thevalues for pis and qis. The most natural way to obtain thesevalues is to compute them with the help of actual genes andnon-coding ORFs. We can calculate the actual seqlet occurrences in the regions annotated as coding in a given training set
and derive the needed pi values; we can then compute eachseqlet’s occurrences in ORFs that are not designated as coding
in a training set and derive the qi values. The values for the fis
and f ′i s can be computed in a similar manner.
But how can these values be obtained in the absence of a
training set? For the pis we can use the probabilities computedwith the help of the protein database from which the2714 Nucleic Acids Research, 2002, Vol. 30, No. 12

Bio-Dictionary is derived. Alternatively, we can compute themusing very long ORFs instead of actual coding regions. For theqis we can use non-ORF regions, or we can estimate theprobability of random occurrence based on an appropriately
chosen amino acid bias.
Once we have attached a corresponding coding-quality
measure to each ORF we can decide which ORFs correspondto putative genes by appropriately setting a threshold value.
The higher the value of the measure we associate with a given
ORF the more likely it is that the ORF is a coding one.
Removing encapsulated genes coded in different frames
Sometimes one encounters ORFs whose span completely
encapsulates other ORFs in one or more of the remaining fivereading frames. Occasionally, our method will give comparable, high scores to a pair of ORFs where one of the ORFs
completely includes the other. It is believed that not both
members of such pairs of coding regions can correspond toactual genes. In these situations, we use the ORF score to subselect, and report the one with the higher score. Note that mostof the time, the ORFs selected and reported in this mannercorrespond to the longer member of the pair. A similar
approach is also employed by Glimmer (2).
Matching seqlets to a sequence of amino acids
The Bio-Dictionary we use in our experiments contains∼30 000 000 seqlets (see the following section for moredetails). All of these seqlets need to be checked against a large
number of amino acid translations from all ORFs on bothstrands of a given complete genome. It is thus important that
this operation be carried out as efficiently as possible so as toreduce the overall computational requirements of our method.
There exist many efficient linear-time algorithms for searchingexact strings (31), but since seqlets include don’t care characters these algorithms are not applicable here. In this section, we
present and discuss two algorithms for quickly determining
whether a seqlet matches a given amino acid sequence. We canaddress this problem in one of two generic ways: we can either
preprocess the amino acid sequence, or we can preprocess the
Bio-Dictionary seqlets with which we will be searching.
Most of the work that has appeared in the literature revolved
around the preprocessing of the sequence which is assumed tobe fixed (31). Moreover, the number of patterns whose
instances were sought in the sequence was substantiallysmaller than the collection of patterns we are interested in.
Note that in our work, the set of seqlets is fixed, known inadvance and very large. On the other hand, the amino acidsequences that we need to examine are numerous and not
known in advance. Moreover, the average length of such aquery sequence is just several hundred amino acids. It thus
appears more effective to preprocess the contents of the
Bio-Dictionary. One approach would be to build hash indicesout of seqlets so as to reduce the number of seqlets that need to
be examined at each position of the amino acid sequence;
however, it is not immediately clear how to build such indicesfor patterns that have variable numbers of don’t care characters
as is the case of the seqlets.
Similarly to what is done in traditional dictionary searches
without don’t care characters, we can build the index with thehelp of the first few characters of a seqlet. For example, theseqlet G..G.GK[ST]TL could be indexed through G..G if we

consider indices that are built using the first four characters of
the seqlet. If one of the positions used to build the index is
occupied by a bracket expression like [ST] it is much simpler
to replace the contents of that position by a don’t carecharacter. This method is rather effective and can in fact be
improved further. In the case of G..G.GK[ST]TL, it would
seem more appropriate to use one of GK.T or K.TL as theseqlet’s four-character-induced index instead of the first characters G..G. Since GK.T and K.TL contain fewer don’t carecharacters, they should generally appear fewer times than G..Gin an arbitrary protein.
Keeping these observations in mind, we next propose a novel
effective scheme for generating indices out of seqlets. Let us
define an (l,w) subpattern of a seqlet t as follows: first, we
replace all brackets in the seqlet by don’t care characters andlet t′ denote the modified seqlet. The (l,w) subpattern is aminimal substring of t′ such that its length is less than w and it
contains l characters that are not don’t care characters. Noticethat we cannot always find such (l,w) subpatterns; if one suchsubpattern exists, we select it and form the index for t. If sucha subpattern does not exist we find the largest value l′ such thatan (l′,w) subpattern exists in t′ and use this subpattern to form
the seqlet’s index: for example, GK.T is a (3,4) subpattern ofG..G.GK[ST]TL.
If a seqlet t is indexed by an (l,w) subpattern, we have tocheck this seqlet approximately n/|Σ|l times when we consider
a random protein of length n that has uniform amino acid bias.
We can of course reduce searching time by setting l to a large
value. But, at each position, we must examine patterns that areindexed by any of C wl′ –– 11 possible (l′,w) subpatterns for all l′ ≤ l.
The total number of patterns to be checked at any one positionis thus O(2w); if we set w and l to large values, the search will
be slow. In the section entitled ‘Performance of the seqlet
search scheme’ we will describe how to choose appropriate
values for l and w. Note that if there exist seqlets without any fixedcharacter (e.g. [LK]..[ST][GKT]..[AERST][ELK]..[AVL]) wemust check them at each position in addition to the seqlets
hashed by subpatterns as above.
If we know the frequency of appearance of individual aminoacids in query sequences, we can use it to estimate the probability of appearance for each subpattern. In our case, we canuse it to further improve the performance. For example, in thecase of G..G.GK[ST]TL, if we know that the amino acid G
appears less frequently than the amino acid L, we can assume
that GK.T is also rarer than K.TL, and thus we should hash theseqlet with the help of GK.T. There can be cases where some(l′,w) subpattern is estimated to be less frequent than an (l,w)
subpattern although l′ < l. In such cases, we should use theformer, more rare subpattern. We can easily search for such a
subpattern in time that is linear to the size of the seqlet. In theexperiments below, and for given values of l and w, we use as
a hash index the least frequent (l′,w) subpattern, with l′ ≤ l, as
can be estimated from a given known amino acid bias.
EXPERIMENTAL DETAILS AND RESULTS
In this section, we describe and report on the results weobtained with BDGF, the implementation of our gene-finding
algorithm. BDGF was applied to several complete archaeal andbacterial genomes. We begin by describing the building of the
Bio-Dictionary that we use and the parameter choices for ourNucleic Acids Research, 2002, Vol. 30, No. 12 2715

search scheme that determines possible matches of a given
seqlet in the amino acid translation of an ORF.

Table 1. Statistics for the seqlet-derived pattern collections used in our
experiments

Generation of the Bio-Dictionary
With the help of the Teiresias algorithm (21,22), we computed
an instance of the Bio-Dictionary for the June 12, 2000 releaseof the SwissProt/TrEMBL (27) database. The processing was
carried in the manner that is outlined in Rigoutsos et al. (23).
As a matter of fact, we used Teiresias with a setting of L = 6,
W = 15 and K = 2. The justification for this choice of values for
L and W is the result of earlier extensive analysis and was
described previously (26). The Bio-Dictionary that resulted
from this processing contained 29 397 880 seqlets. The
instances of these seqlets accounted for 98.10% of the aminoacid positions in the processed database. It is this collectionthat we used in our experiments.
It should be noted that the SwissProt/TrEMBL release that
we used to build our Bio-Dictionary is >1.5 years old andchronologically preceded the releases of some of the test
genomes that we processed for our experiments. This was anintentional choice on our part and was meant to demonstrate
our method’s generality and ability to extrapolate.
Some of the seqlets in the employed Bio-Dictionary haverather long spans. These are typically important patterns,
generated by putative proteins that are coded by genes indistinct genomes from the same phylogenetic domain. Because
these patterns appear infrequently, it is difficult to compute
their weights in the absence of an extraordinarily large training
set. For our experiments, we handle this situation as follows:
we replace seqlets that have k fixed characters by k – s + 1
seqlets each of which contained s fixed characters, for somechoice of the value s (see below). Let S{i..j} denote a
subpattern of a seqlet S that starts at the ith fixed character andends at the jth fixed character. We replace S by subpatterns
S{1..s}, S{2..s + 1},...,S{k – s + 1..k}. Any duplicate seqlets
that appear in the resulting collection are removed before
further processing. Heretofore, and for simplicity purposes, wewill use the shorthand notation BD-i to refer to the derived
pattern collections that are constructed as above by setting s toi. The collections that we used in our experiments were BD-4 andBD-6.
Table 1 shows statistics for collections BD-i with i assuming
values between 2 and 6 inclusive. In particular, for each BD-ithe table lists the following items: (i) the number of seqletderived patterns contained in BD-i—as expected, and because
identical patterns are removed after splitting, the number issmall for small values of i; and (ii) the average numbers of
instances of a derived pattern per 1 000 000 amino acids,
computed from experiments against coding sequences (column 3)
and non-coding sequences (column 4) from the 17 genomes
we used for our experiments—as expected, the number of
instances is much higher in coding sequences than in noncoding ones.
The total lengths for the coding and non-coding sequence
sets for the 17 genomes that we used in our experiments are∼30 and ∼94 Mb respectively (note that the seemingly large
sizes of theses datasets are due to the existence of six reading
frames). If we make use of all 17 organisms for training, thecorresponding coding and non-coding sequence sets havesufficiently large sizes to permit the computation of representative

See text for additional discussion.

weights, even for the patterns in the BD-6 collection. On theother hand, if only one or a handful of genomes are available,
the corresponding sizes for these two sets are not substantial topermit the generation of representative weights for the BD-5
and BD-6 collections. Thus, a given choice for the training set
indirectly dictates which collections BD-i can be used to carryout any planned gene finding experiments.
Performance of the seqlet search scheme
As we have already mentioned, it is important that we be able
to quickly determine which of the seqlets of the Bio-Dictionary
are contained in the amino acid translation of a given ORF. To
this end, we present experimental results on the performance ofthe search scheme that we described in the previous section,and for different combinations of the parameters l and w.
For the purposes of benchmarking the method’s performance, we used the (i) original unmodified Bio-Dictionary that
contained 29 397 880 patterns and (ii) 100 proteins from
Escherichia coli as queries. The average span of the seqlets inthis Bio-Dictionary is 13.15 positions whereas the average
number of fixed characters in these seqlets is 7.41. The average
length of the 100 query proteins is 340.61 with the shortest andlongest sequences having lengths of 21 and 1073 amino acids
respectively. The goal of the experiment is to find all the
Bio-Dictionary seqlets that have instances in the 100 query
proteins. As it turns out, a total of 302 349 Bio-Dictionary
seqlets appear in the cumulative collection of 100 proteins,with each amino acid position participating in the instances of
8.88 distinct seqlets, on average.
Table 2 shows the computation time required to determine
all these matches and for the parameters l and w assuming
values in the range 2 ≤ l ≤ 10, 2 ≤w ≤ 15 (l ≤ w). All of theexperiments are carried out on a single IBM RS64III processor
with a clock speed of 450 MHz. For these experiments, weconstructed a hash index assuming an amino acid bias that
results from a uniformly distributed random sequence of
nucleotides. The best performance is obtained when l = 5 andw = 11: with these settings, we can process a 500 amino acidquery and determine the subset of the Bio-Dictionary’s∼30 000 000 seqlets that have instances in the query in ∼1 s.
Table S1 (Supplementary Material) shows the average
number of seqlets that are examined for instances in the query
sequence. Note that only a subset of the Bio-Dictionary’s
seqlets will actually match somewhere in the query. As anticipated, the number of seqlets that need to be examineddecreases as l and w grow larger, an expected result. Also, thesearch time does not decrease when l >5 or w >11. The

2716 Nucleic Acids Research, 2002, Vol. 30, No. 12

Table 2. Sample search times (in s) using (l,w)-subpattern-based indexing and for various choices of l and w

See text for additional details.

Table 3. Details of the 17 genomes used in our experiments

explanation can be found in Table S2 (see Supplementary
Material) where we show the maximum number of subpatterns
that need to be checked at each position of a given protein: thisnumber increases as l and w increase. As a matter of fact, theactual search times are roughly proportional to 1.6 × x + y,where x is the number of the checked subpatterns at each position and y is the average number of the actually checked seqlets
at each position. Note that x is sometimes larger than y because
subpatterns exist that must be checked but for which no seqlet
has been hashed. We can easily compute the expected numberof checked seqlets against random sequences with a given
amino acid bias, and thus can estimate appropriate values for l
and w in this manner.
For completeness purposes, we also carried out experiments
where we searched using hash keys derived from the prefixes
of seqlets. The results are shown in Table S3 (see Supplementary
Material), for various values of the prefix length. This scheme
does not perform as well as our (l,w)-based hashing schemeand the best result of 237 s that is obtained for a prefix lengthof 11 is far too slow to be useful; thus, we abandoned thisprefix scheme idea.
Gene finding results on archaeal and bacterial genomes
In this section, we outline and discuss the capabilities of ourgene-finding method by reporting the results we obtain fromprocessing 17 complete genomes with BDGF.
Genome identities. Of the 17 genomes we used in our
experiments, four were archaeal (Archaeoglobus fulgidus,
Methanococcus jannaschii, Methanobacterium thermoautotrophicum and Pyrococcus abyssi) whereas the remaining 13
were bacterial. Table 3 shows relevant information for thesegenomes including the genome length in nucleotides, thenumber of all identifiable ORFs that are >18 nt (i.e. six
amino acids), and the number of annotated coding regionsthat have been reported in the public databases for eachgenome. We should mention at this point that each of thesegenomes may contain additional actual coding regions that todate have remained unreported. Also, one should not lose sightNucleic Acids Research, 2002, Vol. 30, No. 12 2717

of the fact that the annotated (= reported) coding regions arefor the most part putative and have typically been reported
without verification via wet laboratory experiments.

Clearly, one can generate very appealing, large values for thesensitivity of the algorithm simply by lowering the employeddecision thresholds. But this is typically done at the expense of
introducing false positives in the output which will in turn leadto decreased values for specificity. The opposite situation is
also possible: one can choose thresholds in a way that will
result in high specificity values at the expense of sensitivity;
i.e. many actual genes will not be reported. Sensitivity andspecificity are competing goals and, ideally, any proposed
algorithm must aim at achieving simultaneous high values forboth of these measures.. In addition to an algorithm’s specificity and sensitivity, also of interest is the cardinality of thecollection of genes that have been predicted by the algorithm
and satisfy the following two conditions: (i) the predictedgenes are not among the genes that have been reported in the
public databases; and (ii) the predicted genes have substantial
similarity to one or more protein/cDNA sequences containedin the public repositories. Naturally, this collection forms asubset of the results that would otherwise be characterized as‘false positives’.
The existence of genes in several distinct genomes that also
exhibit similarity to a gene predicted by a given algorithm adds
support to the hypothesis that this gene is indeed correctlypredicted. In what follows we use the term ‘hits’ to refer to themembers of the special subset of predicted genes that alsosatisfy conditions (i) and (ii) above. In our experiments, we
determined whether a predicted gene satisfied condition (ii) byusing both the FASTA (32) and the BLAST (33) algorithms:
with default threshold and matrix settings we carried our similarity searches against the release of SwissProt/TrEMBL (27)from September 21, 2001. A query was considered to generate
a hit in the searched database if one or more of the reported
results had associated E(.) values that were ≤1.0e–4 for
FASTA and 1.0e–3 for BLASTP (these specific thresholdvalues were recommended by one of the reviewers). In addition to running FASTA and BLAST, we carried out a CD
search for conserved domains using rpsblast and theConserved Domain Database from February 28, 2002 (34): theE(.) value threshold we used here was equal to 1.0e–4.
In all cases that are described below, we quantified theperformance of our approach by simultaneously reporting thevalues of the following three measures: ‘sensitivity,’ ‘specificity’and ‘hits’.

order to compute weights for all the seqlets in the BioDictionary. The experiments that we carried out and whose
results are presented in detail in the section entitled ‘Prediction
results on the various genomes’ were meant to mimic the very
wide spectrum of situations that a researcher may encounter ina real-world setting.
First, we divided each genome into two equal-length parts:
we used the second half of the genome as a training set and thefirst half as a test set (case 1). In spirit, this is a test similar towhat has been previously reported in the literature (2,5).
However, it should be stressed that what we use in this case toderive weights for our seqlets is a mere 50% of a genome
whereas we test our prediction capability on the remaining
50% of it. The size of the training set we use in this set ofexperiments is much smaller than what has been typically
employed to train previously reported methods. This was anintentional decision meant to showcase our system’s capabilities.
We also carried out experiments without using any a priori
knowledge for the training sets. This is discussed in detail in case
2 below. The purpose of these experiments was to determine how
well our algorithm works in the absence of such information.
The next group of experiments (cases 3a and 3b) was
designed to examine the performance of our method in the casewhere the seqlets’ weights were not genome specific. For eachtest genome, we derived weights for the seqlets by training
with a collection of several complete genomes that did notinclude the genome under consideration; we then applied our
method on this test genome. These jack-knifing experiments
are typically too severe for statistical methods such as those
based on Markov models. It is for this reason that many web
implementations of previously reported, statistics-based
methods often provide several parameter settings derived fromtraining on various genomes: users are asked to select theappropriate settings to be used by the algorithm. The experiments for cases 3a and 3b were carried out using BD-4, i.e. thederived pattern collection was constructed by setting theparameter s to 4 (see section ‘Generation of the Bio-Dictionary’).
In order to determine the impact of the different patterncollections on the results, we repeated the experiments of cases
3a and 3b using BD-6, i.e. the collection constructed by setting theparameter s to 6. Cases 4a and 4b correspond to this set of experiments and are the counterparts of cases 3a and 3b respectively.
Finally, for our last group of experiments (case 5) we usedall 17 genomes to assign weights to the seqlets, then tested theresulting non-genome-specific system by processing each of
the genomes in turn.
As we noted in the Introduction, in this first incarnation of
our system we do not make use of any additional information(e.g. promoter information that can either be computed or
retrieved from the public databases) to constrain the discovery
of genes and start sites. Nonetheless, as evidenced by the
results that we report in cases 4 and 5 (see section ‘Prediction
results on the various genomes’), our start site prediction rates
are comparable in quality, and at times superior to those thathave been previously reported in the literature (5,35–37). Whatis more, our results are achieved using a system that is genericand not genome specific. We will revisit the topic of start-site
prediction at the end of this section.

How we built and used our training and test sets. As we
explained previously, an appropriate training set is needed inFormat of the reported results. In the tables that follow, and foreach processed genome, we report the total number of genes

Quantifying the quality of our predictions. There exist several
ways in which one can evaluate the performance of a gene
finding algorithm. But the algorithm’s sensitivity andspecificity remain the most important measures. Sensitivity,
often referred to as the prediction rate, is defined as the ratio of
the number of genes predicted by the algorithm over thenumber of genes that have been reported in the public
databases. Specificity is defined as the ratio of the number ofpredicted genes that are also reported in the public databases
over the number of all genes that the algorithm has predicted..

2718 Nucleic Acids Research, 2002, Vol. 30, No. 12

predicted by our algorithm in column 2; we also indicate how
many of these genes are <300 nt in column 3. In our studies, wethreshold at that score value for which the number of predictedgenes is equal to the number of annotated genes in the public
databases—clearly this threshold value is different for eachgenome. Note that in this case, the exhibited sensitivity andspecificity are equal; this common value is shown in column 4.
For a subset of the genes that are predicted by our algorithm
there is no corresponding database entry characterizing themas such. Column 5 shows how many of the predicted genes fall
into this category, whereas column 6 indicates how many ofthese genes are <300 nt. With the help of FASTA, BLAST and
CD-search we report how many of these genes are in fact hits
(recall the definition from section ‘Quantifying the quality ofour predictions’) in columns 7, 9 and 11 respectively; thenumber of hits which correspond to gene predictions that are<300 nt is listed in columns 8, 10 and 12 respectively. Finally,in each of the result tables we also report on our ability tocorrectly predict the start sites for the reported genes throughcomparison with the existing database annotations. In column
13, we show the number of genes whose start site is correctlypredicted; and in column 14 we list the same figure as apercentage of reported genes.
Prediction results on the various genomes. We now report on
the results of our method in five experimental settings. The
experiments were designed so as to mimic the type ofsituations that a real-world researcher is likely to encounter,and showcase the performance of our approach across a wide
spectrum of settings.
Case 1. BD-4 and weights derived from the second half ofeach genome only. BD-4 was used in this case. For eachgenome, the seqlets’ weights were obtained by training on thesecond half of the genome. Gene prediction was carried out on
the first half. Table S4 (see Supplementary Material) shows theresults for this experiment: in all 17 cases, the sensitivity/
specificity value ranged between 90.1 and 95.1%. Also, ∼33%of the additional putative genes reported by our method corresponded to hits, i.e. we could find statistically significant similarities with proteins in the September 21, 2001 installment of
the SwissProt/TrEMBL database. The implication of this isthat the actual gene prediction rate is likely to be even higher
than what we report here. With respect to the start site prediction rate, the rates of correctly predicted start sites in this caserange from 55.6 to 84.2%. Recall that we currently make no
use of any promoter information.
Case 2. BD-4 and weights derived from using long ORFs
only. The previous experiment assumes the availability ofannotations for at least some of the actual coding genes of atarget genome. But what if such information is not available?
In such a situation the only recourse is to derive the seqlet
weights by restricting ourselves to the very long ORFs that canbe identified in the genome which is being processed; theimplicit assumption here is that long ORFs are more likely thanshort ORFs to be coding for genes and can thus be used astraining sets. Similar heuristics have been employed by othergroups as well (2,38).
In this case, we used as a training set for the coding regions
all the ORFs that were >600 nt and which were not included
within other longer ORFs. As a training set for the non-coding
regions we used all the ORFs that were <200 nt and whichoccasionally (and incorrectly) include bona fide coding ORFs.
We again used BD-4 as the collection of patterns for which toderive weights. As will become evident after we havedescribed our complete set of experiments, it is not necessary
to carry out this kind of training when dealing with a new
genome—we have simply included case 2 for the purpose ofcompleteness of description.
Table S5 (see Supplementary Material) shows the results ofthis experiment. In this case, the sensitivity and specificity
value ranged from 89.9 to 95.6%. Similarly to case 1, ∼32% of
the additional predicted genes have significant similarities withother database entries. It is notable that despite the fact that theamount of information we used for training purposes was
substantially less than that we used in case 1, the performancelevels remained essentially unchanged.
Case 3a. BD-4 and weights derived from fixed 4 of 17
genomes (leave-many-out). Another realistic situation is theone where users will carry out gene prediction on newly
sequenced genomes for which little or no information is yet
available. The situation can be facilitated if a phylogeneticallysimilar genome already exists in the public databases but this is
not always going to be the case. This particular experiment ismeant to simulate the situation where the genomes that are
already available in the public databases are few and rather
distant from the ones that are being examined. To this end, weused the union of the full CDS lists that were reported for
Bacillus subtilis, Campylobacter jejuni, Helicobacter pylori
and Rickettsia prowazekii to derive the weights for BD-4’s
patterns and subsequently tested our method by predicting
genes for the remaining 13 genomes of our genome collection.
Case 3b. BD-4 and weights derived from 16 of 17 genomes(jack-knifing or leave-one-out). Here we study the jack-knife
variation of case 3a: prior to carrying out gene prediction foreach of the 17 genomes of the collection, we used the union of
the full CDS lists from the remaining 16 genomes to derive theweights for BD-4’s patterns.
Tables S6 (see Supplementary Material) and 4 show theresults for cases 3a and 3b, respectively. As one would intuitively
expect, the weights that are derived from 16 genomes are more
representative (case 3b). Consequently, the prediction rates
reported in Table 4 are superior to those reported in Table S6.
The sensitivity/specificity value was in the range of 86.3–94.7%,
exceeding the 90% mark for the majority of the genomes.
Approximately 23% of the additional putative genes reported
by our algorithm correspond to hits.
This set of experiments (i.e. cases 3a and 3b) is particularlyinteresting in light of observations made previously in the literature according to which the performance of statistical methods
typically deteriorates if the method is used with parameters
derived from a phylogenetically distant genome. This deterioration reflects merely the overall difference in the statisticalproperties of the respective genomes.
Unlike statistical approaches, our method depends largely on
the statistical properties of verified and putative proteins andexhibits resilience to statistical variability. Thus, the conclusion of this set of experiments is that if the weights of theseqlets have been appropriately derived we expect to be able toreach very high prediction levels in a manner that will be
relatively independent of the studied organism. This is indeed
corroborated by the results shown in Tables S6 and 4.

Nucleic Acids Research, 2002, Vol. 30, No. 12 2719

Table 4. Gene prediction results for case 3b

BD-4 was used and the weights were derived from training using the remaining 16 genomes (‘jack-knifing’ or‘leave-one-out’ scenario).
Table 5. Gene prediction results for case 4a

Unlike case 3a, we used BD-6 and the weights were derived from training using only four of the 17 genomes(‘leave-many-out’ scenario). The four genomes that were used for training were: Bacillus subtilis, Campylobacter jejuni,
Helicobacter pylori and Rickettsia prowazekii.

Case 4a. BD-6 and weights derived from fixed 4 of 17
genomes (leave-many-out). In this experiment, we repeat theexperiment of case 3a but now using the BD-6 collection. In
general, with smaller values for s one expects that the derived
collection of patterns will be more sensitive but will captureless ‘structure’. On the other hand, larger s values will give rise
to a situation where potentially fewer seqlets match the putative amino acid translations and there is an associatedincreased difficulty to appropriately compute the seqlets’weights. Table 5 shows the results of this experiment.
Case 4b. BD-6 and weights derived from 16 of 17 genomes(jack-knifing). Here we repeat the experiment of case 3b but
this time using the seqlets from the BD-6 collection: results areshown in Table 6.
When we compare the results from cases 3a/3b with thosefrom cases 4a/4b we conclude that for most of the genomes inour collection, and for the training carried out as describedabove, the BD-4 collection will result in better performance; anotable exception is represented by the Chlamydophila andChlamydia species for which BD-6 gives better results.
An additional observation is that the number of genomes forwhich BD-6 performs better than BD-4 increases as the size ofthe training set increases. The very important ramification ofthis is that, if we have access to a rather large training set, ourmethod will exhibit better prediction performance when usedin conjunction with BD-s sets corresponding to larger s values.
Case 5. BD-6 and weights derived from all 17 genomes. The
four experimental cases above provided sufficient informationthat permitted us to generate optimal results with our method.
In particular, (i) we should use the BD-6 collection since it isbound to perform better than BD-4 as more and more information is deposited in the public databases; and (ii) during2720 Nucleic Acids Research, 2002, Vol. 30, No. 12

Table 6. Gene prediction results for case 4b

Unlike case 3b, we used BD-6 and the weights were derived from training using the remaining 16 genomes (‘jackknifing’ or ‘leave-one-out’ scenario).
Table 7. Gene prediction results for case 5This is the optimal setting where we used BD-6 and the weights were derived from training using all 17 genomes.

training, the weights should be derived from all of the previously published coding/non-coding information for all available genomes. With these two observations at hand, we carryout this last set of experiments noting that it is indicative of thelevels of prediction quality we can expect when we use a goodtraining set, i.e. one derived from many complete genomes.
The results for this experiment are shown in Table 7. The
achieved sensitivity/specificity value is in the 93.9–97.0%
range with most of the genomes exceeding the 95% mark.
Invariably, and similarly to all of the above experiments, a
substantial percentage of the additional putative genes that our
method reports correspond to hits, i.e. they can be corroborated
with the help of sequence similarities to entries in the public
databases.

Another very notable result that is reported in Table 7 has to
do with our ability to correctly predict the start sites of genes.
As can be seen, the ratio of correctly predicted start sites ranges
between 80 and 90% across almost all studied genomes. More
importantly, this ratio is achieved simultaneously with very
high specificity and sensitivity values and without making useof any promoter information.
We conclude by stressing a very important point: our method
derives and uses a single set of weights for the seqlets in the
Bio-Dictionary and these weights are not genome specific.
On the prediction of start sites. As we mentioned already, theaccurate prediction of start codon sites is a notoriously difficult
problem. To derive the various ratios for our experiments, weNucleic Acids Research, 2002, Vol. 30, No. 12 2721

ones from case 5 above and correspond to what we consider tobe an optimal setting. The true start site of the coding region as
well as alternative start sites are shown in this plot. Note how
the local sum of the seqlets’ weights is very low just prior to
the true start codon then ‘jumps’ abruptly to a much higher
value immediately after it. These score jumps can be exploited
to predict the start sites of predicted genes. And as the results
for case 5 have shown, we can achieve high ratios of correct
start site prediction simultaneously with high levels of specificity
and sensitivity. We are in the process of incorporating informationfrom promoter regions to our system and expect that the overall
prediction capability of our method will improve further.

Figure 2. Example of start site prediction using a coding sequence from
E.coli—see text for more details.

made use of the gene starts that are reported in the annotateddatabase entries for each processed genome. Although thiskind of information is generally correct, errors are known toexist. Thus, we also carried out a verification step of our startsite prediction rates using experimentally validated genes. In
particular, we focused on the 1248 experimentally validated
genes from B.subtilis (36) and computed the ratio of start sitesthat were correctly predicted by our algorithm when using thepattern collections and weights described in cases 4b and 5
above [in other words, we used the BD-6 collection of patterns
with weights derived from (i) jack-knifing and (ii) optimal
computations from all 17 genomes]. Of the 1248 genes, we
correctly determined the start sites for 797 (case 4b) and 898
genes (case 5) respectively, or 63.9 and 72.0% of the processed
set respectively. These numbers closely match the corresponding entries for the entire B.subtilis genome in Tables 6(case 4b) and 7 (case 5) respectively; in fact, they are slightly
higher than what is listed in these tables. This verification lends
more support to the correctness of our start site predictions.
We conclude the discussion on start sites by noting that ourgene prediction algorithm can also be used for start site localization. In Figure 2 we show a graph that depicts for eachposition i in the neighborhood of an E.coli coding sequence thelocal sum of the weights for the seqlets that match starting at
position i. To obtain the cumulative score corresponding toposition i, the local sums from i through the stop codon position
must be added together. The seqlet weights used here are the

Scheme 1.

Some notes on the quality of our ‘additional gene’ predictions.
In order to demonstrate that our approach does make a tangible
advance to the gene finding problem, we discuss next two
specific examples.
The first example comes from the processing of the Haemophilus
influenzae genome with BDGF. One of the ORFs that we
predict is found on the reverse strand of the genome in theregion 1476557..1477183; the alleged corresponding geneproduct is 206 amino acids long. Quick analysis shows that thissequence is essentially identical to PSTB_ECOLI, a phosphatetransport, atp-binding protein from E.coli. The following
MUSCA-derived (39) alignment shows the similarity—
conserved amino acids have been colored based on theirhydropathy (Scheme 1).
The ability of our method to correctly identify and report thisparticular gene is very important for the following reason: until
the late part of 2001 the list of genes and their positions in the
genome that accompanied the GenBank entry for H.influenzae
included no gene prediction for this region (recall that the input
database out of which we build our pattern collection dates
back to June 12, 2000) (Scheme 2).
The proteins pstS, pstC, pstA and pstB constitute an intramembrane complex that acts as a high-affinity transportsystem for inorganic phosphates. This system permits thetransport of the phosphates through the inner membrane andinto the cell. The four proteins pstS, pstC, pstA and pstB form
an operon; similarly, phoR and phoB form a second, regulatory
operon. The expression of the genes in these two operons iscontrolled by the level of inorganic phosphates in the cell.
Note that the gene that we predicted is coded by the same
strand (= reverse) as the remaining components of the two
operons. Moreover, the localization and hypothesized identity

2722 Nucleic Acids Research, 2002, Vol. 30, No. 12

Scheme 2.

Scheme 3.

Scheme 4.

Scheme 5.

for our predicted gene are further supported by the fact that
all the members of these two operons are present, appear in
the same order, and have been reported in the genome of
Pastereulla multocida; in fact, the relevant segment of the listof the predicted genes for P.multocida reads as in Scheme 3.
As can be seen, in the P.multocida instance of the two operons,
the gene that codes for pstB is located between the genes
coding for phoB and pstA, similarly to the H.influenzae genome.

Nucleic Acids Research, 2002, Vol. 30, No. 12 2723

Figure 3. Result page of the web-based graphical user interface. The reported ORFs are color coded according to their scores. Selecting an ORF with the cursor
permits the user to see its starting and ending positions, the score that the ORF has received, length and nucleotide composition as well as its amino acid translation.
ORFs can be annotated interactively and regions of the genome can be further explored by zooming in or out using the corresponding buttons (http://
cbcsrv.watson.ibm.com/Tgi.html).

The entry for pstB was added to the list of H.influenzae
genes prior to the October 19, 2001 update of the Entrez entry
for this genome. The value of this example lies in showcasing
our ability to make a correct gene prediction for H.influenzae
as early as June 12, 2000, i.e. the date on which we compiled
the knowledge base that was used to derive our patterns. With
time, the database from which we build our pattern collection
will increase in size and become more rich further improving
our gene prediction ability.
Our second example comes from the genome of A.fulgidus.
One of our gene predictions involves an ORF that appears on
the reverse strand in region 282847..284262; based on the listof genes that were reported with the original GenBank submissionof A.fulgidus, this particular region is gene-empty (Scheme 4).
A straightforward analysis reveals the existence of ahypothetical protein, Q9PEZ3, in Xylella fastidiosa whichshares extensive similarities with our predicted gene. A ClustalW (40) alignment between our predicted gene and Q9PEZ3
is shown in Scheme 5.
The very high degree of similarity and the fact that thecorresponding involved organisms belong to distinct phylogenetic domains (archaea and bacteria respectively) provides
strong support for the correctness of this gene prediction.
Web server for this gene finding algorithm. BDGF, theimplementation of our gene finding algorithm, has been madeavailable via the World Wide Web. The implementation uses
the patterns of the BD-6 collection and optimal weight settings

derived from 17 genomes. The server can be accessed byvisiting the URL http://cbcsrv.watson.ibm.com/Tgi.html and isoperational around the clock; the server runs on a single IBM
RS64III processor with a 450 MHz clock, and can process
250 000 nt in all six reading frames in a little over 60 s. Uponcompletion of the computation, the results are presented to theuser via a graphical user interface an instance of which isshown in Figure 3. The predicted genes are color-coded
depending on the score they have been assigned. The interface
allows the user to navigate around the processed DNA
sequence and to zoom in/out of the regions of interest. Once an
ORF has been selected with the help of the mouse, its locationon the processed sequence is reported together with its nucleotide composition, length and amino acid translation. The
amino acid translation of a selected ORF can be annotatedinteractively. Also, the minimum length of an ORF that will bereported is controlled by the user: its default value is 50 aminoacids (i.e. 150 nt). The complete list of ORFs that have been
predicted by BDGF together with their position, length andassociated score can also be downloaded from the same page.
DISCUSSION AND CONCLUSION
In this paper, we described a new method for solving the gene
identification problem. Our method begins with the
Bio-Dictionary, a collection of patterns that is generated fromprocessing very large public databases with the help of theTeiresias algorithm. The collection accounts completely for2724 Nucleic Acids Research, 2002, Vol. 30, No. 12

the processed input, and discovers genes by making use of thisset of patterns alone. The method is augmented by associatingeach of the used patterns with automatically derived weights.
These weights are genome independent and thus remain fixed
across genomes.
Through a series of carefully designed experiments we
extensively explored various settings that mimicked real-world
situations, and determined the optimal settings for our gene
finding approach. As evidenced by reported experimental
results from 17 archaeal and bacterial genomes, our method
can predict genes very accurately. The method achieves
sensitivity and specificity values that are simultaneously very
high while at the same time achieving a high rate of correctlypredicted start sites. Notably, no promoter or other informationis brought to bear during our determination of the genes and/or
start sites.
We demonstrated the capabilities of our method to extrapolate by intentionally relying upon a Bio-Dictionary that was
built from the June 12, 2000 release of SwissProt/TrEMBL,
i.e. a public collection of sequences that is by now >1.5 years
old. We nonetheless applied the resulting system to genomes
whose ORF translations were included in SwissProt/TrEMBL
either only in part or not at all, with exceptional results.
We should note that in addition to correctly discovering andreporting those ORFs that have already been listed in the
public databases as putative genes, our method determines
additional candidate genes in essentially all of the genomes
that were used in the experiments: for a substantial fraction ofthese previously unreported genes, and with the help of
FASTA, BLAST and CD-search, we determined similarities
with amino acid sequences contained in a very recent release of
SwissProt/TrEMBL. Such similarities further support the
hypothesis that these ORFs ought to have been reported as
putative genes in the first place.
We are currently in the process of pursuing several related
topics that include the determination of the performance
impact using Bio-Dictionaries built from selected collectionsof proteins and not from full-size datasets; the determination of
a compact collection of seqlets for the express purpose of
efficient gene identification; extending our strategy to the caseof eukaryotic genomes, etc.
SUPPLEMENTARY MATERIAL
Supplementary Material is available at NAR Online.
ACKNOWLEDGEMENTS
The authors would like to thank the anonymous reviewers fortheir insightful comments and feedback on the manuscript. The
authors would also like to thank Stephen Chin-Bow and Tien
Huynh for designing and implementing the graphical userinterface of the Web server for BDGF.
REFERENCES
1. Bird,A. (1987) CpG islands as gene markers in the vertebrate nucleus.
Trends Genet., 3, 342–347.
2. Delcher,A.L., Harmon,D., Kasif,S., White,O. and Salzberg,S.L. (1999)
Improved microbial gene identification with GLIMMER. Nucleic Acids
Res., 27, 4636–4641.

3. Salzberg,S.L., Delcher,A.L., Kasif,S. and White,O. (1998) Microbial gene
identification using interpolated Markov models. Nucleic Acids Res., 26,
544–548.
4. Borodovsky,M. and McIninch J. (1993) GeneMark: parallel generecognition for both DNA strands. Comput. Chem., 17, 123–133.
5. Lukashin,A.V. and Borodovsky,M. (1998) GeneMark.hmm: newsolutions for gene identification. Nucleic Acids Res., 26, 1107–1115.
6. Frishman,D., Mironov,A., Mewes,H.-W. and Gelfand,M. (1998)
Combining diverse evidence for gene recognition in completely
sequenced bacterial genomes. Nucleic Acids Res., 26, 2941–2947.
7. Badger,J.H. and Olsen,G.J. (1999) CIRITICA: coding region
identification tool invoking comparative analysis. Mol. Biol. Evol., 16,
512–524.
8. Bafna,V. and Huson,D.H. (2000) The conserved exon method for gene
finding. Proc. Int. Conf. Intell. Syst. Mol. Biol., 8, 3–12
9. Gelfand,M.S., Mironov,A.A. and Pevzner,P. (1996) Gene recognition via
spliced alignment. Proc. Natl Acad. Sci. USA, 93, 9061–9066.
10. Gish,W. and States,D.J. (1993) Identification of protein coding regions by
database similarity search. Nature Genet., 3, 266–272.
11. Robinson,K., Gilbert.W. and Church.G. (1994) Large-scale bacterial genediscovery by similarity search. Nature Genet., 7, 205–214.
12. Burge,C. and Karlin,S. (1998) Finding the genes in genomic DNA.
Curr. Opin. Struct. Biol., 8, 346–354.
13. Burset,M. and Guigo,R. (1996) Evaluation of gene structure prediction
programs. Genomics, 34, 353–367.
14. Claverie,J.M. (1997) Computational methods for the identification ofgenes in vertebrate genomic sequences. Hum. Mol. Genet., 6, 1735–1744.
15. Claverie,J.M. (1998) Computational methods for exon detection. Mol.
Biotechnol., 10, 27–48.
16. Fickett,J.W. (1996) The gene identification problem: an overview for
developers. Comput. Chem., 20, 103–118.
17. Fickett,J.W. and Hatzigeorgiou,A.G. (1997) Eukaryotic promoter
recognition. Genome Res., 7, 861–878.
18. Kehoe,M.A., Kapur,V., Whatmore,A.M. and Musser,J.M. (1996)
Horizontal gene transfer among group A streptococci: implications forpathogenesis and epidemiology. Trends Microbiol., 4, 436–443.
19. Nielsen,K.M., Bones,A.M., Smalla,K. and van Elsas,J.D. (1998)
Horizontal gene transfer from transgenic plants to terrestrial bacteria—a rare event? FEMS Microbiol. Rev., 22, 79–103.
20. Fleischmann,R.D., Adams,M.D., White,O., Clayton,R.A., Kirkness,E.F.,
Kerlavage,A.R., Bult,C.J., Tomb,J.F., Dougherty,B.A., Merrick,J.M.
et al. (1995) Whole-genome random sequencing and assembly of
Haemophilus influenzae. Science, 269, 496–512.
21. Rigoutsos,I. and Floratos,A. (1998) Combinatorial pattern discovery inbiological sequences: the Teiresias algorithm. Bioinformatics, 14, 55–67.
22. Rigoutsos,I. and Floratos,A. (1998) Motif discovery without alignment orenumeration. Proceedings of the 2nd ACM International Conference on
Computational Molecular Biology (RECOMB’98). ACM Press,
New York, NY, 221–227.
23. Rigoutsos,I., Floratos,A., Ouzounis,C., Gao,Y. and Parida,L. (1999)
Dictionary building via unsupervised hierarchical motif discovery.
Proteins, 37, 264–277.
24. Rigoutsos,I., Floratos,A., Parida,L., Gao,Y. and Platt,D. (2000) The
emergence of pattern discovery techniques in computational biology.
Metab. Eng., 2, 159–177.
25. Rigoutsos,I., Gao,Y., Floratos,A. and Parida,L. (1999) Building
dictionaries of 1D and 3D motifs by mining the unaligned 1D sequencesof 17 archaeal and bacterial genomes. Proc. Int. Conf. Intell. Syst. Mol.
Biol., 223–233.
26. Floratos,A., Rigoutsos,I., Parida,L. and Gao,Y. (1999) Sequencehomology detection through large scale pattern discovery. Proceedings ofthe 3rd ACM International Conference on Computational Molecular
Biology (RECOMB’99). ACM Press, New York, NY, pp.164–169.
27. Bairoch,A. and Apweiler,R. (2000) The SWISS-PROT protein sequence
database and its supplement TrEMBL in 2000. Nucleic Acids Res., 28,
45–48.
28. Hofmann,K., Bucher,P., Falquet,L. and Bairoch,A. (1999) The PROSITE
database, its status in 1999. Nucleic Acids Res., 27, 215–219.
29. Bateman,A., Birney,E., Durbin,R., Eddy,S.R., Howe,K.L. andSonnhammer,E.L. (2000) The Pfam protein families database. Nucleic
Acids Res., 28, 263–266.
30. Attwood,T.K., Flower,D.R., Lewis,A.P., Mabey,J.E., Morgan,S.R.,
Scordis,P., Selley.J.N. and Wright,W. (1999) PRINTS prepares for the
new millennium. Nucleic Acids Res., 27, 220–225.

Nucleic Acids Research, 2002, Vol. 30, No. 12 2725

31. Gusfield,D. (1997) Algorithms on Strings, Trees and Sequences:
Computer Science and Computational Biology. Cambridge University
Press, Cambridge, UK.
32. Pearson,W.R. and Lipman,D.J. (1988) Improved tools for biologicalsequence comparison. Proc. Natl Acad. Sci. USA, 85, 2444–2448.
33. Altschul,S.F., Madden,T.L., Schaffer,A.A., Zhang,J., Zhang,Z., Miller,W.
and Lipman,D.J. (1997) Gapped BLAST and PSI-BLAST: a newgeneration of protein database search programs. Nucleic Acids Res., 25,
3389–3402.
34. Marchler-Bauer,A., Panchenko,A.R., Shoemaker,B.A., Thiessen,P.A.,
Geer,L.Y. and Bryant,S.H. (2002) CDD: a database of conserved domainalignments with links to domain three-dimensional structure. Nucleic
Acids Res., 30, 281–283.
35. Frishman,D., Mironov,A. and Gelfand,M. (1999) Starts of bacterial genes:
estimating the reliability of computer predictions, Gene, 234, 257–265.

36. Hannenhalli,S.S., Hayes,W.S., Hatzigeorgiou,A.G. and Fickett,J.W.
(1999) Bacterial start site prediction. Nucleic Acids Res., 27, 3577–3582.
37. Shmatkov,A.M., Melikyan,A.A., Chernousko,F.L. and Borodovsky,M.
(1999) Finding prokaryotic genes by the ‘frame-by-frame’ algorithm:
targeting gene starts and overlapping genes. Bioinformatics, 15, 874–886.
38. Audic,S. and Claverie,J.-M. (1998) Self-identification of protein-coding
regions in microbial genomes. Proc. Natl Acad. Sci. USA, 95,
10026–10031.
39. Parida,L., Floratos,A. and Rigoutsos,I. (1999) An approximation
algorithm for alignment of multiple sequences using motif discovery.
Journal of Combinatorial Optimization, 3, 247–275.
40. Thompson,J.D., Higgins,D.G. and Gibson,T.J. (1994) CLUSTAL W:
improving the sensitivity of progressive multiple sequence alignment
through sequence weighting, position-specific gap penalties and weightmatrix choice. Nucleic Acids Res., 22, 4673–4680.

2726 Nucleic Acids Research, 2002, Vol. 30, No. 12