Combining diverse evidence for gene recognition in completely sequenced bacterial genomes [published erratum appears in Nucleic Acids Res 1998 Aug 15;26(16):following 3870]

doi:10.1093/nar/26.12.2941

This Article

	Abstract
	Print PDF (144K)
	Alert me when this article is cited
	Alert me if a correction is posted

Services

	Email this article to a friend
	Similar articles in this journal
	Similar articles in ISI Web of Science
	Similar articles in PubMed
	Alert me to new issues of the journal
	Add to My Personal Archive
	Download to citation manager
	Search for citing articles in: ISI Web of Science (104)
	Request Permissions
	Commercial Re-use Guidelines for Open Access NAR Content

Google Scholar

	Articles by Frishman, D.
	Articles by Gelfand, M.
	Search for Related Content

PubMed

	PubMed Citation
	Articles by Frishman, D.
	Articles by Gelfand, M.

Social Bookmarking

	What's this?

Nucleic Acids Research Pages 2941-2947

Combining diverse evidence for gene recognition in completely sequenced bacterial genomes
Introduction
Materials And Methods
   Data
   Outline of the algorithm
   Similarity search and the set of seed ORFs
   Coding potential and the complete set of candidate ORFs
   RBS weight matrix and assignment of gene starts
   Choice of the minimal allowed ORF length and handling of short ORFs
   Implementation and availability
Results And Discussion
Acknowledgements
References

Combining diverse evidence for gene recognition in completely sequenced bacterial genomes

Dmitrij Frishman*, Andrey Mironov¹, Hans-Werner Mewes, Mikhail Gelfand²

Munich Information Center for Protein Sequences (MIPS) of the German National Center for Health and Environment (GSF), Am Klopferspitz 18a, 82152 Martinsried, Germany, ¹Laboratory of Mathematical Methods, National Centerfor Biotechnology Information NIIGENETIKA, Moscow 113545, Russia and ²Institute of Protein Research,Russian Academy of Sciences, Pushchino 142292, Russia

Recieved March 3, 1998; Revised and Accepted May 4, 1998

ABSTRACT

Analysis of a newly sequenced bacterial genome starts with identification of protein-coding genes. Functional assignment of proteins requires the exact knowledge of protein N-termini. We present a new program ORPHEUS that identifies candidate genes and accurately predicts gene starts. The analysis starts with a database similarity search and identification of reliable gene fragments. The latter are used to derive statistical characteristics of protein-coding regions and ribosome-binding sites and to predict the complete set of genes in the analyzed genome. In a test on Bacillus subtilis and Escherichia coli genomes, the program correctly identified 93.3% (resp. 96.3%) of experimentally annotated genes longer than 100 codons described in the PIR-International database, and for these genes 96.3% (83.9%) of starts were predicted exactly. Furthermore, 98.9% (99.1%) of genes longer than 100 codons annotated in GenBank were found, and 92.9% (75.7%) of predicted starts coincided with the feature table description. Finally, for the complete gene complements of B.subtilis and E.coli, including genes shorter than 100 codons, gene prediction accuracy was 88.9 and 87.1%, respectively, with 94.2 and 76.7% starts coinciding with the existing annotation.

INTRODUCTION

