Bioinformatics Advance Access originally published online on October 25, 2006
Bioinformatics 2006 22(24):3003-3008; doi:10.1093/bioinformatics/btl544
Hypothesis testing approaches to the exon prediction problem
1 Statistics Department, University of Barcelona Barcelona, Spain
2 Genetics Unit, Department of Experimental and Health Sciences, Pompeu Fabra University Barcelona, Spain
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
Motivation: Many gene identification methods assign scores to gene elements prior to their assembly into predicted genes. The scoring system is often based on log-likelihood ratios. These methods usually perform well but it is difficult to interpret how significant a score is.
Results: We have developed several tests of significance for the scores: (1) a sum-of-scores test (SST), (2) an intersection-union test (IUT), based on a multiple hypothesis testing interpretation of an exon's score and (3) a meta-analytical approach (MA), which combines several P-values, corresponding to the exon's parts, to yield a global P-value. We performed simulation studies, which show that the MA has better sensitivity and specificity than other methods and is easier to interpret by non-expert users. This is an improvement over other methods and is especially relevant for users who would like to predict incomplete gene sequences.
Contact: asanchez{at}ub.edu
Supplementary information: Supplementary data are available at Bioinformatics online.
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
The gene identification problem can be formulated as the deduction of the amino acid sequences encoded in a given DNA genomic sequence (Guigo, 1997; http://genetics.mgh.harvard.edu/doc/genefinder.doc.html). This is an important but difficult problem, especially in eukaryotes, where genes are often split into exons separated by introns. The beginning or end of these fragments is specified in the genome sequence by a few types of signals encoded in the primary DNA sequence (see Table 1), which instruct the cellular machinery in the pathway from DNA to protein. Signals involved in gene specification are all encoded in the primary DNA sequence. However, they are ill-defined, lack generality and are highly unspecific. Using currently available detection methods, it is usually impossible to distinguish signals truly processed by the cellular machinery from those actually not functional. Attempting to predict gene structure by processing only DNA sequence signals often results in a computationally intractable combinatorial explosion of potential products. As a result, any gene prediction method that relies on these signals has to be able to distinguish between false exons, i.e. sequences delimited by the signals but not producing any transcript, and true exons, which have both the signals and the capacity to produce transcripts.
|
The main objective of this study was to derive and compare methods to improve the performance of some gene prediction methods, such as geneid (Parra et al., 2000), which are based on a log-likelihood ratio scoring system. These methods have some drawbacks: scores are not bounded and their values depend on various factors such as exon type or exon length (Zhang, 1998). As a consequence, a high score is often not related to a high significance.
All the methods discussed in this study view the exon prediction step as a problem of hypothesis testing. This leads to the possibility of introducing a probabilistic scoring system based on P-values, which could be used as an additional input in the dynamic programming algorithms used to build the final gene prediction from the candidate exons. Other dynamic algorithms, such as genetic algorithms, could also be implemented using P-values. These latter could also be useful in detecting alternative splicing.
|
2 METHODS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
A commonly used (ab initio) approach to gene-prediction is to construct a set of potential exons, which are scored using a statistical model. These potential exons are highly overlapping and only a few are assumed to form part of a gene. The subset of exons that will be assembled into predicted genes are selected using some type of dynamic programming algorithm, which seeks to maximize the total score assigned to the gene.
The geneid (Parra et al., 2000) or GeneFinder [Green,P., (2001) unpublished data, http://genetics.mgh.harvard.edu/doc/genefinder.doc.html] programs use this approach and score potential exons using log-likelihood ratio scores. HMM based programs, such as GENSCAN (Burge and Karlin, 1997) also rely on the use of this type of scores (see Barret et al., 1997).
In simple terms, a gene, G, can be considered as a sequence formed by a series of one (G = e) or more exons separated by zero or more introns 
. As shown in Table 1, any exon is formed, as indicated, by three parts: Begin (B) site (start or acceptor site), Coding (M) part and End (E) site (donor or terminal site):
|
|
2.1 Geneid scoring model
The scoring system used by geneid assigns values to each part separately and then defines the exon score as a weighted sum of the log-likelihood of the scores of each component, then it can also be considered as a log-likelihood ratio score [LLRs, see Parra et al., 2000).
 | (1) |
- Begin and End sites are scored by assigning values to each base at each position of the sequence according to their match to the ‘consensus’ sequence (see Table 1) for its signal type. This is done by using positional weight matrices [(PWM, see Staden, 1984) whose elements represent the logarithms of the ratios between two probabilities (or two likelihoods): the probability
of finding a given base i at a given position j in a site A (either Begin or End) versus the probability 
of finding base i at position j in a random sequence. LB and LE represent the sum of these logarithms for Begin and End sites, respectively, and can also be seen as LLRs, as they are sum of log-likelihood ratios.
- Begin and End sites are scored assuming independence between successive bases. The scoring of coding parts of the exon uses a 5th order Markov chain to model the relation between successive bases in the sequence (see Ewans and Grant, 2001). Log-likelihood ratios are computed, analogously to Begin and End sites, by dividing the probability computed according to the Markov chain by the corresponding random probability (see Parra et al., 2000 for details).
- The final score is a weighted sum of sites and coding potential scores (see Parra et al., 2000) where wB, wM and wE represent the weights given to each score. This ensures that when the coding part is too short, the scores corresponding to the initial and end parts do not have too great an influence on the score. The weights are obtained empirically from a set of training sequences, (which change from one genome to another) and are chosen to ensure the maximum discrimination between known true and known false exons (see www1.imim.es/software/geneid/docs/chapter8/index.htmlwww1.imim.es/software/geneid/docs/chapter8/index.html). The values of the weights for this study, corresponding to the Drosophila melanogaster genome, were 0.6 for Begin and End site scores and 0.4 for Coding potential scores.
2.2 Tests of exonicity
Statistically, the LLRs scoring method has an obvious drawback: a high score does not necessarily mean a significant score. Although, there are theoretical results on the assessment of the significance of LLR-based scores [see Huang et al., 2004] these cannot be easily applied to complex score systems, such as Equation (1).
This study considers three different ways of deciding whether, given an exon, e, which is scored as L(e), this value should be considered statistically as greater than zero i.e. significant.
All three methods require the sampling distribution of the test statistic under the null hypothesis to be known. This has been approximated using two different approaches:
- Resampling. The empirical distribution is approximated by taking samples with replacement from a large list of scores computed on a set of known false exons. The list was obtained from a larger sequence used in a study by Guigo et al. (2000) and is available at http://www1.imim.es/databases/gpeval2000/http://www1.imim.es/databases/gpeval2000/. It contained 22 041 starts, 76 728 acceptors, 64 566 donors, 54 919 stops and 6 811 740 coding fragments. Guigo et al. (2000) also describe how the set of potential exons for these models is constructed.
- Monte Carlo. The empirical distribution is approximated by scoring a large list of sequences generated using estimated PWM for Begin and End sites and 5th order Markov chain models for the Coding parts, as described above. Both PWM and Markov chain models were estimated from the same reference dataset used for resampling.
2.2.1 Sum-of-scores tests (SST)
The SST is a straightforward approach that uses the scores as test statistics. It requires the distribution of the sum of the scores under the null hypothesis (i.e. in case the scored potential exon is false) to be known. This method yields approximate P-values of the significance test, allowing standard scores to be converted into probabilistic scores, and is easy to interpret. Although LLR depends on the length of the sequence, the SST works without exon length differentiation. We propose SST2, which assigns P-values accounting for exon length.
The hypothesis test can be stated as:
 | (2) |
, where 
is the 1 – 
-quantile of the null (reference) distribution to be approximated by simulation. The SST test does not consider the lengths of the exons. This problem has been overcome by introducing SST2 where null hypothesis and test statistics are the same as in SST. When simulating the null hypothesis distribution, exon length is accounted for by resampling each exon from a set of exons of similar length (see more details in section 2.3).
The Monte Carlo method is not applied to SST because it builds exons with the same length as the candidate.
2.2.2 Intersection-union tests (IUT)
In IUT, the components of the score are first considered separately and a global test is then performed using the intersection-union method [see Casella and Berger, 2002]. This seems to be the natural translation of the question ‘is this an exon?’ into a test, but has the drawback that a unique P-value cannot be obtained and therefore no probabilistic scoring is available. We consider an exon to have three components: Begin site, Coding part and End site. Given a three part sequence to be tested, the decision that it is an exon is taken if it is decided that the first part is a Begin site, the second part is Coding and the third part is an End site. Each decision can be seen as the alternative hypothesis for a null hypothesis of a random sequence, which gives:
 | (3) |
 | (4) |
 | (5) |
, 
and 
are the respective 1– quantiles of the null hypothesis distributions, which can be approximated by simulation. As a consequence of the properties of intersection union-tests [see Casella and Berger, 2002] the global significance level of the IUT test, with every single test of level , is, at most, .
2.2.3 Meta-analysis based tests (MA)
Appealing as the IUT test may be, it seems reasonable to assume that it may lack sensitivity (as shown in the next section) due to the strict requirements that must be fulfilled to reject the null hypothesis. An alternative approach consists of combining P-values obtained separately for each exon's component using MA methods, which were originally developed to test the statistical significance of combined results [see Hedges and Olkin, 1985]. MA are based on the idea that relatively small P-values, not necessarily significant, do not appear at random. It is known that, under the null-hypothesis, P-values are distributed uniformly. If P-values are consistently small, even if they are not significant, the uniform law could be discarded and this could be interpreted as evidence favoring the rejection of the null-hypothesis.
Of the various methods for combining P-values those included in this study are among the most widely used. They are based on known reference distributions (implying a low computational cost as the distribution is tabulated) and all satisfy the monotonicity principle [see Hedges and Olkin, 1985] and will thus be optimal in most cases. In addition, all three methods account for exon length. The null hypothesis is the same for all three tests considered:
 | (6) |
- MA1: Fisher's Chi-square method The test statistic is defined as:
and the null hypothesis is rejected if:
(7)
where
(8) 
is the 1-alpha quantile of a chi-squared distribution with 2k degrees of freedom.
- MA2: Logit method The test statistic is defined as:
where
(9) 
and 
. The null hypothesis is rejected if:
where
(10) 
is the 1-alpha quantile of a Student's t-distribution with 5k + 4 degrees of freedom.
- MA3: Inverse gaussian method The test statistic is defined as:
(11)
, where z is the 1-alpha quantile of a standard Normal [N(0,1)] distribution. All tests considered in this study were computed using an alpha value = 0.10, instead of the usual 0.05 in order to decrease the probability of false negatives, i.e. in order to avoid rejecting true exons unnecessarily.
2.3 Simulation experiments
Simulation experiments were designed to study and compare the performance of all tests. First, the number of simulation replicates was determined empirically, relying on the situation with the highest variability, which occurs when resampling coding regions between 5 and 104 bp long. The empirical procedure to determine the number of replicates consisted of:
- Selecting a sample of (205) real exons with a length between 5 and 104 bp.
- Repeatedly take 100 resamples of sizes 100, 200, 500, 1000, 2000, 5000, 10 000, 15 000, 20 000 and 50 000, from the false coding set having the same length as the corresponding real exon.
- For each resample R of size N we approximated the P-value by the percentage of false exons scoring equal or more than the corresponding true exon.
The criterion was to choose a size such that P-values did not change appreciably from this size on and this led us to use 15 000 as the appropriate resampling size for this study.
Ideally any gene identification method should have a type II error as low as possible in order to yield high sensitivity. In a real application, there are many more negative samples than real exons and this was dealt with by using 10 times as many false exons as real ones in all comparisons.
2.3.1 Experimental design
Test performance was evaluated by considering the following factors for each method: test (SST1, SST2, IUT, MAi), exon type (First, Internal, Terminal, Single), exon length (arbitrarily split into: <5, 5–104, 105–504, >504), P-value generation method (Resampling and Monte Carlo) and testing data (known false exons and known true exons). More details of the simulations can be found in Vilardell and Sanchez-Pla (2004), http://www.imub.ub.es/publications/preprints/pdf/preprint353.pdf.
The results are presented as tables containing the percentage of times that the null hypothesis that a given sequence is not an exon is rejected during each simulation and in each combination of factors, i.e. the percentage of times that the sequence is identified as an exon. This percentage corresponds to the sensitivity for real exons and to specificity for false exons. A good test should have a high percentage in tables for real exons (i.e. high-sensitivity) and for false exons (i.e. high-specificity).
The general appropriateness of all tests was evaluated using receiver operator characteristic (ROC) curves. ROC were calculated to compare resampling and Monte Carlo methods. Sensitivity and specificity values for detecting sites and the coding part at the same alpha level were calculated (see Table 6). ROC curves corresponding to Tables 2–6 can be found in the Supplementary material. We also compared the resampling distribution of sites and the coding part with the simulated ones using variability graphs, which are shown in the Supplementary material.
All the methods considered here can be applied to any genome. This study used only Drosophila sequences. However, a wider study using other genomes is underway.
|
3 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
The Monte Carlo simulation as a method of generating P-values is similar to the resampling method for site models but differences appeared in coding part models (see Supplementary material). It systematically yielded lower values in most combination of factors affecting sensitivity and specificity. ROC curves (see Figs 1 and 2) illustrate this when different tests are compared. For this reason, only results obtained with the resampling method are presented.
|
|
The main results of the simulation study are presented in Tables 2–6 and can be summarized as follows:
- The IUT test performed worst, with a lower sensitivity than the SSTs (1,2) and MA tests (1,2,3).
- SST seems to be less specific than SST2. SST tends to lose specificity with exon length (Table 3 or Fig. 1) tending to classify large potential exons as true exons. SST2 is more consistent and yields similar specificity values for all lengths (Table 3). Therefore, we propose using SST2 rather than SST.
- The SST2 test performs very similarly to all three MA tests in terms of sensitivity. However globally it is less specific than any of the MA tests. SST2 has a significant loss of specificity in classes >504 bp. In general, we observed that, as the exon length increased, the sensitivity of SST2 tended to be slightly higher than that of any of the MA tests, while its specificity tended to be lower (cf. Tables 2 and 3).
|
|
The performance of the different tests in the various conditions analyzed showed some well-defined trends:
- Most tests performed better for internal exons and worse for single ones. The internal exon was the exon type with the strongest signals and the single exon that with the weakest (see Table 6). In single exons, signals were globally consistent with random distribution in P-values. Here, MA test fails by definition. An alternative model for coding parts based on Markov Chain Models should be considered for single exons.
- Fisher's test (MA1) performed slightly better than MA2 and MA3 globally, although it was less robust to the lack of adjustment than MA2 and MA3 (see single exons in Table 4). If there were an alternative model for coding part in single exons MA1 would be more powerful and more differences among MA1, MA2 and MA3 could be found.
|
In summary, MA-based tests performed similarly to the SST test, with the advantage that no additional parameter, such as weights in SST tests, needed to be estimated. Fisher's chi-square (MA1) test was slightly better than the Logit (MA2) or Gauss (MA3) tests in terms of higher sensitivity and specificity. It performed better than the other for most combinations of factors and would thus be the test of choice to obtain a probabilistic score.
|
4 DISCUSSION AND CONCLUSIONS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
This study tries to formulate an approach to the meaning of the scores as test statistics (see Vilardell and Sanchez-Pla, 2004), which led us to introduce the IUT. This test allows the question ‘Is a sequence, S, an exon?’ to be recast as a test of hypothesis, which we have called a test of exonicity and seems to be a natural choice because the question is broken down into three individual tests: (1) ‘Is the initial fragment a Begin site?’ (2) ‘Is the internal sequence Coding?’ (3) ‘Is the final fragment an End site?’ The final decision comes from rejecting the three null hypotheses. Appealing as it may be, the test has been shown to perform badly in simulation studies. The nature of this test implies that all three hypotheses must be rejected in order to decide that the sequence is an exon, and the fact that the signals are so weak (see Table 6) suggests that it is very easy for this not to occur in at least one of the three tests.
Although a testing approach method based on IUT seemed reasonable, the weakness of the signals meant that it did not work well enough. One further possibility was to try to combine the different P-values instead of using them separately. Techniques for combining P-values have been developed by several authors and are often used in meta-analysis [Hedges and Olkin, (1985)]. We chose three of these techniques and carried out three different tests based on combining P-values for tests (1), (2) and (3) of the previous paragraph. All three tests performed well high-sensitivity and high-specificity in simulation studies.
However, we would recommend Fisher's method to improve the selection of exon candidates as it is easier to apply and performed slightly better than the other MA tests. Since exact P-values could be calculated from background models with low computational costs [see Huang et al. (2004)] and although better background models in the coding part are needed to avoid the use of resampling methods, Fisher's method could be used to develop methods to assign exact P-values for each potential exon. In addition, background models take exon length into account, and are not biased towards large exons, as occurs with SST. Exons longer than 504 bp present a classification problem, which is directly related to the difficulty in classifying single type exons, as this class tends to have longer exons, as shown in Table 5. We suggest studying single exons in more detail in order to build specific models for them.
|
|
In summary, we have introduced a new method for detecting coding sequences. This method performs at least as well as LLR methods but does not require estimation of additional parameters, such as in Parra et al., (2000). In addition, this method is easier to interpret than LLR because it yields P-values, which can be interpreted as probabilistic scores ranging from 0 to 1, making its interpretation more intuitive than is the case with the arbitrary scales used in LLR scores.
This method could be extended in several directions. For example, the resulting P-values could be used with genetic algorithms to detect alternative splicing or to compare genomes.
|
Acknowledgments |
|---|
The authors thank Roderic Guigo, Genis Parra and Enrique Blanco from Genome Bioinformatics for their support and helpful comments. The authors also thank Lauro Sumoy, Roberto Castello and Josep Lluis Mosquera for reviewing the manuscript.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Dmitrij Frishman
Received on May 18, 2006; revised on October 18, 2006; accepted on October 18, 2006
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSION AND CONCLUSIONSREFERENCES |
|---|
Barrett, C., et al. (1997) Scoring Hidden Markov Models. Comput. Appl. Biosci. Apr, . 13, 191–199.
Burge, C. and Karlin, S. (1997) Prediction of complete gene structures in human genomic DNA. J. Mol. Biol, . 268, 78–94[CrossRef][Web of Science][Medline].
Casella, R. and Berger, G. Statistical Inference, (2002) 2nd edn , Pacific Grove, CA Duxbury Press.
Ewans, W.J. and Grant, G.R. Statistical Methods in Bioinformatics: An Introduction, (2001) , NY Springer.
Green, P. (2001) GENEFINDER documentation.
Guigó, R. (1997) Computational gene identification. J. Mol. Med, . 75, 389–397[CrossRef][Web of Science][Medline].
Guigó, R., et al. (2000) An assessment of gene prediction accuracy in large DNA sequences. Genome Res, . 10, 1631–1642
Hedges, L.V. and Olkin, I. Statistical Methods for Meta-Analysis, (1985) , CA Academic Press.
Huang, H., et al. (2004) Determination of local statistical significance of patterns in Markov sequences with application to promoter element identification. J. Comput. Biol, . 11, 1–14[CrossRef][Web of Science][Medline].
Parra, G., et al. (2000) geneid in Drosophila. Genome Res, . 10, 511–515
Staden, R. (1984) Computer methods to locate signals in nucleic acid sequences. Nucleic Acids Res, . 12, 505–519
Stormo, G. (1990) Consensus patterns in DNA sequences. Meth. Enzymol, . 183, 211–221[Web of Science][Medline].
Vilardell, M. and Sánchez-Pla, A. The Gene Identification Problem from a Hypotheses Test Perspective, (2004) Institute of Mathematics, University of Barcelona 353, .
Zhang, M.Q. (1998) Statistical features of human exons and their flanking regions. Hum. Mol. Genet, . 7, 919–932
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||