Significativité d’un alignement pour BLAST

On a vu il y a quelque temps comment estimer la significativité d’un alignement en utilisant la formule de Bayes. Mais pratiquement, ce n’est pas la méthode la plus utilisée. Il est temps d’introduire BLAST…

Derrière cet acronyme percutant se cache l’algorithme le plus utilisé en alignement de séquences biologiques. Il correspond à une heuristique d’alignement local deux-à-deux présenté dans l’article « Basic local alignment search tool » dont il tire son nom. Cet article d’Altschul, Gish, Miller, Myers et Lipman, paru en 1990, a fait date: à ce jour il a été cité dans au moins 30000 publications, un record! Il est très utilisé parce qu’il est rapide et permet de trouver les alignements significatifs entre une sequence inconnue et une (très) grande banque de séquences connues.

Cet algorithme fonctionne en plusieurs étapes. Étant donné deux séquences à comparer, il commence par rechercher des segments de longueur égale, sans gap,  identiques entre eux. On appelle ça les seeds, ce sont des k-mers avec k=11 pour des séquences nucléotidiques. Puis l’algorithme essaie d’étendre ces seeds en autorisant les gaps: il obtient alors des HSPs (pour high-scoring segment pairs). Finalement, pour chaque HSP (donc pour chaque alignement local entre les deux séquences),  l’algorithme choisit ceux dont le score est considéré comme étant significatif.

Optimisation of the Smith-Waterman algorithm

Supposons que l’on aligne une séquence de fonction inconnue, la query, de taille m, avec une séquence de fonction connue, la subject, de taille n, via BLAST. On obtient N HSPs, chacun ayant un score H. Pour estimer la significativité d’un HSP, on a besoin de savoir à quel point c’est probable d’obtenir un HSP de score H par hasard, sans que les deux séquences aient quoi que ce soit à voir entre elles. Cela revient à calculer la probabilité que le score maximal des N HSPs obtenus dépasse un seuil S donné.

Maintenant considérons que les H_i sont des variables aléatoires indépendamment distribuées selon une densité de probabilité f(H).  Si on prend leur maximum X = \max_i ( H_i ), on peut écrire:

P_N ( X \le S ) = P(H_1 \le S) ... P(H_N \le S)

P_N ( X \le S ) = \left( \int_{-\infty}^{S} p(H) \, dH \right)^N

P_N ( X \le S ) = \left( 1 - \int_{S}^{+\infty} p(H) \, dH \right)^N

Si les deux séquences ont beaucoup d’HSPs entre elles (N grand) et comme S est vraisemblablement élevé, il est dans la queue de la distribution et donc la surface sous l’intégrale est petite. On peut donc appliquer l’approximation bien connue (1-x)^k \approx \exp^{-kx} quand x est petit:

P_N ( X \le S ) \approx \exp \left( - N \int_{S}^{+\infty} p(H) \, dH \right)

Maintenant, supposons pour l’instant que p(H) tombe exponentiellement sur ces extrêmes (les queues de distribution, quand S est suffisamment grand):

\int_{S}^{+\infty} p(H) dH = \int_{S}^{+\infty} a \exp^{-\lambda H} dH = \frac{a}{\lambda} \exp^{-\lambda S}

Ainsi: P_N ( X \le S ) \approx \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Finalement, la probabilité que le score maximal X des N HSPs obtenus entre nos deux séquences dépasse un seuil S donné vaut:

P_N ( X > S ) \approx 1 - \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Cette formule nous permet donc de dire quand le score X tombe dans les 1% ou 5% de la distribution. C’est bien ce qu’on cherchait, non ?

Pour la petite histoire, X suit une distribution dite « des valeurs extrêmes ». Vous pouvez trouver ici l’article de Gumbel, en français, décrivant les premiers résultats. On se sert de cette théorie dans de nombreux domaines, notamment en finance et assurance, pour estimer la variabilité d’un système et donc les risques d’atteindre des valeurs extrêmes. Les lois des valeurs extrêmes ont l’une de leur queues de distribution plus épaisse que la loi normale, donc la probabilité d’occurrence d’une valeur extrême, c’est-à-dire tombant dans cette queue de distribution plus épaisse, est plus forte que pour une loi normale. Mais certains aiment quand même bien la loi normale, et un prix Nobel ça inspire confiance, quoique après ça, certains vont réviser leur modèles…

Loi Normale et loi de GumbelLa 1e image est tirée de « The estimation of statistical parameters for local alignment score distribution », Althschul et al., Nucleic Acids Research, 2001. Quant à la 2e image, voici le code R pour la générer:


require( distrEx )
require( Cairo )
Cairo( "normal_vs_gumbel.png", width=700, height=500, bg="white" )
x <- seq(-4,6,0.01)
par( mar=c(5,5,3,2), font=2, font.axis=2, font.lab=2, cex=1.5, lwd=2 )
plot( x, dnorm(x), type="l", col="black", main="Loi Normale et loi de Gumbel", ylab="densité de probabilité", ylim=c(0.0,0.5), xlim=c(-4,6) )
points( x, dgumbel(x), type="l", col="red" )
abline( v=0, lty=2, lwd=1 )
legend( 0.2, 0.5, c( "densité de la loi Normale", "densité de la loi de Gumbel" ), col=c("black","red"), lwd=2, bty="n" )
dev.off()

Ajout: un tutoriel sur BLAST vient de sortir dans PLoS Biology, avec une figure qui résume très bien les différentes étapes.

2 commentaires pour Significativité d’un alignement pour BLAST

  1. Leïla dit :

    merci du petit cours très clair
    une question tout le monde utilise la même méthode dans le monde? a t-elle été comparée à d’autres méthodes?

    • walrus dit :

      Quand tu veux trouver LE meilleur alignement entre une séquence de taille M et une autre séquence de taille N, l’algorithme exact (de Needleman et Wunsch) a une complexité en temps O(MN). C’est le nombre d’instructions que le programme doit faire pour trouver le meilleur alignement. Pour des séquences de longue taille, ça peut durer quelque temps.
      Mais quand on doit comparer x séquences inconnues avec les y séquences connues stockées dans les bases de données, comme x et y sont très grands (x ~= 10^4; y ~= 10^9), c’est impossible en pratique d’utiliser l’algorithme exact. Il faut donc utiliser une heuristique, cad un algorithme qui te donne une solution approchée. BLAST est justement utilisée pour ça, tout le monde l’utilise donc. Beaucoup de travaux théoriques ont été menés pour savoir avec quelle probabilité on pouvait observer une séquence inconnue s’alignant avec une séquence connue. Des simulations avec des vraies données ont également été faites.
      Mais attention, BLAST ne fait qu’aligner une séquence avec une autre ! Un alignement de bonne qualité est important mais on s’est vite rendu compte que dans certains cas, par exemple avec les séquences d’ARN non traduites en protéines, la structure 2D et 3D est aussi importante que la séquence, et que les unes sont indirectement corrélées avec l’autre. Il a donc fallu mettre au point d’autres algorithmes, beaucoup plus sophistiqués, à base de modèles de Markov cachés par exemple… Mais ce sera (peut-être) pour un autre billet.

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion / Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion / Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion / Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion / Changer )

Connexion à %s

%d blogueurs aiment cette page :