The principal goal of large-scale genome sequencing is to obtain new insights into physiological and biochemical processes in living organisms. An essential step in this process is gene identification with subsequent computer-based annotation of the corresponding gene products. Although bacterial genomic sequences are devoid of introns, gene recognition in bacteria is far from being simple. It is easy to extract all possible open reading frames (ORFs) from a given DNA sequence; it is much less trivial to decide which of them correspond to genes that are actually expressed and code for proteins. The following features are important indicators of protein coding regions in DNA: (i) sufficient ORF length. Long ORFs rarely occur by chance; (ii) specific patterns of codon usage that are different from triplet frequencies in non-coding regions (`coding potential'); (iii) the presence of ribosome binding sites (RBS) in the (-20)...(-1) region upstream of the start codon that help to direct ribosomes to the correct translation start positions (1). A part of the RBS is formed by the purine-rich Shine-Dalgarno (SD) sequence which is complementary to the 3[prime] end of the 16S rRNA (2); (iv) similarity to known, especially experimentally characterized, gene products.

Correspondingly, the approaches to gene recognition are traditionally divided into two broad categories (3,4). Intrinsic, or ab initio methods utilize statistic, linguistic or pattern recognition algorithms to find genes in DNA through detection of specific nucleic acid motifs or global statistical patterns, whereas extrinsic methods take into account information about other known proteins.

There exist numerous algorithms for ab initio recognition of protein-coding regions and functional sites (reviewed in 5,6). The most popular gene prediction program for prokaryotes, GeneMark (6), utilizes non-homogeneous Markov models to find DNA regions that code for proteins or are complementary to them. Non-coding regions are described by homogeneous Markov models. A Bayesian decision rule is applied to deduce the coding capacity of sliding windows. GeneMark has been used in several genome sequencing projects (e.g. 7-9). Recently this algorithm was extended to take into account information about candidate ribosomal binding sites (10). The recently developed GLIMMER (11) has been reported to provide very high gene prediction accuracy in Haemophilus influenzae and Helicobacter pylori genomes. GLIMMER relies on interpolated Markov models to take into account DNA oligomers of varying length, dependent on the local composition of the sequence. Another program, EcoParse (12), utilizes hidden Markov models to find the maximum likelihood parse of a DNA sequence into coding and non-coding regions without the use of sliding windows. A program for gene recognition relying solely on ORF length and RBS was described by Hatzigeorgiou and Fickett (13).

Extrinsic analysis involves similarity searches with candidate gene products against protein sequence databanks. The most popular program of this class, BLASTX (14), performs six-frame translation of the query DNA and compares the resulting amino acid sequence to known proteins. Search results are represented in an integrated report, with hits from different reading frames combined to produce one statistically meaningful similarity score. Robinson et al. (15) used BLASTX to detect 450 new bacterial genes missed in original publications, including several genes previously known only in eukaryotes. Another DNA-protein search program, DPS (16), is the only currently available software tool that allows us to compare a complete genome sequence (3-5 mbp and more) with the total protein sequence databank in one pass.

Neither extrinsic nor intrinsic methods taken separately can ensure successful prediction. Practical experience in prokaryotic genome analysis as well as the recent trends in gene recognition in higher organisms (17,18) show that it is necessary to incorporate all available evidence in order to achieve reliable results. In real life, putative coding regions predicted by intrinsic methods are verified by similarity searches. Finding a related protein serves as a decisive supporting evidence. Pearson et al. (19) studied the ORFs predicted with GeneMark in H.influenzae, Methanococcus jannaschii and Mycoplasma genitalium. In many cases they were able to correct the length of genes based on comparative analysis with known proteins. Additionally they found many short genes not identified by GeneMark. The overall conclusion of this work is that a sizeable amount of genes annotated within the framework of large-scale sequencing projects are fully or partially wrong.

Experience coming from many computational genome analysis efforts (20-22) shows that 60-80% of genes in newly sequenced organisms have known counterparts in other species. In many cases the similarity is only marginal, partial or to gene products without known function. However, in at least 30% of the cases, reliable global alignments with well characterized proteins can be obtained. We were thus tempted to invert the usual procedure in which genes predicted by ab initio statistical methods are accepted or rejected based on subsequent similarity searches. In this work we use DNA regions significantly related to known proteins to extract codon usage statistics and other intrinsic recognition parameters that are further applied to unexplored parts of a genome. The leading idea of this work is that extrinsic evidence should be given higher priority than intrinsic information.

We also pay specific attention to assignment of gene starts. This is important since 5[prime] ends of genes often are not conserved, whereas they carry important functional and structural information. In particular, a signal peptide may provide information about protein localization (23). The N-terminus can contain information about the life span of a protein (24). The estimated strength of the RBS can be an indicator of the efficacy of the translation initiation (25). Thus correct determination of the gene start can be as important as identification of the gene itself.

MATERIALS AND METHODS

Data

Complete nucleotide sequences of the Bacillus subtilis (9) and Escherichia coli (8) genomes were downloaded from the SubtiList WWW Server at the Pasteur Institute (http://www.pasteur.fr/Bio/SubtiList.html ) and the E.coli genome project resource at the University of Wisconsin (http://ecoliftp.genetics.wisc.edu/ ), respectively. In addition, we also obtained full sets of the protein sequences encoded in these two genomes (4099 for B.subtilis and 4277 for E.coli) as assigned by the genome authors based on the application of various computational techniques as well as manual analysis.

Independently, B.subtilis and E.coli sequences were extracted from the PIR-International protein sequence database using the Sequence Retrieval System (SRS; 26). Special care was taken to select only sequences determined by individual researchers in different laboratories not associated with large-scale sequencing projects. Sequences submitted after 1995, plasmid sequences, fragments and proteins described in PIR as hypothetical, as well as the PIR entries containing the names of the main researchers involved in the B.subtilis and E.coli genome projects were discarded. These two sets were compared with the full sets of gene products from the two genomes using the BLAST2 software (W.Gish, unpublished; 27,14). Only the PIR sequences at least 98% identical to their counterparts in complete genomes and having the same N-terminal sequence were retained. This selection procedure resulted in 219 E.coli and 346 B.subtilis proteins.

Throughout this work, the full sets of gene products from complete genomes as determined by the authors will be referred to as SUBGEN (for B.subtilis) and ECOGEN (for E.coli), and the sequence sets extracted from the PIR database as SUBPIR and ECOPIR, respectively.

For similarity searches we created a non-redundant protein databank by merging the PIR-International (28), SWISS-PROT, SWISS-NEW, TREMBL and TREMBLNEW (29) sequence collections using the NRDB2 software developed by W.Gish (unpublished). Sequences from all species for which genome sequencing projects have been completed were excluded. The resulting databank currently contains 208 660 protein sequences.

Outline of the algorithm

Our algorithm is based on the assumption that information about coding regions derived from similarity searches is in principle more reliable than statistical data. We use the term `seed ORF' to describe the minimal, most reliable possible ORF that can be inferred. In the case of similarity searches, a seed ORF is obtained by extending the reliably aligned region in the upstream direction until the first start codon occurs and in the downstream direction until a stop codon is encountered. These similarity-derived seed ORFs are used to calculate coding potential parameters. For ORFs predicted ab initio a seed ORF results from extending a DNA region of a given minimal length (e.g., 300 nt) possessing sufficiently high coding potential (see below) in the same fashion. At the next step of analysis the algorithm tries to extend the seed ORFs by including additional upstream DNA fragments encompassing the next available start codon provided that the DNA region between the old and new candidate starts satisfies conditions imposed on coding potential. The sample of ORFs with a single possible start codon is used to derive the RBS recognition matrix. Finally, in ORFs with multiple candidate starts, the leftmost start codon having sufficiently strong RBS is selected.

Similarity search and the set of seed ORFs

We used the DPS program (16) to compare complete genomic sequences with the complete non-redundant protein sequence databank. DPS performs mapping of all protein sequences from the database onto the query genomic sequence. The DPS output contains full information about a DNA-protein match, including the start and end positions, reading frame, similarity score and alignment of the high-scoring DNA segment with the corresponding protein sequence fragment represented in three-letter code. The alignment may be split into several high-scoring fragments, in which case reading frames, coordinates and similarity scores are given for each such fragment, and the aggregate similarity for the entire alignment is indicated separately. We took into account only DPS hits with sufficiently high aggregate scores (typically >750) involving only one reading frame; cases involving more than one reading frame are subject to a separate procedure aimed at detecting frame-shifts.

Coding potential and the complete set of candidate ORFs

Seed ORF sets produced by the similarity search were utilized to calculate the codon usage tables and the average and standard deviation of the coding potential. To do that, the significant DNA-protein alignments were extended until the first stop codon downstream and the first start codon upstream occurred. The obtained DNA sequences are the most reliable representatives of the coding parts of the genome that can be extracted automatically.

Let F(abc) be the genomic frequency of the codon abc. Statistical weight of abc is defined as W(abc) = log F(abc). Primary coding potential of a DNA segment of length n codons is its log-likelihood:

To account for DNA fragments of different length, we will use the normalized potential measured in the standard deviation units:

where µ is the average codon weight which depends on the base composition of the genome and [sigma] is the standard deviation:

Finally, to avoid the influence of the local base composition and gene shadows, and to set the strand and reading frame, we define the coding quality of a DNA fragment as

Upon derivation of the statistical parameters above, the DPS output was screened again to extract all similarity-based seed ORFs. The parts of the genome not covered by the similarity-based seed ORFs were subsequently analyzed for the presence of other protein-coding seed ORFs. A seed ORF was accepted if its length exceeded a given threshold (100 codons) and its coding quality [Omega] was sufficiently high.

All seed ORFs were then extended in the 5[prime] direction as far as possible. Short overlaps between genes (up to 6 nt) were allowed. Each extension piece of DNA started with ATG, GTG or TTG. The extension was accepted if the coding quality [Omega] of the DNA segment of length 99 nt starting with the new candidate start codon was acceptable; otherwise the extension of a given ORF was interrupted. The window length of 99 was chosen to ensure sufficient statistical significance of the calculated coding potential (30).

This procedure resulted in the complete set of `open-start' candidate ORFs. The start codons for this set were assigned at the final step.

RBS weight matrix and assignment of gene starts

Candidate ORFs with only one possible start codon and not having neighbors closer than 30 bases upstream were selected. Regions (-20)...(+3) of these ORFs were aligned at start codons. These sequences were used to derive the RBS weight matrix.

Let L be the expected length of the SD box (L = 6). Denote by F(b,j) positional nucleotide frequencies in the initial alignment [j = (-20)...(-1); b = T,C,A,G]. Positional information content is:

G(b) is the genomic frequency of the base b (31).

Initially the RBS signal was assumed to reside in positions having the maximum total information content:

Then the position of the SD box in each individual sequence was determined using the following two-stage procedure.

We start with some definitions. Denote by N(b,k) positional nucleotide counts in the SD profile at a given iteration (k = 1...L). Positional nucleotide weights are:

max N(b,k) is the frequency of the consensus nucleotide in the position k. The SD signal score is calculated by the formula:

The first re-alignment stage involves the iteration until convergence of the following two steps: (i) find in each sequence the segment of length L with the highest score; (ii) re-calculate the nucleotide weight matrix.

A distinctive feature of our algorithm is that at each optimization step only the top scoring fraction (usually top 80%) of sequences are used to produce the current weight matrix.

At the second stage the preferences for the distance between the SD box and the start codon are taken into account. Let M be a possible position of the SD box within the RBS region, and let this position occur N(M) times. Denote by N_max the count of the most frequent position. The positional weights are calculated using the standard formula:

Now the strength of the SD signal is defined as:

Thus the RBS profile is the nucleotide weight matrix and the vector of position weights. The two step iterative procedure is used again until convergence.

The final step of the genome annotation is the assignment of start codons to `open start' candidate ORFs. If a candidate ORF contains start codons with sufficiently strong RBS, the 5[prime]-proximal of these starts was accepted. Otherwise the initial start generated at the previous stage was used, thus taking into account the possibility of translation re-initialization from an upstream gene.

Choice of the minimal allowed ORF length and handling of short ORFs

First versions of our program worked with a fixed minimal ORF length, typically 100 codons. The conflicts between overlapping ORFs were resolved based on the strength of the coding potential, as described above. The disadvantage of this approach was that quite often short ORFs defeated much longer competing ORFs, and the genome regions vacated by the latter would be returned to the pool of unoccupied space, giving rise to additional abundant short ORFs. This problem became especially severe when the minimal allowed ORF length was set to values under 100 codons, which led to a large number of predicted short ORFs at the expense of the longer ones.

To resolve this problem, we modified the final stage of gene prediction process as follows. The minimal allowed ORF length is first set to a very high value (2000 nt), after which all genes are predicted as described above. Then the minimal ORF length is reduced step-wise, and the gene prediction process repeated. At that, the genes predicted at the previous step remain unchanged. Thus, longer ORFs get higher priority, and the next pool of ORFs is derived from the genome regions that are unoccupied after completion of the previous step. This allows us to avoid the explosion of short ORFs while preserving high overall prediction accuracy.

Table 1. RBS weight matrix for B.subtilis and E.coli

Nucleotide Nucleotide position in the window

1 2 3 4 5 6

Bacillus subtilis

A 0.000 -0.909 -0.845 0.000 -0.909 -0.511

C -0.856 -1.000 -0.943 -0.923 -0.999 -0.748

G -0.804 0.000 0.000 -0.709 0.000 0.000

T -0.765 -0.897 -0.980 -0.868 -0.962 0.725

Consensus A G G A G G

Escherichia coli

A 0.000 -0.035 0.000 -0.995 -0.936 0.000

C -0.299 0.000 -0.804 -0.984 -0.987 -0.814

G -0.386 -0.115 -0.903 0.000 0.000 -0.710

T -0.027 -0.426 -0.752 -0.937 -1.009 -0.749

Consensus A/T C/A A G G A

Implementation and availability

The program to calculate weight matrices based on a multiple alignment of putative RBS regions is called STARTERand is written in C programming language. All other computational steps described in this paper are implemented as a Perl 5 script called ORPHEUS. Both programs as well as all data mentioned (protein sequence sets, weight matrices, DPS search results, etc.) are freely available to academic users; see http://pedant.mips.biochem.mpg.de/frishman/orpheus_home.html

Nucleic Acids Research