BACKGROUND
at the time of this writing nearly  <dig> complete prokaryotic genomes, including  <dig> archaeal ones, have been deposited in the genbank database.
driven by the acceleration in genome sequencing, several successful gene prediction programs have been designed for prokaryotic genomes, such as the genemark series  <cit> , glimmer  <cit>  and zcurve  <cit> .
the first two algorithms employ an inhomogeneous markov model for short dna segments , from which an estimate of the likelihood for the segment to belong to a protein coding sequence is derived after a training with existing gene data.
zcurve is based on a z-curve representation of a segment of dna sequence, which is a specific statistic of the whole sequence.
this category of methods are capable of an ab initio prediction of genes for newly sequenced genomes, and they integrate the information from sequence statistics with signal identification  <cit> .
another broad category of gene prediction methods developed earlier are based on similarity search, which comprises such programs as blastx and fasta  <cit> .
orpheus  <cit>  is a typical system that utilizes similarity-based algorithms.
larsen and krogh  <cit>  have developed an easygene system that combines the two approaches, using the classical hidden markov model  in combination with protein similarity search blastp.
however, not all of genes  in newly sequenced prokaryotic genomes show significant similarity with function-known genes, it is still a practical goal to develop the ab initio gene-finding algorithms that require no prior knowledge of the genes of the prokaryotic genome of interest  <cit> .
despite the success of the existing prokaryotic gene finders, systematic bias is significantly due to the necessity of pre-training the prediction program with existing gene data  <cit> .
this is particularly serious for gene starts and short genes whose computational annotation is still quite suspicious  <cit> .
the systematic bias is also notable for gc-rich genomes to be annotated, since the sequence patterns of gc-rich genomes appear to be different to those with lower gc content  <cit> .
it is argued that the number of potential errors in the annotation may be much higher than what is usually believed  <cit> .
a serious case is the archaeal genome aeropyrum pernix.
initially, all orfs longer than  <dig> bps were annotated as coding genes by the original authors who submitted the data to genbank, but significant disagreements arose later from several computational prediction groups  <cit> .
moreover, both recent experiments and in silico analyses have shown that the genomic patterns and hence the mechanisms of translation initiation process are proved to be diversified in archaea  <cit> .
thus the prediction of archaeal genes is far from being a solved problem.
it is commonly believed that a better understanding of the structures of genes both in bacterial and archaeal genomes with a wider variety of gc content is an essential component for improvements of gene prediction.
in this paper, we present a new non-supervised gene prediction algorithm for bacterial and archaeal genomes.
the algorithm aims to develop a comprehensive statistical model with a clear picture of the architecture of prokaryotic genes, in which the biological understanding is explicitly presented.
it is based on the model of protein coding open reading frames  and translation initiation sites .
the former is based on a linguistic "entropy density profile"  model  <cit>  of coding dna sequence and the latter comprises several relevant features related to the initiation of a translation process.
they are combined to form a so-called multivariate entropy distance  algorithm, med  <dig> , that incorporates several strategies in the iterative program.
the iterations enable an non-supervised learning process and to obtain a set of genome-specific parameters before making the prediction of genes.
the main advantage of the algorithm is that it gives efficient and accurate prediction of genes particularly for genomic gc-rich species and archaea, and of gene-related properties such as the divergent translational initiation signals and genome-specific usage of start codon atg, gtg and ttg without any training data.
this feature of our work, namely a close link between the gene prediction method and the biological understanding of gene structure, is believed to be helpful in understanding the genomic comparison and evolution for bacteria and archaea.
the model developed here includes an edp model for the coding potential of an orf, and a tis model for the gene starts.
we begin by briefly reviewing edp model used in our previous work  <cit> .
we then emphasize on describing the tis model, and the implementation of the gene prediction system med  <dig>  based on an iterative learning algorithm.
the edp model for coding potential of orf
the edp model is a global statistical description for a dna sequence, which employs a shannon's artificial linguistic description for a dna sequence of finite length like an orf.
it is different from the markov model-based local approach that usually selects a few local features of the class of sequences.
instead of the amino acid composition {pi}  for an orf, an edp vector s = {si} inferred from {pi} is used to represent the sequence with an emphasis on the information content, where i is the index of the twenty amino acids.
the edp {si} is defined by  <cit> .
si=−1hpilog⁡pi,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbwcdawgaawcbagaemyaakgabeaakiabg2da9iabgkhitmaalaaabagaegymaedabagaemisageaaiabdchawnaabaaaleaacqwgpbqaaeqaaogagiibawmaei4ba8maei4zacmaemicaa3aasbaasqaaiabdmgapbqabagccqggsaalcawljagaaczcamaabmaabagaegymaedacagloagaayzkaaaaaa@4267@
where h=−∑j=120pjlog⁡pj
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgibascqgh9aqpcqghsisldaaewaqaaiabdchawnaabaaaleaacqwgqbgaaeqaaaqaaiabdqgaqjabg2da9iabigdaxaqaaiabikdayiabicdawaqdcqghris5aogagiibawmaei4ba8maei4zacmaemicaa3aasbaasqaaiabdqgaqbqabaaaaa@40ee@ is the shannon entropy, and j =  <dig>  ...,  <dig> also represents the index of the twenty amino acids.
for a set of sequences such as all coding and noncoding orfs from a given genome, the vectors {si} scatter widely in the 20-dimensional edp phase space.
our basic working hypothesis is that the edp vectors for coding orfs form separate clusters from those for noncoding orfs, for a possible reason that they moved along separate paths during the evolution due to different selection pressure.
if this hypothesis is true, one can then develop a clustering method in the edp phase space that discriminates the coding orfs from the noncoding ones.
the success of the method confirms the validity of this hypothesis.
our previous study shows that there exists indeed one universal center representative of coding orfs and one universal center of noncoding orfs in the edp phase space for dozens of analysed species  <cit> .
for gc-rich genomes, our current research based on principal component analysis  shows that orfs form six clusters in the edp phase space due to the genomic gc content bias , one for coding orfs while other five for non-coding orfs, which is also consistent with a previous study  <cit> .
thus the noncoding orfs of gc-rich genomes are represented with five universal centers in the edp phase space.
following the zcurve algorithm  <cit> , we refer to the gc-rich genome as a genome with gc content higher than 56%.
consequently, the universal centers allow us to construct an iterative procedure to define a refined coding potential  <cit> .
this procedure begins by finding a set of "root" coding and noncoding orfs coding} and {snoncoding}), which have been demonstrated with high reliability based on the universal edp centers  <cit> , and then builds a coding potential of a dna sequence s like an orf as the ratio of the distance to each of the edp clusters
Φedp=min⁡kdcoding)min⁡ldnoncoding),
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaahaawcbeqaaiabdweafjabdseaejabdcfaqbaakiabcicaoiabdofatjabcmcapiabg2da9maalaaabagagiyba0maeiyaakmaeioba42aasbaasqaaiabdugarbqabagccqwgebardaqadaqaaiabdofatjabcycasiabdofatjabcicaoiabdugarjabcmcapmaacaaaleqabagaem4yammaem4ba8maemizaqmaemyaakmaemoba4maem4zacgaaagccagloagaayzkaaaabagagiyba0maeiyaakmaeioba42aasbaasqaaiabdygasbqabagccqwgebardaqadaqaaiabdofatjabcycasiabdofatjabcicaoiabdygasjabcmcapmaacaaaleqabagaemoba4maem4ba8maemoba4maem4yammaem4ba8maemizaqmaemyaakmaemoba4maem4zacgaaagccagloagaayzkaaaaaiabcycasiaaxmaacawljawaaewaaeaacqaiyagmaiaawicacaglpaaaaaa@6c5f@
where k and l are the indices of two clusters of root coding and noncoding orfs, d can be defined by the euclidean distance between two edp vectors s = {si} and s' = {s′i
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgzbwcgaqbamaabaaaleaacqwgpbqaaeqaaaaa@2fae@} as
d = 2)1/ <dig>  i =  <dig>  ...,  <dig>
a noteworthy feature of this coding potential is that it is able to generate a cluster with non-trivial boundary and thus achieve a classification of the two classes of sequences with remarkable accuracy  <cit> .
for this reason, the edp model has been successfully used to produce training orf set in the newest version of glimmer method  <cit> .
the tis model for translation start site of prokaryotic genes
the tis model in med  <dig> , which is a further extension of the rbs model described earlier  <cit> , makes a conscious effort to model the following prominent features for prokaryotic gene structure:  upstream to the tis there exists one or a few common motifs that act as the binding sites for the initiation of translation or transcription, these common motifs usually have a position-specific distribution in upstream region from the start sites;  among various genomes, the usage of start codon atg, gtg and ttg is different;  for orfs with multiple candidate start codons, the leftmost start codons are usually the tiss.
however, those start codons located after the leftmost ones are also chosen to be tiss according to different probability;  the  distance between the neighbor genes has a characteristic probability distribution;  for gc-rich genomes, there exists a high gc content at the codon first and third position due to the gc content bias  <cit> .
these features are included in the tis model, allowing us to iteratively find the true tis after the determination of the right 3' end coding orfs.
it is notable that no training data is required to fix the model parameters.
the various probabilities mentioned above are obtained by non-supervised learning during the iteration.
it is intriguing to note also that these probabilities are genome-specific quantities and their differentiation between genomes may likely be interesting and significant for comparative genomic studies.
for the rbs model in our previous work  <cit> , two likelihood functions are defined, one being associated with the shine-dalgarno  motif  <cit> , Φmotif, and another being directly related to the start codon, Φstart.
the total likelihood reads Φtotal = Φmotif + Φstart.
we have made several crucial extensions of that model as follows.
first, many archaea  are detected to have both sd signals of bacterial type in near upstream region and box a motifs of eukaryotic type in the further upstream region from tis.
so, med  <dig>  defines two motif searching regions, vup =  and vlow = , and carries out search for candidate motifs separately in each region.
secondly, a new likelihood function is defined to describe the probability that two orfs overlap, which is given by a distribution of the distance between a start codon and the stop codon of the immediate upstream orf.
thirdly, for gc-rich genomes we use a parameter cp defined by formula  to describe the difference between the coding region vc =  <cit>  downstream and the noncoding region vnc =  upstream to tis.
it has long been known that the sequence patterns of genomes with gc-rich content show differences to those with normal or low gc content.
as a result, high gc content at the codon third and first position has been reported for gc-rich genomes  <cit> .
similar to nishi et al.
<cit> , we refer to gc as the g+c occurrence at the ith codon position, where i =  <dig>   <dig> or  <dig>  then parameter θ=gc+gcgc+gc+gc
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf4oqccqgh9aqpdawcaaqaaiabdeeahjabdoeadjabcicaoiabigdaxiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabiodaziabcmcapaqaaiabdeeahjabdoeadjabcicaoiabigdaxiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabikdayiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabiodaziabcmcapaaaaaa@4a17@ is used to describe the nucleotide usage at the codon first and third position.
for both two regions vc and vnc, we treat them as a set of codon sequences regardless of the fact that the vnc are noncoding sequences.
our analysis shows that the probability distribution of usage pc for coding region vc is significantly different from that of usage pnc of noncoding region vnc.
thus for a dna sequence with a usage θ, the cp parameter is designed to describe the probability of the sequence belonging to coding region as follow
cp = pc/ + pnc),
where  <dig> ≤ cp ≤  <dig>  to evaluate the likelihood of a start codon as tis, we calculate the nucleotide usage θc and θnc for its downstream region  <cit>  and upstream region  separately.
thus the likelihood function is defined by the formula
Φgcstart
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ =  * cp,
where  <dig> ≤ Φgcstart
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ ≤  <dig>  note that the score Φgcstart
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ is calculated only for the genomes with gc content higher than 56%, and this item is added into the start codon likelihood function Φstart developed in our previous work  <cit> .
fourthly, med  <dig>  introduces a measure of signal strength in order to counter-balance the genomic nucleotide content bias.
let t be the ensemble of all l-mers in the set v of upstream sequences of all candidate coding orfs, and for each l-mer t ∈ t, denote by p^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqcaaaa@2e25@ the observed frequency in the region v. denote also by p¯
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ the frequency of the same l-mer if it is a random combination of its nucleotides; thus, p¯
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ can be estimated from single nucleotide composition measured from the whole genome.
similarly to fuglsang  <cit> , we refer to x =  - p¯
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@)2/p¯
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ as a measure of the over-representation of t ∈ t. when summing over the ensemble t, we have a measure, χ <dig> for deciding whether tupor tlow is more likely to contain a strong signal in vup or vlow, where tup  means the ensemble of l-mers in the region vup .
such measure is also used to refer to a weight for each region, as mup = χ2/ + χ2), and mlow = χ2/ + χ2) which allows us to obtain a combined measure for a given candidate start codon
Φtis = mupΦup + mlowΦlow + mstartΦstart,
here, mstart = max.
in above formula, score Φup reflects the likelihood of a motif associated with transcription initiation signal in region vup = , while Φlow means the likelihood of a motif associated with translation initiation signal in region vlow = , both two scores include a weight matrix scoring function and an occurrence probability in upstream region to tis.
for score Φstart, it includes weight matrix scoring function of sequences around tis, the probability of start codon as tis compared with the leftmost start codon, and the orf overlap scoring function.
if the genomic gc content is higher than 56%, the scoring function given by formula  will be added in.
note that the new likelihood function Φtis is able to describe both bacterial and archaeal genomes in a unified way.
finally, a distinction is made between transcript unit internal  genes and transcript unit leader  genes.
following tolstrup et al.
<cit>  and torarinsson et al.
<cit> , we refer to tui genes as those with a start codon separated from the stop codon of the nearest upstream gene by <  <dig> bps, and tul genes otherwise.
our analysis in s. solfataricus shows that only sd-like motifs are found for tui genes while box a motifs are usually found for tul genes.
in med  <dig> , each candidate start codon then definitely belongs to one of the two above categories, each of which conducts the search for motifs separately.
furthermore, med  <dig>  forms a measure
λ=−χ2)−χ2),
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf7oabcqgh9aqpcqggoaakcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgjbqsaeaacqwg1bqdcqwgwbacaagccqggpaqkcqghsislcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgjbqsaeaacqwgsbabcqwgvbwbcqwg3bwdaagccqggpaqkcqggpaqkcqggoaakcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgmbataeaacqwg1bqdcqwgwbacaagccqggpaqkcqghsislcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgmbataeaacqwgsbabcqwgvbwbcqwg3bwdaagccqggpaqkcqggpaqkcqggsaalcawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@6d64@
to decide which of the two translation initiation mechanisms is predominant in a genome: λ >  <dig> for either completely sd signals or completely box a consensus, and λ ≤  <dig> for a mixture of the two.
tests show that most of bacteria and part of archaea belong to the former, whereas many of archaea belong to the latter.
in the latter case, med  <dig>  uses the z-score to form a normalized tis score which allows to compare and select the right tis.
therefore, our method allows to describe more than one type of initiation signals which are located at different upstream positions.
the iterative learning algorithm
the gene finding system med  <dig>  predicts genes in two stages: the coding orf detection stage and the tis refinement stage.
the outline of sequence processing and model self-training for two stages are shown in fig.
<dig>  we give a summary of this algorithm as follows.
in the first step of coding orf detection stage, all orfs longer than  <dig> bps are extracted in both strands and in all three reading frames, thus to determine the edp coding potential.
as stated above, using the universal coding and noncoding edp centers and with an iteration of consecutive root orf finding, we decide a set of root coding/noncoding orfs for each genome being analysed, which typically cover more than 60% of all genes with a very high reliability over 99%  <cit> .
note that we obtain both "seed" coding orfs and non-coding sequences, instead of generating artificial non-coding orfs as other methods do  <cit> .
these seed sequences form a reliable learning set for the further analysis.
the system then goes to determine the remaining orfs excluded from the root orfs.
we search for conservative motifs only in the near upstream region vlow of the root coding orfs, namely setting mup =  <dig>  and mlow = mstart =  <dig> in calculating the tis score in formula .
the candidate start site with the highest tis score is selected, then the parameters are recalculated using the updated start sites.
this iteration process ends when no start site needs to be relocated.
the parameters are then obtained and used to determine the remaining orfs.
this consists in calculating the tis score of the remaining orfs, and then evaluating their coding potential using a fisher discriminant algorithm  in a two-dimensional plane of the edp coding score and the tis score.
for gc-rich genomes, the fisher discriminant algorithm is applied in a three-dimensional space of the edp coding score, the tis score and the θ score , which is trained using p. aeruginosa.
our studies show that the discriminant coefficients are universal across prokaryotic genomes.
thus the remaining orfs are classified into coding orfs and noncoding ones.
to reduce the false positive in the detected coding orfs determined from the remaining orfs, a strategy to resolve overlap of coding orfs is then applied.
we define a quantity ξ to be the percentage of nucleotides overlapped with other coding orfs in either the same or the opposite strand.
the iteration goes as follows:  extend all the detected coding orfs to the longest orfs by assigning the leftmost start codon as the start of each orf;  remove those orfs with ξ > ξi ;  recalculate the overlapping percentage for all orfs retained from the last step; and  reduce ξi by 1%, then repeats from step .
this procedure iterates until ξ reaches a threshold .
after the stage of coding orf detection, med  <dig>  goes to the next stage of tis refinement for all coding orfs.
the detected coding orfs are first classified into two sets, tui and tul.
the system considers searching for either sd motifs or box a motifs in both regions vup =  and vlow = .
then, one forms the combined tis score in formula , which is calculated together with all three coefficients mup, mlow and mstart described above.
for gc-rich genomes, the score Φstart includes the item of Φgcstart
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ in formula .
then the global genomic measure λ is used to decide which of the initiation mechanisms is more appropriate for the genome studied.
if λ >  <dig>  the self-training of tis model parameters runs in both vup and vlow for all detected coding orfs, we select the sd-like motif as the main mechanism if χ <dig> <χ <dig>  or the box a motif otherwise.
for each coding orf, the tis score in formula  is calculated for all candidate start codons, then the start codon with the highest score is selected as the predicted tis.
if λ ≤  <dig>  the self-training is run with vup and vlow together, but for tui and tul genes independently.
for each coding orf, the normalized scores of both tui type and tul type are calculated, and then compared to select one with higher score as the most likely start codon to relocate the gene start.
RESULTS
genome sequences and reliable gene datasets
the bacterial and archaeal genomes and their annotations used in this paper were downloaded from the genbank release  <dig>  in  <dig>
we have selected experimentally confirmed genes as benchmarks.
for e. coli, we have used two datasets, the ecogene  <cit>  and the link dataset  <cit> , with their newest version .
for b. subtilis, following the standard by besemer et al.
<cit> , we have taken it as one of a few annotated complete genomes that can be used to evaluate the performance of exact gene prediction.
thus we denoted all  <dig>  genes annotated in genbank as bsub_all dataset.
for verifying short genes, we have chosen bsub <dig>  bsub <dig> and bsub <dig> datasets, which are selected from the b. subtilis and verified by protein similarity search  <cit> .
we have extracted  <dig> short genes from the  <dig> genes of ecogene, denoted as ecogene_short, to evaluate short gene prediction for e. coli.
for gc-rich genomes, we have built two data sets, mtub <dig> and paer <dig>  the former includes  <dig> reliable genes with confirmed tiss in m. tuberculosis with gc content of  <dig> %, while the latter has  <dig> such genes in p. aeruginosa with gc content of  <dig> %.
for archaeal genomes, we have built a data set solfgene including  <dig> reliable genes in s. solfataricus with confirmed tiss.
the three data sets mtub <dig>  paer <dig> and solfgene, as well as the process of determining these genes with their tiss confirmed by n-terminal protein sequencing or inferred from experimental evidences, can be accessed through  <cit> .
other benchmark programs
in order to benchmark med  <dig> , we have tested four gene-finding methods , glimmer , genemarks  <cit> , zcurve  <cit>  and easygene  <cit> , to make a comparison of the performance of gene prediction.
the former three methods belong to the same category of our system as the ab initio gene prediction methods, while the easygene belongs to another category of methods using the extrinsic information by means of similarity search  <cit> .
glimmer  <dig>  was downloaded from  <cit>  and installed locally.
we ran it following the instructions given in the distribution file.
a post-processor rbsfinder  <cit>  has been designed to further improve the tis prediction by glimmer  <dig> .
thus, rbsfinder was used to process the original output of glimmer  <dig>  and the refined tiss are then taken as predicted ones by glimmer system.
rbsfinder is accessible via the website at  <cit> .
herein, rbsfinder was run repeatedly until over 99% of gene starts remain unchanged.
recently, the newest version of glimmer, i.e.
v <dig> , has been released with several algorithmic changes to reduce the number of false positive predictions and to improve the accuracy of tis predictions.
although the article describing glimmer  <dig>  has not yet been found from publication, a locally executable program is current available from  <cit> .
in order to have a comprehensive comparison with glimmer method, glimmer  <dig>  was further included herein as benchmark program.
we downloaded the predictions by glimmer  <dig>  from its website for all genomes studied in this paper.
the software zcurve  <dig>  is freely available at  <cit> .
it was executed on our pc following the instruction file.
genemarks provides only an online service of genomic annotation, instead of a local executable program.
the newest results on several test sets mentioned above for e. coli and b. subtilis by the current genemarks have been greatly improved than that of its initial version  <cit> .
similarly, there is still no new literature to report the improved version of genemarks.
however we adopted all the unpublished results returned to us via email with the current online analysis behind the interface  <cit> , to compare with our program.
similar to genemarks, easygene is only accessible via its web interface, and the pre-trained models are available for only  <dig> genomes  <cit> .
however, the easygene web server has provided predictions for hundreds of chromosomes.
therefore, for comparison purpose, easygene predictions for a total of  <dig> genomes were downloaded from its website  <cit>  available at the time we prepared our paper.
accuracy of gene detection and interpretation of genome-specific model parameters
to illustrate the prediction accuracy of med  <dig> , we first present the comparison against the genbank annotation for the 3' end match.
two independent quantities, sn  and sp , are defined to evaluate the performance of a gene finder at gene level as:
sn = tp/, sp = tp/.
here, tp, fp, and fn are the number of true positive, false positive, and false negative, respectively.
the med  <dig>  program has been run on all the complete bacterial and archaeal genomes currently available on genbank.
additional file  <dig> includes the prediction results on forty representative genomes, comprising  <dig> bacteria and  <dig> archaea.
selected organisms cover most of the taxonomic groups, including crenarchaeota, euryarchaeota, nanoarchaeota for archaea and α/β/γ/ε-proteobacteria, actinobacteria, bacteroidetes/chlorobi, for bacteria  <cit> .
for the forty species, med  <dig>  achieves an average sensitivity of  <dig> % and specificity of  <dig> %.
the results of glimmer  <dig> , glimmer  <dig> , zcurve  <dig>  and the newest genemarks for the forty genomes as a comparison are not shown in additional file  <dig>  but may be found in our website.
on the average, the newest version of genemarks gives the same level of sensitivity  and the highest specificity , while glimmer  <dig>  and zcurve  <dig>  give the sensitivity of  <dig> % and  <dig> %, and the specificity of  <dig> % and  <dig> %, respectively.
as reported by easygene's website, the average sensitivity and specificity are  <dig> % and  <dig> % respectively for  <dig> of the  <dig> genomes.
compared with glimmer  <dig>  on average, glimmer  <dig>  significantly reduces the false positives  with almost no decreasing of the sensitivity of  <dig> %.
the average sensitivity and specificity of genemarks reported in the original version  <cit> , which was run on eight genomes including e. coli, b. subtilis and a. fulgidus, are  <dig> % and  <dig> %, respectively.
a full comparison for  <dig> genomes has also been listed at our web page.
for the  <dig> genomes, med  <dig>  gives the average sensitivity of  <dig> % and specificity of  <dig> %, while glimmer  <dig>  gives  <dig> % and  <dig> %, zcurve  <dig>  gives  <dig> % and  <dig> %, glimmer  <dig>  gives  <dig> % and  <dig> %, genemarks gives  <dig> % and  <dig> %, and finally easygene gives  <dig> % and  <dig> % averaged over  <dig> genomes, respectively.
however, we argue that the higher specificity of the latter three programs is not unrelated to the fact that they are highly optimized with the genbank annotation.
it is clear that med  <dig>  is competitive with glimmer  <dig>  and zcurve  <dig> , while the specificity of prediction is a little higher than both of them.
both glimmer  <dig>  and the newest genemarks present rather high accuracies of sensitivity and specificity against the current genbank annotation for the 3' end match, although there are not yet publications to report them.
using the extrinsic information of similarity search, easygene is able to extremely raise the specificity at the cost of its sensitivity slightly lower than the ab initio gene prediction methods.
it should be noted that easygene has included many a confirmed gene of the query-genome as training set by using blastp to search for significant protein matches in swiss-prot  <cit> , which makes an essential distinction between such a method and the ab initio gene prediction method such as med  <dig> .
for instance, there are about 65% of all  <dig>  genes in e. coli predicted by easygene have been taken as the training set.
therefore a comparison against all genes in the benchmark would over-estimate easygene's prediction performance.
to discuss the predicting performance on the current genbank annotation, we should point out that it can not be free from the bias since both glimmer and genemark series have been widely used or involved in the genbank annotation pipeline.
the statistical analysis shows that there are over  <dig> genomes from genbank before  <dig> are in this case.
while the comparison against the genbank annotation has played a role during the early development of the computational prediction programs, this comparison becomes increasingly suspicious as its bias becomes evident .
since the genbank annotation is not fully accurate, further evaluation is performed based on the function-known genes which have more reliable annotations.
for each genome, the function-known genes are selected from genbank by excluding those with product descriptions with any of the key words as "-like", "conserved", "hypothetical", "homolog", "probable", "possible", "predicted", "putative", "similarity" and "unknown".
the ratio of the number of these genes predicted correctly by med  <dig>  is also listed in additional file  <dig>  note that the specificity is meaningless in this case.
as we can see, the average accuracy for med  <dig>  is  <dig> %, while the average accuracy for glimmer  <dig> , glimmer  <dig> , zcurve  <dig> , easygene  and the newest genemarks is  <dig> %,  <dig> %,  <dig> %,  <dig> % and  <dig> % , respectively.
for  <dig> genomes analysed, the average accuracy for med  <dig> , glimmer  <dig> , glimmer  <dig> , zcurve  <dig> , easygene and the newest genemarks may be calculated, which gives  <dig> %,  <dig> %,  <dig> %,  <dig> %,  <dig> % and  <dig> % respectively, clearly, they are all at the same level.
these results can be found on our website.
therefore the gene-finding accuracy for function-known genes of med  <dig>  and other five gene finders is well matched.
a distinctive feature of med <dig>  is that the algorithm builds a prokaryotic gene structure model based on a set of genome-specific parameters, which are calculated by several self-learning iterations without any prior knowledge or training data.
comparing the model parameters among different species, e.g.
those of gene-related properties as tis-upstream signals and usage of start codons atg, gtg and ttg, would shed light on the study of genomic comparison and evolution for bacteria and archaea  <cit> .
to this end, we present in additional file  <dig> the parameters of start codon usage and motifs found upstream to tis for each species, although a more detailed list is provided in the output file by our program for the user.
genomes in additional file  <dig> are listed as two groups of bacteria  and archaea  and sorted alphabetically.
almost for all analysed species, atg is the most common start codon  compared with ttg and gtg.
an exception is the m. kandleri genome, for which the three start codons are used nearly equally .
in contrast, the usage of ttg to gtg is different for different species.
for example, gtg in c. glutamicum is much more favored than ttg ; while ttg in b. burgdorferi is more favored than gtg .
for  <dig> archaeal genomes, the start codon usage predicted by med  <dig>  is very close to an earlier result annotated by easygene including the protein similarity search method  <cit> .
with the automatic control coefficient λ in formula  and the resulting predicted motifs obtained by a non-supervised learning process, med  <dig>  clearly defines three types of translation initiation mechanisms and ascertains which of the types each genome belongs to.
in additional file  <dig>  we go further and report the predicted motifs, meaning the most significant consensuses associated with translation initiation signals in two searching regions  and  upstream to tiss.
for all of  <dig> analysed bacterial genomes, motifs are presented only in region , while most of them belong to sd signal type, which means the typical mechanism of translation initiation with the sd signals for transcripts with leaders.
however, at-rich motifs, "aattt" and "aaatt" listed in additional file  <dig>  are detected in region  for the synechocystis sp.
genome, indicating the existence of a leadless translation initiation mechanism in cyanobacteria.
a similar case is to the b. thetaiotaomicron genome, where a-rich motifs are detected.
for  <dig> analysed archaeal genomes, three types of mechanisms are detected.
six species m. jannaschii, m. kandleri, m. maripaludis, m.
thermoautotro., p. abyssi and p. furiosus are reported with motifs of sd signal type only in region , meaning the mechanism similar to the bacterial.
in contrast, three genomes halobacterium sp., n. equitans and p. aerophilum are shown that the motifs of box a type are predicted only in region , which suggest a different mechanism of translation initiation with a high level of genes producing leaderless transcripts.
for three genomes a. fulgidus, p. torridus and s. solfataricus, both box a motifs in region  for leaderless transcripts and sd motifs in  for transcripts with leaders are detected.
it implies that a mixed mechanism of leaderless translation and leadered translation is used in these genomes.
the result is well matched with that of previous studies  <cit> .
a notable result is that most of the motifs of searching region  found in bacteria contain the tetramer "ggag" in contrast with that of "ggtg" in most archaea, it is also exactly consistent with that reported in the earlier study  <cit> .
however, note that for the thermotoga maritima genome, the med  <dig>  analysis leads to an interesting finding of sd-like motifs very richly observed in archaea.
a phylogenetic tree based on 16s rrna sequence has revealed that t. maritima is a deep branching species  <cit> .
moreover, it has been reported that over 20% of genes in t. maritima were acquired through horizontal gene transfer from archaea  <cit> .
thus, the archaea-like feature in t. maritima might be interpreted as either an ancient feature or, alternatively, a consequence of horizontal gene transfer.
accuracy of exact prediction on reliable genes
in order to perform an accuracy test on the exact gene prediction, reliable gene datasets must be used since the genbank annotations have systematic bias in the 5' end prediction  <cit> .
before showing the results in this section, we should like to stress that all six gene finders, including med  <dig> , are compared to analyze each genome with the same sensitivity and specificity against the genbank annotation as mentioned above.
as a gene finder using a similarity search by alignment in swiss-prot to build tls-confirmed orfs, many a protein coding gene has been included as training set of tis for easygene.
for example, there are  <dig> of the  <dig>  genes in bsub_all dataset, and  <dig> of the  <dig> genes in ecogene, have been employed as its training set for b. subtilis and e. coli respectively.
thus for a fair comparison with the ab initio methods, the accuracies of easygene on reliable test sets were calculated by excluding the number of genes included into its training set.
first we present the prediction accuracy on reliable test sets bsub_all, ecogene and link for bacterial genomes.
as indicated in table  <dig>  for only 3' end match, all six programs have nearly the same performance with accuracy above  <dig> %.
for both 5' and 3' end matches, the current version of genemarks achieves the highest performance with average accuracy of  <dig> %, while med  <dig>  has a comparative performance with average accuracy of  <dig> %.
note that for bsub_all dataset, genemarks has reported its accuracy of  <dig> % for 3' end match and  <dig> % for both end matches on  <dig>  genes annotated in genbank  <cit> , which is lower than the newest result  by the current genemarks version, as well as lower than that of med  <dig>   .
zcurve  <dig>  shows slightly lower average accuracy of  <dig> % for both end matches, and easygene does  <dig> %, while glimmer  <dig>   gives the lowest average accuracy of  <dig> %.
glimmer  <dig>  has significantly improved the average accuracy to  <dig> %, reaching a similar performance as med  <dig> .
one of the remaining challenges for prokaryotic gene finding is the identification of short genes 'buried' in an enormous pile of false-positive short orfs.
following the standard adopted by many gene finders, we hereby define short genes as those with length between  <dig> and  <dig> bps.
we have carried out a test on  <dig> short genes from  <dig> confirmed ones in ecogene data set .
the test indicates that med  <dig>  has an accuracy of  <dig> % for 3' end match and also a high accuracy of  <dig> % for both end matches.
for the bsub <dig>  bsub <dig> and bsub <dig> datasets  <cit>  from b. subtilis, med  <dig>  detects  <dig> %,  <dig> % and  <dig> % of 3' end match and  <dig> %,  <dig> % and  <dig> % of both end matches, respectively.
similarly, note that the corresponding published results of genemarks on the same three datasets are  <dig> %,  <dig> % and  <dig> % of 3' end match as well as  <dig> %,  <dig> % and  <dig> % of both end matches  <cit> , which have been greatly improved given by the newest online version.
however as listed in table  <dig>  comparing with other five gene finders shows that med  <dig>  gives a higher performance in the short gene prediction than glimmer  <dig>  , glimmer  <dig> , easygene and zcurve  <dig> , while is comparable with genemarks.
aprograms med, gl <dig> , gl <dig> and zcv were run locally, while gmk was run online, as described in the text.
predictions for eg were downloaded from  <cit> .
bexperiment confirmed tiss data sets: the first three represent two well-studied genomes: b. subtilis  and e. coli ; the fourth to seventh represent short genes for e. coli  and b. subtilis ; psaer <dig> and mtub <dig> are selected for two gc rich genomes, m. tuberculosis  and p. aeruginosa ; solfgene corresponds to the archaeal s. solfataricus.
cnumbers in parentheses indicate that the results of genemarks have been reported in literature,  means no data reported.
to evaluate the performance of exact gene prediction for gc-rich genomes, two reliable sets paer <dig> and mtub <dig> are used here for p. aeruginosa with gc content of  <dig> % and m. tuberculosis with gc content of  <dig> %.
test on paer <dig> shows that med  <dig>  has an accuracy of  <dig> % for 3' end match and  <dig> % for both end matches, while on mtub <dig> gives an accuracy of  <dig> % for 3' end match and  <dig> % for both end matches.
as we can see from the eighth and ninth lines in table  <dig>  compared with glimmer  <dig>  , glimmer  <dig> , zcurve  <dig> , easygene and the newest online version of genemarks, med  <dig>  demonstrates the highest performance of exact gene prediction for gc-rich genomes.
it is still problematic to evaluate gene prediction for archaeal genomes due to the insufficient number of genes confirmed by independent  methods  <cit> .
to this end, we have built the data set solfgene comprising  <dig> reliable genes with confirmed tiss.
all six programs detect all genes at the 3' end in solfgene.
as for both ends prediction, both med  <dig>  and easygene give the highest accuracy of  <dig> % vs glimmer  <dig>   , glimmer  <dig>  , genemarks , and zcurve  <dig>   .
in summary, upon test of the reliable gene sets built for two well studied genomes e. coli and b. subtilis, we show that the total prediction performance of med  <dig>  can be matched with that of genemarks in its publication, and higher than that of glimmer, easygene and zcurve  <dig> .
while for the newly analysed genomes, especially for gc-rich genomes and archaea, med  <dig>  shows an evident advantage, because test on the reliable gene sets clearly indicates that med  <dig>  has given a higher accuracy than other five gene finders, at least the same level of the best one.
analysis of upstream region of predicted tiss for archaea
med  <dig>  gives a comprehensive modeling of the tis for both bacterial and archaeal genomes, in particular for archeal genomes in which translation initiation is more complex than in bacteria  <cit> .
we report herein the results using the extended sequence logos  <cit>  for the region  upstream from predicted tiss by med  <dig> .
the sequence logos facilitate detection of sequence patterns that are conserved in both content and position.
we present the logos of three representative archaeal genomes: m. jannaschii, p. abyssi, and n. equitans in fig.
<dig>  the logos indicate that several mechanisms associated with the translation initiation have been detected by med  <dig> .
the genome of m. jannaschii displays the logos with a remarkable sd signals in region  like bacteria, which is well compatible with the notion that this genome produces nearly no leaderless transcripts  <cit> .
on the other hand, the n. equitans genome that carries a high level of genes producing leaderless transcripts  <cit>  produces a very rich set of signals.
an evident a/t-rich signal, the so-called box a promoter motif, is observed near the position - <dig>  in addition, a tfb recognition element  motif constituting  <dig> to  <dig> a/ts, which interact with the archaeal transcription factor tfb, is also recognized about  <dig> bps upstream from the box a motif.
furthermore, an a/t peak is detected near the position - <dig> to - <dig> upstream to the tis, which has been reported being functionally important  <cit> .
finally, the p. abyssi genome is found to have a strongly conserved sd signal in - <dig> to - <dig> upstream region together with a rather weak box a motif upstream from the sd signal.
this is well consistent with the fact that the p. abyssi genome produces fewer leaderless than leadered transcripts  <cit> .
inspection on the upstream sequence patterns of predicted tiss would be helpful for the assessment of the overall quality of gene annotation or prediction.
herein we choose an example, the s. solfataricus genome, for a detailed comparison between med  <dig>  and other prediction programs, as well as the genbank annotation.
first we contrast med  <dig>  with genbank annotation.
as we can calculate from additional file  <dig>  med  <dig>  detects  <dig> genes with common 3' end against genbank annotation, among which  <dig> genes have common tiss and  <dig> genes do not .
fig.
3b shows the logos of the upstream region  from tiss of the common  <dig> genes, where a typical sd signal can be observed at  and an evident box a motif at .
the logos for the other  <dig> tiss predicted by med  <dig>  are shown in fig.
3c, where similar sd signal and box a motif are found, while those tiss from the genbank annotation shown in fig.
3d show no informational structure in the logos.
this comparison suggests that med  <dig>  can give a more relevant annotation of tis for archaeal genomes than the current genbank annotation.
similar comparison is made with glimmer  <dig> , genemarks and zcurve  <dig> , as shown in fig.
<dig>  we have shown that glimmer  <dig>  has achieved an overall improvement of tis prediction than that of its earlier version, thus have no need for the comparison with glimmer  <dig> .
the sets of common tis predictions are slightly different in each case , but in all comparisons the tiss predicted by med  <dig>  show the sequence logos similar to that in fig.
3b , while those of glimmer  <dig>  , genemarks  and of zcurve  <dig>   show almost no sd signal or box a motif structures.
thus to a great extent, med  <dig>  makes a substantially more accurate prediction than glimmer  <dig> , genemarks and zcurve  <dig>  for archaeal genomes.
CONCLUSIONS
in this paper, we present a comprehensive model to describe a set of properties about the coding potential and the translation initiation mechanisms for both bacteria and archaea.
based on the model, with a design of multiple iterations, an non-supervised ab initio gene prediction system med  <dig>  is developed.
generally speaking, the system is able to adapt to any newly sequenced prokaryotic genome with no need for any data training or prior knowledge, and can predict the divergent translation initiation mechanisms and the resulting signals upstream from the tis.
the model seems to be biologically sound, since the program yields genome-specific model parameters such as various probabilities associated with the translation initiation signals and the start codon usage, which are matched with the current knowledge from earlier works.
thus the model may provide a good tool for comparative genomic studies.
upon test of a set of reliable gene sets, the total prediction performance of med  <dig>  for the well studied bacteria such as e. coli and b. subtilis with usual gc content can be matched with that of the existing published methods.
while for the newly analysed genomes, especially for gc-rich genomes and archaea, test on the reliable gene sets indicates that med  <dig>  outperforms or at least gives the same level of the best of the current gene finders in exact gene prediction for both 3' and 5' end matches.
furthermore, med  <dig>  adapts to a broad range of archaeal genome as well as to bacterial genomes.
for archaeal genomes with more complex mechanisms of translational initiation, our method has a more accurate prediction of tiss compared to the existing gene finders and the current genbank annotation.
as many a prokaryotic gene prediction program, such as glimmer, genemark and easygene, has been widely used in raw genomic sequence annotation, the predicting performance for a new gene prediction method evaluated by data extracted from genbank should be in caution.
while the comparison against the genbank annotation has played a role during the early development of the computational prediction programs, this comparison becomes increasingly suspicious as its bias becomes evident .
therefore the consistency with the genbank annotation should not be considered to be the only factor, and should not even be the essential one for some archaeal genomes.
with hundreds of prokaryotic genomes have been sequenced and made publicly available, subsequent studies will focus more on the development of the comprehensive model of gene structure that allow capture of the biological evidence of prokaryotic genomes, e.g.
the translation initiation mechanism, which is essential to understanding comparative genomics and biology evolution.
at present the problems of being biased and erroneous have not yet been solved in the genbank annotation for prokaryotes, especially for newly sequenced archaeal genomes and gc-rich genomes, and the further promotion of comparative genomics, it is hoped that our system is shown to meet the demand and be qualified as an alternative tool of gene computational prediction.
the expert annotators, who are concerned with more accurate and complete annotation of sequenced genomes, should find our method useful when used independently or alongside with other gene finding tools.
availability and requirements of the system
project name: med project;
project homepage: ;
operating systems: microsoft windows  <dig>  operating system is recommended for the pre-complied program downloaded from our web site.
in addition, we also provide the source code on linux with a c++ compiler;
programming language: c++;
licence: the sourc code is freely available  under the gnu gpl license;
any restrictions to use by non-academics: none.
authors' contributions
zss and hqz co-supervised the whole development of the work and had the main idea of the algorithm, gqh and hqz designed and implemented the gene-finder system, drafted most parts of the manuscript.
yfy and jw contributed biological expertise and performed part of the experimental test.
all authors read and approved the final manuscript.
supplementary material
additional file 1
med  <dig>  prediction accuracy and genomic-specific model parameters.
contains data for med prediction accuracy and genomic-specific model parameters for  <dig> bacterial genomes and  <dig> archaeal genomes.
the file is edited by microsoft excel  <dig>
click here for file
additional file 2
source code.
contains source code of med  <dig>  written for windows and linux/unix operation systems.
click here for file
acknowledgements
we thank zhengqing ouyang, yongchu liu, zongxiao sun, xiaobin zheng and others for beneficial discussions and helps to the work.
we thank also j. o'brien for kindly reading and revising the manuscript.
we also thank prof. chun-ting zhang for providing zcurve  <dig>  program.
the work received partial support by the national natural science foundation  of china; the present study was also supported by the national basic research program of china  under grant no.
2003cb <dig>

 at the time of this writing nearly complete prokaryotic genomes, including archaeal ones, have been deposited in the genbank database. driven by the acceleration in genome sequencing, several successful gene prediction programs have been designed for prokaryotic genomes, such as the genemark series, glimmer and zcurve. the first two algorithms employ an inhomogeneous markov model for short dna segments, from which an estimate of the likelihood for the segment to belong to a protein coding sequence is derived after a training with existing gene data. zcurve is based on a z-curve representation of a segment of dna sequence, which is a specific statistic of the whole sequence. this category of methods are capable of an ab initio prediction of genes for newly sequenced genomes, and they integrate the information from sequence statistics with signal identification. another broad category of gene prediction methods developed earlier are based on similarity search, which comprises such programs as blastx and fasta. orpheus is a typical system that utilizes similarity-based algorithms. larsen and krogh have developed an easygene system that combines the two approaches, using the classical hidden markov model in combination with protein similarity search blastp. however, not all of genes in newly sequenced prokaryotic genomes show significant similarity with function-known genes, it is still a practical goal to develop the ab initio gene-finding algorithms that require no prior knowledge of the genes of the prokaryotic genome of interest. despite the success of the existing prokaryotic gene finders, systematic bias is significantly due to the necessity of pre-training the prediction program with existing gene data. this is particularly serious for gene starts and short genes whose computational annotation is still quite suspicious. the systematic bias is also notable for gc-rich genomes to be annotated, since the sequence patterns of gc-rich genomes appear to be different to those with lower gc content. it is argued that the number of potential errors in the annotation may be much higher than what is usually believed. a serious case is the archaeal genome aeropyrum pernix. initially, all orfs longer than bps were annotated as coding genes by the original authors who submitted the data to genbank, but significant disagreements arose later from several computational prediction groups. moreover, both recent experiments and in silico analyses have shown that the genomic patterns and hence the mechanisms of translation initiation process are proved to be diversified in archaea. thus the prediction of archaeal genes is far from being a solved problem. it is commonly believed that a better understanding of the structures of genes both in bacterial and archaeal genomes with a wider variety of gc content is an essential component for improvements of gene prediction. in this paper, we present a new non-supervised gene prediction algorithm for bacterial and archaeal genomes. the algorithm aims to develop a comprehensive statistical model with a clear picture of the architecture of prokaryotic genes, in which the biological understanding is explicitly presented. it is based on the model of protein coding open reading frames and translation initiation sites. the former is based on a linguistic "entropy density profile" model of coding dna sequence and the latter comprises several relevant features related to the initiation of a translation process. they are combined to form a so-called multivariate entropy distance algorithm, med, that incorporates several strategies in the iterative program. the iterations enable an non-supervised learning process and to obtain a set of genome-specific parameters before making the prediction of genes. the main advantage of the algorithm is that it gives efficient and accurate prediction of genes particularly for genomic gc-rich species and archaea, and of gene-related properties such as the divergent translational initiation signals and genome-specific usage of start codon atg, gtg and ttg without any training data. this feature of our work, namely a close link between the gene prediction method and the biological understanding of gene structure, is believed to be helpful in understanding the genomic comparison and evolution for bacteria and archaea. the model developed here includes an edp model for the coding potential of an orf, and a tis model for the gene starts. we begin by briefly reviewing edp model used in our previous work. we then emphasize on describing the tis model, and the implementation of the gene prediction system med based on an iterative learning algorithm. the edp model for coding potential of orf the edp model is a global statistical description for a dna sequence, which employs a shannon's artificial linguistic description for a dna sequence of finite length like an orf. it is different from the markov model-based local approach that usually selects a few local features of the class of sequences. instead of the amino acid composition {pi} for an orf, an edp vector s = {si} inferred from {pi} is used to represent the sequence with an emphasis on the information content, where i is the index of the twenty amino acids. the edp {si} is defined by. si=−1hpilog⁡pi,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbwcdawgaawcbagaemyaakgabeaakiabg2da9iabgkhitmaalaaabagaegymaedabagaemisageaaiabdchawnaabaaaleaacqwgpbqaaeqaaogagiibawmaei4ba8maei4zacmaemicaa3aasbaasqaaiabdmgapbqabagccqggsaalcawljagaaczcamaabmaabagaegymaedacagloagaayzkaaaaaa@4267@ where h=−∑j=120pjlog⁡pj mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgibascqgh9aqpcqghsisldaaewaqaaiabdchawnaabaaaleaacqwgqbgaaeqaaaqaaiabdqgaqjabg2da9iabigdaxaqaaiabikdayiabicdawaqdcqghris5aogagiibawmaei4ba8maei4zacmaemicaa3aasbaasqaaiabdqgaqbqabaaaaa@40ee@ is the shannon entropy, and j =..., also represents the index of the twenty amino acids. for a set of sequences such as all coding and noncoding orfs from a given genome, the vectors {si} scatter widely in the 20-dimensional edp phase space. our basic working hypothesis is that the edp vectors for coding orfs form separate clusters from those for noncoding orfs, for a possible reason that they moved along separate paths during the evolution due to different selection pressure. if this hypothesis is true, one can then develop a clustering method in the edp phase space that discriminates the coding orfs from the noncoding ones. the success of the method confirms the validity of this hypothesis. our previous study shows that there exists indeed one universal center representative of coding orfs and one universal center of noncoding orfs in the edp phase space for dozens of analysed species. for gc-rich genomes, our current research based on principal component analysis shows that orfs form six clusters in the edp phase space due to the genomic gc content bias, one for coding orfs while other five for non-coding orfs, which is also consistent with a previous study. thus the noncoding orfs of gc-rich genomes are represented with five universal centers in the edp phase space. following the zcurve algorithm, we refer to the gc-rich genome as a genome with gc content higher than 56%. consequently, the universal centers allow us to construct an iterative procedure to define a refined coding potential. this procedure begins by finding a set of "root" coding and noncoding orfs coding} and {snoncoding}), which have been demonstrated with high reliability based on the universal edp centers, and then builds a coding potential of a dna sequence s like an orf as the ratio of the distance to each of the edp clusters Φedp=min⁡kdcoding)min⁡ldnoncoding),      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaahaawcbeqaaiabdweafjabdseaejabdcfaqbaakiabcicaoiabdofatjabcmcapiabg2da9maalaaabagagiyba0maeiyaakmaeioba42aasbaasqaaiabdugarbqabagccqwgebardaqadaqaaiabdofatjabcycasiabdofatjabcicaoiabdugarjabcmcapmaacaaaleqabagaem4yammaem4ba8maemizaqmaemyaakmaemoba4maem4zacgaaagccagloagaayzkaaaabagagiyba0maeiyaakmaeioba42aasbaasqaaiabdygasbqabagccqwgebardaqadaqaaiabdofatjabcycasiabdofatjabcicaoiabdygasjabcmcapmaacaaaleqabagaemoba4maem4ba8maemoba4maem4yammaem4ba8maemizaqmaemyaakmaemoba4maem4zacgaaagccagloagaayzkaaaaaiabcycasiaaxmaacawljawaaewaaeaacqaiyagmaiaawicacaglpaaaaaa@6c5f@ where k and l are the indices of two clusters of root coding and noncoding orfs, d can be defined by the euclidean distance between two edp vectors s = {si} and s' = {s′i mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgzbwcgaqbamaabaaaleaacqwgpbqaaeqaaaaa@2fae@} as d = 2)1/ i =..., a noteworthy feature of this coding potential is that it is able to generate a cluster with non-trivial boundary and thus achieve a classification of the two classes of sequences with remarkable accuracy. for this reason, the edp model has been successfully used to produce training orf set in the newest version of glimmer method. the tis model for translation start site of prokaryotic genes the tis model in med, which is a further extension of the rbs model described earlier, makes a conscious effort to model the following prominent features for prokaryotic gene structure: upstream to the tis there exists one or a few common motifs that act as the binding sites for the initiation of translation or transcription, these common motifs usually have a position-specific distribution in upstream region from the start sites; among various genomes, the usage of start codon atg, gtg and ttg is different; for orfs with multiple candidate start codons, the leftmost start codons are usually the tiss. however, those start codons located after the leftmost ones are also chosen to be tiss according to different probability; the distance between the neighbor genes has a characteristic probability distribution; for gc-rich genomes, there exists a high gc content at the codon first and third position due to the gc content bias. these features are included in the tis model, allowing us to iteratively find the true tis after the determination of the right 3' end coding orfs. it is notable that no training data is required to fix the model parameters. the various probabilities mentioned above are obtained by non-supervised learning during the iteration. it is intriguing to note also that these probabilities are genome-specific quantities and their differentiation between genomes may likely be interesting and significant for comparative genomic studies. for the rbs model in our previous work, two likelihood functions are defined, one being associated with the shine-dalgarno motif, Φmotif, and another being directly related to the start codon, Φstart. the total likelihood reads Φtotal = Φmotif + Φstart. we have made several crucial extensions of that model as follows. first, many archaea are detected to have both sd signals of bacterial type in near upstream region and box a motifs of eukaryotic type in the further upstream region from tis. so, med defines two motif searching regions, vup = and vlow =, and carries out search for candidate motifs separately in each region. secondly, a new likelihood function is defined to describe the probability that two orfs overlap, which is given by a distribution of the distance between a start codon and the stop codon of the immediate upstream orf. thirdly, for gc-rich genomes we use a parameter cp defined by formula to describe the difference between the coding region vc = downstream and the noncoding region vnc = upstream to tis. it has long been known that the sequence patterns of genomes with gc-rich content show differences to those with normal or low gc content. as a result, high gc content at the codon third and first position has been reported for gc-rich genomes. similar to nishi et al., we refer to gc as the g+c occurrence at the ith codon position, where i = or then parameter θ=gc+gcgc+gc+gc mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf4oqccqgh9aqpdawcaaqaaiabdeeahjabdoeadjabcicaoiabigdaxiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabiodaziabcmcapaqaaiabdeeahjabdoeadjabcicaoiabigdaxiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabikdayiabcmcapiabgucariabdeeahjabdoeadjabcicaoiabiodaziabcmcapaaaaaa@4a17@ is used to describe the nucleotide usage at the codon first and third position. for both two regions vc and vnc, we treat them as a set of codon sequences regardless of the fact that the vnc are noncoding sequences. our analysis shows that the probability distribution of usage pc for coding region vc is significantly different from that of usage pnc of noncoding region vnc. thus for a dna sequence with a usage θ, the cp parameter is designed to describe the probability of the sequence belonging to coding region as follow cp = pc/ + pnc),     where ≤ cp ≤ to evaluate the likelihood of a start codon as tis, we calculate the nucleotide usage θc and θnc for its downstream region and upstream region separately. thus the likelihood function is defined by the formula Φgcstart mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ = * cp,     where ≤ Φgcstart mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ ≤ note that the score Φgcstart mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ is calculated only for the genomes with gc content higher than 56%, and this item is added into the start codon likelihood function Φstart developed in our previous work. fourthly, med introduces a measure of signal strength in order to counter-balance the genomic nucleotide content bias. let t be the ensemble of all l-mers in the set v of upstream sequences of all candidate coding orfs, and for each l-mer t ∈ t, denote by p^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqcaaaa@2e25@ the observed frequency in the region v. denote also by p¯ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ the frequency of the same l-mer if it is a random combination of its nucleotides; thus, p¯ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ can be estimated from single nucleotide composition measured from the whole genome. similarly to fuglsang, we refer to x = - p¯ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@)2/p¯ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgwbacgaqeaaaa@2e2d@ as a measure of the over-representation of t ∈ t. when summing over the ensemble t, we have a measure, χ for deciding whether tupor tlow is more likely to contain a strong signal in vup or vlow, where tup means the ensemble of l-mers in the region vup. such measure is also used to refer to a weight for each region, as mup = χ2/ + χ2), and mlow = χ2/ + χ2) which allows us to obtain a combined measure for a given candidate start codon Φtis = mupΦup + mlowΦlow + mstartΦstart,     here, mstart = max. in above formula, score Φup reflects the likelihood of a motif associated with transcription initiation signal in region vup =, while Φlow means the likelihood of a motif associated with translation initiation signal in region vlow =, both two scores include a weight matrix scoring function and an occurrence probability in upstream region to tis. for score Φstart, it includes weight matrix scoring function of sequences around tis, the probability of start codon as tis compared with the leftmost start codon, and the orf overlap scoring function. if the genomic gc content is higher than 56%, the scoring function given by formula will be added in. note that the new likelihood function Φtis is able to describe both bacterial and archaeal genomes in a unified way. finally, a distinction is made between transcript unit internal genes and transcript unit leader genes. following tolstrup et al. and torarinsson et al., we refer to tui genes as those with a start codon separated from the stop codon of the nearest upstream gene by < bps, and tul genes otherwise. our analysis in s. solfataricus shows that only sd-like motifs are found for tui genes while box a motifs are usually found for tul genes. in med, each candidate start codon then definitely belongs to one of the two above categories, each of which conducts the search for motifs separately. furthermore, med forms a measure λ=−χ2)−χ2),      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf7oabcqgh9aqpcqggoaakcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgjbqsaeaacqwg1bqdcqwgwbacaagccqggpaqkcqghsislcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgjbqsaeaacqwgsbabcqwgvbwbcqwg3bwdaagccqggpaqkcqggpaqkcqggoaakcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgmbataeaacqwg1bqdcqwgwbacaagccqggpaqkcqghsislcqwfhpwydaahaawcbeqaaiabikdayaaakiabcicaoiabdsfaunaadaaaleaacqwgubavcqwgvbqvcqwgmbataeaacqwgsbabcqwgvbwbcqwg3bwdaagccqggpaqkcqggpaqkcqggsaalcawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@6d64@ to decide which of the two translation initiation mechanisms is predominant in a genome: λ > for either completely sd signals or completely box a consensus, and λ ≤ for a mixture of the two. tests show that most of bacteria and part of archaea belong to the former, whereas many of archaea belong to the latter. in the latter case, med uses the z-score to form a normalized tis score which allows to compare and select the right tis. therefore, our method allows to describe more than one type of initiation signals which are located at different upstream positions. the iterative learning algorithm the gene finding system med predicts genes in two stages: the coding orf detection stage and the tis refinement stage. the outline of sequence processing and model self-training for two stages are shown in fig. we give a summary of this algorithm as follows. in the first step of coding orf detection stage, all orfs longer than bps are extracted in both strands and in all three reading frames, thus to determine the edp coding potential. as stated above, using the universal coding and noncoding edp centers and with an iteration of consecutive root orf finding, we decide a set of root coding/noncoding orfs for each genome being analysed, which typically cover more than 60% of all genes with a very high reliability over 99%. note that we obtain both "seed" coding orfs and non-coding sequences, instead of generating artificial non-coding orfs as other methods do. these seed sequences form a reliable learning set for the further analysis. the system then goes to determine the remaining orfs excluded from the root orfs. we search for conservative motifs only in the near upstream region vlow of the root coding orfs, namely setting mup = and mlow = mstart = in calculating the tis score in formula. the candidate start site with the highest tis score is selected, then the parameters are recalculated using the updated start sites. this iteration process ends when no start site needs to be relocated. the parameters are then obtained and used to determine the remaining orfs. this consists in calculating the tis score of the remaining orfs, and then evaluating their coding potential using a fisher discriminant algorithm in a two-dimensional plane of the edp coding score and the tis score. for gc-rich genomes, the fisher discriminant algorithm is applied in a three-dimensional space of the edp coding score, the tis score and the θ score, which is trained using p. aeruginosa. our studies show that the discriminant coefficients are universal across prokaryotic genomes. thus the remaining orfs are classified into coding orfs and noncoding ones. to reduce the false positive in the detected coding orfs determined from the remaining orfs, a strategy to resolve overlap of coding orfs is then applied. we define a quantity ξ to be the percentage of nucleotides overlapped with other coding orfs in either the same or the opposite strand. the iteration goes as follows: extend all the detected coding orfs to the longest orfs by assigning the leftmost start codon as the start of each orf; remove those orfs with ξ > ξi; recalculate the overlapping percentage for all orfs retained from the last step; and reduce ξi by 1%, then repeats from step. this procedure iterates until ξ reaches a threshold. after the stage of coding orf detection, med goes to the next stage of tis refinement for all coding orfs. the detected coding orfs are first classified into two sets, tui and tul. the system considers searching for either sd motifs or box a motifs in both regions vup = and vlow =. then, one forms the combined tis score in formula, which is calculated together with all three coefficients mup, mlow and mstart described above. for gc-rich genomes, the score Φstart includes the item of Φgcstart mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqhmogrdaqhaawcbagaem4rackaem4qameabagaem4camnaemidaqnaemyyaemaemocainaemidaqhaaaaa@3782@ in formula. then the global genomic measure λ is used to decide which of the initiation mechanisms is more appropriate for the genome studied. if λ > the self-training of tis model parameters runs in both vup and vlow for all detected coding orfs, we select the sd-like motif as the main mechanism if χ <χ or the box a motif otherwise. for each coding orf, the tis score in formula is calculated for all candidate start codons, then the start codon with the highest score is selected as the predicted tis. if λ ≤ the self-training is run with vup and vlow together, but for tui and tul genes independently. for each coding orf, the normalized scores of both tui type and tul type are calculated, and then compared to select one with higher score as the most likely start codon to relocate the gene start.  genome sequences and reliable gene datasets the bacterial and archaeal genomes and their annotations used in this paper were downloaded from the genbank release in we have selected experimentally confirmed genes as benchmarks. for e. coli, we have used two datasets, the ecogene and the link dataset, with their newest version. for b. subtilis, following the standard by besemer et al., we have taken it as one of a few annotated complete genomes that can be used to evaluate the performance of exact gene prediction. thus we denoted all genes annotated in genbank as bsub_all dataset. for verifying short genes, we have chosen bsub bsub and bsub datasets, which are selected from the b. subtilis and verified by protein similarity search. we have extracted short genes from the genes of ecogene, denoted as ecogene_short, to evaluate short gene prediction for e. coli. for gc-rich genomes, we have built two data sets, mtub and paer the former includes reliable genes with confirmed tiss in m. tuberculosis with gc content of %, while the latter has such genes in p. aeruginosa with gc content of %. for archaeal genomes, we have built a data set solfgene including reliable genes in s. solfataricus with confirmed tiss. the three data sets mtub paer and solfgene, as well as the process of determining these genes with their tiss confirmed by n-terminal protein sequencing or inferred from experimental evidences, can be accessed through. other benchmark programs in order to benchmark med, we have tested four gene-finding methods, glimmer, genemarks, zcurve and easygene, to make a comparison of the performance of gene prediction. the former three methods belong to the same category of our system as the ab initio gene prediction methods, while the easygene belongs to another category of methods using the extrinsic information by means of similarity search. glimmer was downloaded from and installed locally. we ran it following the instructions given in the distribution file. a post-processor rbsfinder has been designed to further improve the tis prediction by glimmer. thus, rbsfinder was used to process the original output of glimmer and the refined tiss are then taken as predicted ones by glimmer system. rbsfinder is accessible via the website at. herein, rbsfinder was run repeatedly until over 99% of gene starts remain unchanged. recently, the newest version of glimmer, i.e. v, has been released with several algorithmic changes to reduce the number of false positive predictions and to improve the accuracy of tis predictions. although the article describing glimmer has not yet been found from publication, a locally executable program is current available from. in order to have a comprehensive comparison with glimmer method, glimmer was further included herein as benchmark program. we downloaded the predictions by glimmer from its website for all genomes studied in this paper. the software zcurve is freely available at. it was executed on our pc following the instruction file. genemarks provides only an online service of genomic annotation, instead of a local executable program. the newest results on several test sets mentioned above for e. coli and b. subtilis by the current genemarks have been greatly improved than that of its initial version. similarly, there is still no new literature to report the improved version of genemarks. however we adopted all the unpublished results returned to us via email with the current online analysis behind the interface, to compare with our program. similar to genemarks, easygene is only accessible via its web interface, and the pre-trained models are available for only genomes. however, the easygene web server has provided predictions for hundreds of chromosomes. therefore, for comparison purpose, easygene predictions for a total of genomes were downloaded from its website available at the time we prepared our paper. accuracy of gene detection and interpretation of genome-specific model parameters to illustrate the prediction accuracy of med, we first present the comparison against the genbank annotation for the 3' end match. two independent quantities, sn and sp, are defined to evaluate the performance of a gene finder at gene level as: sn = tp/, sp = tp/.     here, tp, fp, and fn are the number of true positive, false positive, and false negative, respectively. the med program has been run on all the complete bacterial and archaeal genomes currently available on genbank. additional file includes the prediction results on forty representative genomes, comprising bacteria and archaea. selected organisms cover most of the taxonomic groups, including crenarchaeota, euryarchaeota, nanoarchaeota for archaea and α/β/γ/ε-proteobacteria, actinobacteria, bacteroidetes/chlorobi, for bacteria. for the forty species, med achieves an average sensitivity of % and specificity of %. the results of glimmer, glimmer, zcurve and the newest genemarks for the forty genomes as a comparison are not shown in additional file but may be found in our website. on the average, the newest version of genemarks gives the same level of sensitivity and the highest specificity, while glimmer and zcurve give the sensitivity of % and %, and the specificity of % and %, respectively. as reported by easygene's website, the average sensitivity and specificity are % and % respectively for of the genomes. compared with glimmer on average, glimmer significantly reduces the false positives with almost no decreasing of the sensitivity of %. the average sensitivity and specificity of genemarks reported in the original version, which was run on eight genomes including e. coli, b. subtilis and a. fulgidus, are % and %, respectively. a full comparison for genomes has also been listed at our web page. for the genomes, med gives the average sensitivity of % and specificity of %, while glimmer gives % and %, zcurve gives % and %, glimmer gives % and %, genemarks gives % and %, and finally easygene gives % and % averaged over genomes, respectively. however, we argue that the higher specificity of the latter three programs is not unrelated to the fact that they are highly optimized with the genbank annotation. it is clear that med is competitive with glimmer and zcurve, while the specificity of prediction is a little higher than both of them. both glimmer and the newest genemarks present rather high accuracies of sensitivity and specificity against the current genbank annotation for the 3' end match, although there are not yet publications to report them. using the extrinsic information of similarity search, easygene is able to extremely raise the specificity at the cost of its sensitivity slightly lower than the ab initio gene prediction methods. it should be noted that easygene has included many a confirmed gene of the query-genome as training set by using blastp to search for significant protein matches in swiss-prot, which makes an essential distinction between such a method and the ab initio gene prediction method such as med. for instance, there are about 65% of all genes in e. coli predicted by easygene have been taken as the training set. therefore a comparison against all genes in the benchmark would over-estimate easygene's prediction performance. to discuss the predicting performance on the current genbank annotation, we should point out that it can not be free from the bias since both glimmer and genemark series have been widely used or involved in the genbank annotation pipeline. the statistical analysis shows that there are over genomes from genbank before are in this case. while the comparison against the genbank annotation has played a role during the early development of the computational prediction programs, this comparison becomes increasingly suspicious as its bias becomes evident. since the genbank annotation is not fully accurate, further evaluation is performed based on the function-known genes which have more reliable annotations. for each genome, the function-known genes are selected from genbank by excluding those with product descriptions with any of the key words as "-like", "conserved", "hypothetical", "homolog", "probable", "possible", "predicted", "putative", "similarity" and "unknown". the ratio of the number of these genes predicted correctly by med is also listed in additional file note that the specificity is meaningless in this case. as we can see, the average accuracy for med is %, while the average accuracy for glimmer, glimmer, zcurve, easygene and the newest genemarks is %, %, %, % and %, respectively. for genomes analysed, the average accuracy for med, glimmer, glimmer, zcurve, easygene and the newest genemarks may be calculated, which gives %, %, %, %, % and % respectively, clearly, they are all at the same level. these results can be found on our website. therefore the gene-finding accuracy for function-known genes of med and other five gene finders is well matched. a distinctive feature of med is that the algorithm builds a prokaryotic gene structure model based on a set of genome-specific parameters, which are calculated by several self-learning iterations without any prior knowledge or training data. comparing the model parameters among different species, e.g. those of gene-related properties as tis-upstream signals and usage of start codons atg, gtg and ttg, would shed light on the study of genomic comparison and evolution for bacteria and archaea. to this end, we present in additional file the parameters of start codon usage and motifs found upstream to tis for each species, although a more detailed list is provided in the output file by our program for the user. genomes in additional file are listed as two groups of bacteria and archaea and sorted alphabetically. almost for all analysed species, atg is the most common start codon compared with ttg and gtg. an exception is the m. kandleri genome, for which the three start codons are used nearly equally. in contrast, the usage of ttg to gtg is different for different species. for example, gtg in c. glutamicum is much more favored than ttg; while ttg in b. burgdorferi is more favored than gtg. for archaeal genomes, the start codon usage predicted by med is very close to an earlier result annotated by easygene including the protein similarity search method. with the automatic control coefficient λ in formula and the resulting predicted motifs obtained by a non-supervised learning process, med clearly defines three types of translation initiation mechanisms and ascertains which of the types each genome belongs to. in additional file we go further and report the predicted motifs, meaning the most significant consensuses associated with translation initiation signals in two searching regions and upstream to tiss. for all of analysed bacterial genomes, motifs are presented only in region, while most of them belong to sd signal type, which means the typical mechanism of translation initiation with the sd signals for transcripts with leaders. however, at-rich motifs, "aattt" and "aaatt" listed in additional file are detected in region for the synechocystis sp. genome, indicating the existence of a leadless translation initiation mechanism in cyanobacteria. a similar case is to the b. thetaiotaomicron genome, where a-rich motifs are detected. for analysed archaeal genomes, three types of mechanisms are detected. six species m. jannaschii, m. kandleri, m. maripaludis, m. thermoautotro., p. abyssi and p. furiosus are reported with motifs of sd signal type only in region, meaning the mechanism similar to the bacterial. in contrast, three genomes halobacterium sp., n. equitans and p. aerophilum are shown that the motifs of box a type are predicted only in region, which suggest a different mechanism of translation initiation with a high level of genes producing leaderless transcripts. for three genomes a. fulgidus, p. torridus and s. solfataricus, both box a motifs in region for leaderless transcripts and sd motifs in for transcripts with leaders are detected. it implies that a mixed mechanism of leaderless translation and leadered translation is used in these genomes. the result is well matched with that of previous studies. a notable result is that most of the motifs of searching region found in bacteria contain the tetramer "ggag" in contrast with that of "ggtg" in most archaea, it is also exactly consistent with that reported in the earlier study. however, note that for the thermotoga maritima genome, the med analysis leads to an interesting finding of sd-like motifs very richly observed in archaea. a phylogenetic tree based on 16s rrna sequence has revealed that t. maritima is a deep branching species. moreover, it has been reported that over 20% of genes in t. maritima were acquired through horizontal gene transfer from archaea. thus, the archaea-like feature in t. maritima might be interpreted as either an ancient feature or, alternatively, a consequence of horizontal gene transfer. accuracy of exact prediction on reliable genes in order to perform an accuracy test on the exact gene prediction, reliable gene datasets must be used since the genbank annotations have systematic bias in the 5' end prediction. before showing the results in this section, we should like to stress that all six gene finders, including med, are compared to analyze each genome with the same sensitivity and specificity against the genbank annotation as mentioned above. as a gene finder using a similarity search by alignment in swiss-prot to build tls-confirmed orfs, many a protein coding gene has been included as training set of tis for easygene. for example, there are of the genes in bsub_all dataset, and of the genes in ecogene, have been employed as its training set for b. subtilis and e. coli respectively. thus for a fair comparison with the ab initio methods, the accuracies of easygene on reliable test sets were calculated by excluding the number of genes included into its training set. first we present the prediction accuracy on reliable test sets bsub_all, ecogene and link for bacterial genomes. as indicated in table for only 3' end match, all six programs have nearly the same performance with accuracy above %. for both 5' and 3' end matches, the current version of genemarks achieves the highest performance with average accuracy of %, while med has a comparative performance with average accuracy of %. note that for bsub_all dataset, genemarks has reported its accuracy of % for 3' end match and % for both end matches on genes annotated in genbank, which is lower than the newest result by the current genemarks version, as well as lower than that of med. zcurve shows slightly lower average accuracy of % for both end matches, and easygene does %, while glimmer gives the lowest average accuracy of %. glimmer has significantly improved the average accuracy to %, reaching a similar performance as med. one of the remaining challenges for prokaryotic gene finding is the identification of short genes'buried' in an enormous pile of false-positive short orfs. following the standard adopted by many gene finders, we hereby define short genes as those with length between and bps. we have carried out a test on short genes from confirmed ones in ecogene data set. the test indicates that med has an accuracy of % for 3' end match and also a high accuracy of % for both end matches. for the bsub bsub and bsub datasets from b. subtilis, med detects %, % and % of 3' end match and %, % and % of both end matches, respectively. similarly, note that the corresponding published results of genemarks on the same three datasets are %, % and % of 3' end match as well as %, % and % of both end matches, which have been greatly improved given by the newest online version. however as listed in table comparing with other five gene finders shows that med gives a higher performance in the short gene prediction than glimmer, glimmer, easygene and zcurve, while is comparable with genemarks. aprograms med, gl, gl and zcv were run locally, while gmk was run online, as described in the text. predictions for eg were downloaded from. bexperiment confirmed tiss data sets: the first three represent two well-studied genomes: b. subtilis and e. coli; the fourth to seventh represent short genes for e. coli and b. subtilis; psaer and mtub are selected for two gc rich genomes, m. tuberculosis and p. aeruginosa; solfgene corresponds to the archaeal s. solfataricus. cnumbers in parentheses indicate that the results of genemarks have been reported in literature, means no data reported. to evaluate the performance of exact gene prediction for gc-rich genomes, two reliable sets paer and mtub are used here for p. aeruginosa with gc content of % and m. tuberculosis with gc content of %. test on paer shows that med has an accuracy of % for 3' end match and % for both end matches, while on mtub gives an accuracy of % for 3' end match and % for both end matches. as we can see from the eighth and ninth lines in table compared with glimmer, glimmer, zcurve, easygene and the newest online version of genemarks, med demonstrates the highest performance of exact gene prediction for gc-rich genomes. it is still problematic to evaluate gene prediction for archaeal genomes due to the insufficient number of genes confirmed by independent methods. to this end, we have built the data set solfgene comprising reliable genes with confirmed tiss. all six programs detect all genes at the 3' end in solfgene. as for both ends prediction, both med and easygene give the highest accuracy of % vs glimmer, glimmer, genemarks, and zcurve. in summary, upon test of the reliable gene sets built for two well studied genomes e. coli and b. subtilis, we show that the total prediction performance of med can be matched with that of genemarks in its publication, and higher than that of glimmer, easygene and zcurve. while for the newly analysed genomes, especially for gc-rich genomes and archaea, med shows an evident advantage, because test on the reliable gene sets clearly indicates that med has given a higher accuracy than other five gene finders, at least the same level of the best one. analysis of upstream region of predicted tiss for archaea med gives a comprehensive modeling of the tis for both bacterial and archaeal genomes, in particular for archeal genomes in which translation initiation is more complex than in bacteria. we report herein the results using the extended sequence logos for the region upstream from predicted tiss by med. the sequence logos facilitate detection of sequence patterns that are conserved in both content and position. we present the logos of three representative archaeal genomes: m. jannaschii, p. abyssi, and n. equitans in fig. the logos indicate that several mechanisms associated with the translation initiation have been detected by med. the genome of m. jannaschii displays the logos with a remarkable sd signals in region like bacteria, which is well compatible with the notion that this genome produces nearly no leaderless transcripts. on the other hand, the n. equitans genome that carries a high level of genes producing leaderless transcripts produces a very rich set of signals. an evident a/t-rich signal, the so-called box a promoter motif, is observed near the position - in addition, a tfb recognition element motif constituting to a/ts, which interact with the archaeal transcription factor tfb, is also recognized about bps upstream from the box a motif. furthermore, an a/t peak is detected near the position - to - upstream to the tis, which has been reported being functionally important. finally, the p. abyssi genome is found to have a strongly conserved sd signal in - to - upstream region together with a rather weak box a motif upstream from the sd signal. this is well consistent with the fact that the p. abyssi genome produces fewer leaderless than leadered transcripts. inspection on the upstream sequence patterns of predicted tiss would be helpful for the assessment of the overall quality of gene annotation or prediction. herein we choose an example, the s. solfataricus genome, for a detailed comparison between med and other prediction programs, as well as the genbank annotation. first we contrast med with genbank annotation. as we can calculate from additional file med detects genes with common 3' end against genbank annotation, among which genes have common tiss and genes do not. fig. 3b shows the logos of the upstream region from tiss of the common genes, where a typical sd signal can be observed at and an evident box a motif at. the logos for the other tiss predicted by med are shown in fig. 3c, where similar sd signal and box a motif are found, while those tiss from the genbank annotation shown in fig. 3d show no informational structure in the logos. this comparison suggests that med can give a more relevant annotation of tis for archaeal genomes than the current genbank annotation. similar comparison is made with glimmer, genemarks and zcurve, as shown in fig. we have shown that glimmer has achieved an overall improvement of tis prediction than that of its earlier version, thus have no need for the comparison with glimmer. the sets of common tis predictions are slightly different in each case, but in all comparisons the tiss predicted by med show the sequence logos similar to that in fig. 3b, while those of glimmer, genemarks and of zcurve show almost no sd signal or box a motif structures. thus to a great extent, med makes a substantially more accurate prediction than glimmer, genemarks and zcurve for archaeal genomes.  in this paper, we present a comprehensive model to describe a set of properties about the coding potential and the translation initiation mechanisms for both bacteria and archaea. based on the model, with a design of multiple iterations, an non-supervised ab initio gene prediction system med is developed. generally speaking, the system is able to adapt to any newly sequenced prokaryotic genome with no need for any data training or prior knowledge, and can predict the divergent translation initiation mechanisms and the resulting signals upstream from the tis. the model seems to be biologically sound, since the program yields genome-specific model parameters such as various probabilities associated with the translation initiation signals and the start codon usage, which are matched with the current knowledge from earlier works. thus the model may provide a good tool for comparative genomic studies. upon test of a set of reliable gene sets, the total prediction performance of med for the well studied bacteria such as e. coli and b. subtilis with usual gc content can be matched with that of the existing published methods. while for the newly analysed genomes, especially for gc-rich genomes and archaea, test on the reliable gene sets indicates that med outperforms or at least gives the same level of the best of the current gene finders in exact gene prediction for both 3' and 5' end matches. furthermore, med adapts to a broad range of archaeal genome as well as to bacterial genomes. for archaeal genomes with more complex mechanisms of translational initiation, our method has a more accurate prediction of tiss compared to the existing gene finders and the current genbank annotation. as many a prokaryotic gene prediction program, such as glimmer, genemark and easygene, has been widely used in raw genomic sequence annotation, the predicting performance for a new gene prediction method evaluated by data extracted from genbank should be in caution. while the comparison against the genbank annotation has played a role during the early development of the computational prediction programs, this comparison becomes increasingly suspicious as its bias becomes evident. therefore the consistency with the genbank annotation should not be considered to be the only factor, and should not even be the essential one for some archaeal genomes. with hundreds of prokaryotic genomes have been sequenced and made publicly available, subsequent studies will focus more on the development of the comprehensive model of gene structure that allow capture of the biological evidence of prokaryotic genomes, e.g. the translation initiation mechanism, which is essential to understanding comparative genomics and biology evolution. at present the problems of being biased and erroneous have not yet been solved in the genbank annotation for prokaryotes, especially for newly sequenced archaeal genomes and gc-rich genomes, and the further promotion of comparative genomics, it is hoped that our system is shown to meet the demand and be qualified as an alternative tool of gene computational prediction. the expert annotators, who are concerned with more accurate and complete annotation of sequenced genomes, should find our method useful when used independently or alongside with other gene finding tools. availability and requirements of the system project name: med project; project homepage:; operating systems: microsoft windows operating system is recommended for the pre-complied program downloaded from our web site. in addition, we also provide the source code on linux with a c++ compiler; programming language: c++; licence: the sourc code is freely available under the gnu gpl license; any restrictions to use by non-academics: none. authors' contributions zss and hqz co-supervised the whole development of the work and had the main idea of the algorithm, gqh and hqz designed and implemented the gene-finder system, drafted most parts of the manuscript. yfy and jw contributed biological expertise and performed part of the experimental test. all authors read and approved the final manuscript. supplementary material additional file 1 med prediction accuracy and genomic-specific model parameters. contains data for med prediction accuracy and genomic-specific model parameters for bacterial genomes and archaeal genomes. the file is edited by microsoft excel click here for file additional file 2 source code. contains source code of med written for windows and linux/unix operation systems. click here for file acknowledgements we thank zhengqing ouyang, yongchu liu, zongxiao sun, xiaobin zheng and others for beneficial discussions and helps to the work. we thank also j. o'brien for kindly reading and revising the manuscript. we also thank prof. chun-ting zhang for providing zcurve program. the work received partial support by the national natural science foundation of china; the present study was also supported by the national basic research program of china under grant no. 2003cb 
