Bioinformatics Advance Access originally published online on April 19, 2005
Bioinformatics 2005 21(13):3051-3052; doi:10.1093/bioinformatics/bti451
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
S-SPatt: simple statistics for patterns on Markov chains
Laboratoire Statistique et Génome 523 place des terrasses de l'Agora, 91000 Evry, France
 |
Abstract |
|---|
 TOP Abstract 1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
Summary: S-SPatt allows the counting of patterns occurrences in text files and, assuming these texts are generated from a random Markovian source, the computation of the P-value of a given observation using a simple binomial approximation.
Availability: S-SPatt is available at: http://stat.genopole.cnrs.fr/spatt
Contact: spatt{at}genopole.cnrs.fr
|
1 PREVIOUS WORK |
|---|
|
TOPAbstract1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
It is well known that patterns frequencies are highly dependent on their overlapping structures as well as on the composition bias in the considered sequences. In order to compare patterns efficiently, we need more than their simple counts. A common solution to this problem consists in using a Markov chain to modelize the sequence thereby ranking patterns according to their P-values.
Several software propose solutions to compute these P-values: QuickScore using exact computation with generative series as proposed in Régnier (2000) (under development but not yet available), RMES (Schbath, 1997) for Gaussian and Compound Poisson approximations and LD-SPatt for large deviations approximations (Nuel, 2004). Unfortunately, all these methods suffer severe drawbacks: QuickScore can not consider Markovian model of order higher than 1 or 2 and is limited to short sequences (for exact computations), Gaussian approximations are wrong for tail distribution events (which are of course the events of interest), Compound Poisson approximations are designed only for rare patterns (and have some numerical issues in the RMES implementation making its use difficult anyway) and, finally, large deviations techniques are limited to short patterns (less than 10 letters long on a DNA alphabet, for example).
The Regulatory Sequence Analysis tool (RSAT) proposes to use simple binomial approximations to compute these P-values. Simulations have shown (van Helden et al., 1998) that such approximations are wrong but very close to the optimal solution when the considered patterns are not self-overlapping too much.
As the RSAT package is only available on-line, we decided to implement the same method in a stand-alone GPL program called S-SPatt to compare its overall performances with other similar programs.
|
2 METHOD |
|---|
|
TOPAbstract1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
2.1 Statistics
Let X = X1, ..., Xn be an order m stationary Markov chain on a size k finite alphabet
, with 
(km x km dimension) sparse transition matrix and µ (km dimension vector) as stationary distribution.
We consider a pattern 
as a set of r words (of respective length h1, ..., hr).
We count the number of occurrences of a given word W = w1, ..., wh using N(W) as defined by
 | (1) |
According to the model, Yi has a Bernoulli distribution of parameter
 | (2) |
for all 1 
i j h.
Couples (Yi, Yj) are clearly not independent nevertheless, we can use the following heuristic distribution for a single word:
 | (3) |
 | (4) |
2.2 Algorithms
2.2.1 Stationary distribution
According to the theorem of Perron–Frobénius, for any irreducible Markov chain, the stationary distribution can be computed by solving an eigenvalue problem. Growing exponentially with the Markov order, the scale of this problem is usually very large. To solve this problem, we propose to use an explicitly restarted Arnoldi's algorithm (Stewart, 1994) which is known to be very efficient with large, sparse eigenproblems (for an order m Markov model we have a non zero density of k1–m).
2.2.2 Number of occurrences
All words of length smaller than a given L are first computed (n x L in time and memory) and deterministic finite state automata are used for larger patterns (Hopcroft and Ullman, 1979).
2.2.3 P-values
P-value for the observations are computed with the incomplete Beta function (Press et al., 1992).
|
3 RESULTS |
|---|
|
TOPAbstract1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
3.1 Numerical performances
The following table gives the computation time T (in seconds) for N computations of simple words' P-values (all words of length h = 6, ..., 10 in Mycoplasma genitalium with an order 1 Markov model; computations are performed on an Intel Pentium 4 processor at 2.8 Gz).
A simple linear regression gives about 200 000 computations per second which is very fast. Let us add, as a remark, that this result is independent of the Markov model order (thanks to the Arnoldi algorithm). In comparison, for words of length 8 on a DNA alphabet, RMES Gaussian performs about 30 000 computations per second, RMES compound Poisson 
2000 and LD-SPatt 4 (computation time with LD-SPatt grows exponentially with the length of the patterns).
S-SPatt can also be used to compute very efficiently simple word frequencies. For example, counting all words of length 8 on Escherichia coli complete genome requires roughly 15 min with the wordcount program from the popular EMBOSS package while S-SPatt can achieve the same task in less than half a second (hence, S-SPatt is 2000 times faster than EMBOSS).
3.2 Reliability of the heuristic
As Gaussian approximations are expected to be good in the center of the distribution (for high P-values), they are taken as reference for such events while large deviations are used as reference for all tail distribution events (see Nuel, 2004 for more discussion on the subject).
We can see from Table 1 that our simple binomial approximation is the closest to the reference both in terms of relative error and, more important, in terms of rank agreement (see Press et al., 1992 for more details on Kendall's Tau).
|
|
4 CONCLUSION |
|---|
|
TOPAbstract1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
S-SPatt is not only the fastest tool (about 200 000 computations per second) but also the most reliable one (after the reference ones).
Moreover, the proposed implementation have several interesting features:
- Support of any user defined alphabet (regular DNA, purin-pyrimidin, amino acids, group of amino acids, latin, case sensitive, ...) with a simple syntax.
- Markov model parameters are estimated using maximum likelihood or are specified by the user.
- Support for high order Markov chain including computation of the stationary distribution using efficient linear algebra methods.
In conclusion, S-SPatt seems to be a very good heuristic for the computation of pattern statistics on Markov chains.
Received on February 15, 2005; revised on March 29, 2005; accepted on April 13, 2005
|
REFERENCES |
|---|
|
TOPAbstract1 PREVIOUS WORK2 METHOD3 RESULTS4 CONCLUSIONREFERENCES |
|---|
van Helden, J., et al. (1998) Extracting regulatory sites from the upstream region of yeast genes by computational analysis of oligonucleotide frequencies. J. Mol. Biol., 281, 827–842[CrossRef][ISI][Medline].
Hopcroft, J.E. and Ullman, J.D. Introduction to Automata Theory, Languages, and Computation, (1979) , Reading, MA Addison-Wesley.
Nuel, G. (2004) LD-SPatt: large deviations statistics for patterns on Markov chains. J. Comp. Biol., 11, 1023–1033[CrossRef].
Press, W.H., et al. Numerical Recipes in C, (1992) Cambridge University Press.
Régnier, M. (2000) A unified approach to word occurrence probabilities. Discrete Appl. Math., 104, 259–280[CrossRef].
Schbath, S. (1997) An efficient statistic to detect over- and under-represented words in DNA sequences. J. Comp. Biol., 4, 189–192.
Stewart, W.J. Introduction to Numerical Solution to Markov Chains, (1994) Princeton University Press.
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||