results of extensive tests show that med achieves a competitive high performance in the gene prediction for both 5' and 3' end matches, compared to the current best prokaryotic gene finders. in particular, med is shown to reveal divergent translation initiation mechanisms in archaeal genomes while making a more accurate prediction of tiss compared to the existing gene finders and the current genbank annotation. this paper describes a new prokaryotic genefinding algorithm based on a comprehensive statistical model of protein coding open reading frames and translation initiation sites. the iterations enable us to develop a non-supervised learning process and to obtain a set of genome-specific parameters for the gene structure, before making the prediction of genes. furthermore, the genome-specific parameters given by med match with the current understanding of prokaryotic genomes and may serve as tools for comparative genomic studies. they are combined to form a so-called multivariate entropy distance algorithm, med, that incorporates several strategies in the iterative program. despite a remarkable success in the computational prediction of genes in bacteria and archaea, a lack of comprehensive understanding of prokaryotic gene structures prevents from further elucidation of differences among genomes. the former is based on a linguistic "entropy density profile" model of coding dna sequence and the latter comprises several relevant features related to the translation initiation. it continues to be interesting to develop new ab initio algorithms which not only accurately predict genes, but also facilitate comparative studies of prokaryotic genomes. the advantage of the med is particularly evident for gc-rich genomes and archaeal genomes.
 despite a remarkable success in the computational prediction of genes in bacteria and archaea, a lack of comprehensive understanding of prokaryotic gene structures prevents from further elucidation of differences among genomes. it continues to be interesting to develop new ab initio algorithms which not only accurately predict genes, but also facilitate comparative studies of prokaryotic genomes.  this paper describes a new prokaryotic genefinding algorithm based on a comprehensive statistical model of protein coding open reading frames and translation initiation sites. the former is based on a linguistic "entropy density profile" model of coding dna sequence and the latter comprises several relevant features related to the translation initiation. they are combined to form a so-called multivariate entropy distance algorithm, med, that incorporates several strategies in the iterative program. the iterations enable us to develop a non-supervised learning process and to obtain a set of genome-specific parameters for the gene structure, before making the prediction of genes.  results of extensive tests show that med achieves a competitive high performance in the gene prediction for both 5' and 3' end matches, compared to the current best prokaryotic gene finders. the advantage of the med is particularly evident for gc-rich genomes and archaeal genomes. furthermore, the genome-specific parameters given by med match with the current understanding of prokaryotic genomes and may serve as tools for comparative genomic studies. in particular, med is shown to reveal divergent translation initiation mechanisms in archaeal genomes while making a more accurate prediction of tiss compared to the existing gene finders and the current genbank annotation. 