Tests unitaires

27 juillet 2009

Un programme informatique, tant qu’il ne fait qu’imprimer « Hello World! » à l’écran, ça a l’air tout bête. Mais dès qu’on essaie d’en faire un peu plus, ça devient plus délicat… Alors vous imaginez ce qu’il en est lorsqu’on développe un programme informatique totalisant des milliers de lignes de code ! Par exemple voir le cas de BLAST pour les languages C/C++ ou Biopython pour le language Python.

A cela s’ajoute une difficulté de taille: le code évolue, tout le temps. De ce point de vue, un programme informatique est très proche, conceptuellement parlant, d’un organisme vivant. Il manipule de l’information, remplit une fonction bien précise et doit être performant dans un environnement donné. Malheureusement, à chaque fois qu’un développeur (le type qui pianote sur son clavier toute la journée…) met son nez dans le code pour y ajouter ou modifier quelque chose, il risque de casser autre chose (une mutation délétère en quelque sorte), et la plupart du temps sans s’en rendre compte.

crash test

Dans ces conditions, soit vous êtes comme Donald Knuth et vous savez écrire un programme dans les règles de l’Art qui compile du premier coup, soit vous êtes comme tout le monde et vous testez votre code. Les informaticiens étant des gens très pointilleux, vous imaginez bien qu’ils ont peaufiné la chose: une bonne façon de développer est de le faire en Test-Driven Development (TDD).

Pour ce faire, supposons que l’on veuille impimer « METHINKS IT IS LIKE A WEASEL » à l’écran (lire cet article pour comprendre l’origine de la phrase…). Disons que l’on utilise la programmation orientée objet pour développer et qu’on écrit le code en TDD. Pratiquement, on aura une classe « Example » ayant un attribut « message » correspondant à la phrase donnée ci-dessus et ayant une méthode « getMessage() » qui servira à renvoyer le message en attribut.

En bon élève, on commence par écrire le test. Voici ce que ça donne en Python:

import unittest
import Example
class Test_Example( unittest.TestCase ):
    def test_getMessage( self ):
        expected = "METHINKS IT IS LIKE A WEASEL"
        myObject = Example.Example();
        observed = myObject.getMessage();
        self.assertEqual( observed, expected )
test_suite = unittest.TestSuite()
test_suite.addTest( unittest.makeSuite( Test_Example ) )
if __name__ == '__main__':
    unittest.TextTestRunner(verbosity=2).run( test_suite )

On exécute ce bout de code et bien sûr le test ne passe pas puisqu’on n’a pas encore écrit la classe à tester. Ce qu’on s’empresse de faire:

class Example:
    def __init__( self ):
        self.message = "METHINKS IT IS LIKE A WEASEL"
    def getMessage( self ):
        return self.message
if __name__ == '__main__':
    main()

On relance le code de test et maintenant ça passe: il renvoie « ok ». On sait donc que notre méthode « getMessage() » fait bien ce qu’on lui demande de faire.

$ python Test_Example.py
test_getMessage (__main__.Test_Example) ... ok
----------------------------------------------------------------------
Ran 1 test in 0.001s
OK

Vu comme ça, ce billet de blog paraît un peu léger, je m’en doute. Mais il faut savoir que cette approche très « professionnelle » de la programmation (ExtremeProgramming, méthodes Agiles) est bien souvent absente des labos de bioinformatique. Un chercheur peut alors se retrouver à télécharger le dernier programme publié dans Bioinformatics sans que celui-ci ne fonctionne (en gros, on installe, on lance et ça plante). Mais le pire arrive lorsque le programme téléchargé tourne sans bug apparent mais ne fait pas exactement ce que ses auteurs prétendent qu’il fasse…

Pour éviter ça et faire évoluer votre code de manière modulaire et adaptive, faites des tests unitaires ;-). Vous verrez, on s’y met très vite, d’autant plus que cette façon de travailler est grandement simplifiée grâce aux bibliothèques de tests unitaires déjà existantes: unittest en Python et cppunit en C++. En C++, c’est moins évident de comprendre du premier coup comment intégrer cppunit dans son code. Si vous voulez un exemple qui marche (on en trouve plein sur le web mais à qui il manque toujours qqch pour que ça compile comme il faut), demandez-moi, je peux vous en fournir un…

Publicités

Significativité d’un alignement pour BLAST

12 juillet 2009

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.


Significativité d’un alignement en bayésien

14 mai 2009

On a vu précédemment comment aligner (globalement) deux séquences l’une avec l’autre et comment calculer le score de l’alignement résultant étant donné un scoring scheme (match, mismatch, gap). Mais bien que l’on sache obtenir un alignement optimal, cela ne nous dit pas si notre alignement est significatif, c’est-à-dire s’il est biologically meaningful

Replaçons nous dans le contexte; nous venons d’obtenir (par des moyens expérimentaux) une nouvelle séquence. Pour en savoir plus sur l’information qu’elle code, on se dit qu’on peut essayer de l’aligner avec une séquence déjà connue. C’est ce que nous faisons et on trouve l’alignement optimal A_1 de score S_1. C’est bien, mais pas suffisant. En effet, nous n’avons (pour l’instant) aucun moyen de savoir si cet alignement est dû au hasard ou bien si ces deux séquences sont homologues, c’est-à-dire si elles ont un ancêtre commun.

Et oui, on peut très bien avoir un alignement A_2 entre notre séquence inconnue et une 2e séquence que l’on connaît déjà, de score S_2 celui-là. Or on aimerait bien avoir un moyen pour juger que deux séquences sont homologues puisque dans ce cas-là on se permettra (grossièrement) d’inférer que la fonction de notre séquence inconnue est la même que la fonction de la séquence qu’on connaît déjà (à peu de choses près).  Dans notre cas, quel alignement doit-on utiliser pour inférer la fonction ? A_1 ou A_2 ?

Et comme toujours, il existe plusieurs approches pour (tenter de) répondre à cette question, l’approche bayésienne (voir ci-dessous) et la distribution du maximum de scores (prochainement dans un autre billet).

On a vu que le score S d’un alignement était un log de ratio de vraisemblances:

S = log( \frac{P(x,y/M)}{P(x,y/R)} )

Mais en fait, au lieu de P(x,y/M), la probabilité d’obtenir les séquences x et y sachant qu’elles ont un ancêtre commun (le match modèle M), ce que l’on veut vraiment, c’est calculer la probabilité que les deux séquences aient un ancêtre commun connaissant ces deux séquences: P(M/x,y). Grâce à la formule de Bayes, on peut calculer cette probabilité:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y)}

Cependant, au préalable, il faut spécifier les prior: P(M) la probabilité a priori que les deux séquences sont reliées (par un ancêtre commun) et P(R) la probabilité a priori du modèle aléatoire:  P(R) = 1 - P(M).

Ainsi:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y/M)P(M)+P(x,y/R)P(R)}

P(M/x,y) = \frac{P(x,y/M)P(M) / P(x,y/R)P(R)}{1+P(x,y/M)P(M) / P(x,y/R)P(R)}

Finalement, on pose:

S' = S + log( \frac{P(M)}{P(R)} )

On obtient alors:

P(M/x,y) = \sigma(S') avec \sigma(x)=\frac{\exp^x}{1+\exp^x}

\sigma est la fonction logistique. C’est une sigmoide tendant vers 1 quand x tend vers \infty, et vers 0 quand x tend vers -\infty. Cette fonction est bien utile pour convertir des scores faits de sommes en probabilités.

Une fois qu’on a fait tout ça, on peut comparer la valeur finale de P(M/x,y) à 0 pour savoir si nos séquences sont reliées entre elles ou non. Sans rentrer dans les détails, il faut être bien sûr de travailler avec des probabilités (leur somme sur toutes les paires possibles de séquences doit faire 1), ce qui n’est pas toujours le cas dans un scoring scheme construit de façon ad hoc. De plus, les priors prennent beaucoup d’importance quand on cherche un alignement significatif entre une séquence inconnue et une banque de séquences connues. Dans ce cas-là, il faut prendre en compte la taille de la banque puisque la probabilité d’obtenir, par chance, un bon alignement entre deux séquences non reliées entre elles augmente avec la taille de la banque.


Le système de score d’un alignement

20 avril 2009

Supposons que l’on ait deux séquences qui nous intéressent (d’ADN ou de protéines). Nous pouvons aligner globalement ces deux séquences l’une avec l’autre (voir ce billet). Quand on compare comme cela des séquences entre elles, on cherche à savoir si elles ont divergé à partir d’un ancêtre commun sous l’action de processus évolutifs de mutation et sélection.

Les processus mutationnels de base sont les substitutions qui remplacent un résidu par un autre, et les insertion-délétions (indels) qui ajoutent ou enlèvent des résidus (les résidus sont les bases azotées pour l’ADN et les acides aminés pour les protéines). La sélection naturelle joue sur ces variations aléatoires et les crible de telle sorte que certaines modifications vont être plus visibles que d’autres: si une modification est très délétère pour l’individu qui la porte, elle sera contre-sélectionnée, n’apparaîtra pas ou peu dans la génération suivante et disparaîtra (schématiquement…).

Jusqu’à maintenant, pour construire un alignement nous avons utilisé des scores comme par exemple +1 pour le match (deux résidus alignés identiques), 0 pour le mismatch (deux résidus alignés différents) et -1 pour le gap (un résidu aligné avec un gap, dû à un indel). Mais comme nous voulons un système de score qui donne le score le plus élevé à l’alignement le plus plausible biologiquement, il nous faut prendre en compte beaucoup de choses comme l’histoire évolutive de ces séquences, leur structure tridimensionnelle…

Pour cela, nous utilisons un modèle probabiliste, c’est un moyen de simuler l’objet considéré, qui donne différents scénarios possibles avec une probabilité attachée à chacun d’eux. Dans notre cas, les séquences biologiques sont des chaînes de caractères (string en anglais) tirés d’un alphabet de résidus de taille finie (A,T,G,C pour l’ADN). Supposons que le résidu x apparaisse aléatoirement avec la probabilité q_x, indépendamment des autres résidus présents dans la séquence. Si la séquence est notée x_1...x_n, la probabilité d’observer cette séquence est alors le produit q_{x_1}q_{x_2}...q_{x_n}. Dans la suite nous appellerons ce modèle le « modèle de séquence aléatoire », c’est notre hypothèse nulle, contre lequel on comparera les autres modèles.

Concernant notre alignement des deux séquences, le score total correspond à la somme des termes pour chaque paire de résidus alignés plus un terme pour chaque gap. Intuitivement, si nos deux séquences ont un ancêtre commun, on s’attend à ce que les paires de résidus identiques (les identités) soient plus probable que dans le modèle aléatoire et donc contribuent au score par un terme positif, et inversement pour les paires de résidus différents et les gaps. Dans l’interprétation probabiliste, on calcule le logarithme d’un ratio de vraisemblance: la vraisemblance pour les séquences d’avoir un ancêtre commun par rapport à la vraisemblance de ne pas en avoir.

Quelques notations tout d’abord. Considérons deux séquences, x et y, de longueur respective m et n. Le i-ème résidu de x est noté x_i. Les résidus sont tirés d’un alphabet dont les symboles sont représentés par des lettres en minuscule comme a et b.

Vu précédemment, le modèle aléatoire noté R suppose que la probabilité d’avoir les deux séquences correspond au produit des probabilités d’avoir chaque séquence. Formellement, cela donne:

P(x,y/R) = \prod_{i} q_{x_i} q_{y_i}

Dans le modèle alternatif (supposant que les deux séquences ont un ancêtre commun), une paire de résidus alignés est observée avec une probabilité jointe p_{ab}. On peut voir ça comme la probabilité que les résidus a et b dérivent indépendamment tous les deux d’un résidu originel inconnu c qui est leur ancêtre commun. Ainsi:

P(x,y/M) = \prod_{i} p_{x_i y_i}

Les probabilités P(x,y/R) et P(x,y/M) sont des vraisemblances. En probabilité, la vraisemblance est la probabilité d’obtenir les données D sachant le modèle M: P(D/M). Le ratio de vraisemblances est connu sous le nom d’odds ratio:

\frac{P(x,y/M)}{P(x,y/R)} = \prod_{i} \frac{p_{x_i y_i}}{q_{x_i} q_{y_i}}

Pour obtenir un système de score additif, on prend le logarithme de ce ratio (le logarithme d’un produit étant une somme de logarithmes):

S = \sum_{i}{} s(x_i,y_i)

avec: s(x_i,y_i) = \log ( \frac{p_{x_i y_i}}{q_{x_i}q_{y_i}} )

Ici, utiliser un schéma de score additif revient à supposer que les mutations apparaissent indépendamment à chaque site de l’alignement. En probabilité, observer deux événement indépendants revient à multiplier leur probabilité. Et comme le logarithme d’un produit est égal à la somme des logarithmes, on obtient bien notre score additif. Finalement, les scores s(a,b) peuvent être stockés dans une matrice, de taille 4×4 pour l’ADN.

Il est important de voir que lorsqu’un biologiste construit une matrice de substitution ad hoc pour son alignement, il suppose implicitement les probabilités p(a,b) d’obtenir la paire de résidus ab dans son alignement.


L’alignement global deux-à-deux

29 mars 2009

Le titre du billet ne vous évoque probablement rien, mais peut-être en sera-t-il autrement après avoir lu la suite, enfin j’espère parce que l’article est long et qu’il introduit des concepts majeurs, indispensables à la compréhension des méthodes utilisées dans les recherches portant sur l’évolution des génomes !

Comme vu précédemment (ici), l’alignement de séquences biologiques est une tâche fondamentale à l’heure de la génomique. En effet, à notre époque, on peut séquencer le génome d’organismes vivants « à haut débit », de manière quasi-industrielle (le « quasi » est même superflu dans certains cas). Mais obtenir toutes ces séquences n’a pas pour seul but de les stocker dans un ordinateur bien sûr: on séquence de l’ADN dans le but de répondre à des questions scientifiques. Supposons que l’on ait deux séquences d’ADN à notre disposition, ces fragments d’ADN codant chacun pour une protéine permettant de fixer l’oxygène. Or ces deux fragments sont légèrement différents: l’un code pour l’hémoglobine (dans le sang) et l’autre pour la myoglobine (dans les muscles). Comment fait-on pour « voir » ces différences ? Où sont-elles le long de la séquence ? Ces différences au niveau de l’ADN ont-elles des conséquences au niveau de la séquence d’acides aminés ? Ces différences sont-elles importantes pour la fonction de ces protéines ? …

On se rend bien compte qu’il nous faut trouver un moyen pour comparer deux séquences d’ADN. On peut commencer par dire qu’il y a x% de A dans la première séquence alors qu’il y en a seulement y% dans la deuxième, mais ça ne nous avancera pas à grand chose puisque c’est l’enchaînement des nucléotides qui est important. De plus, on sait qu’au cours de l’évolution des mutations apparaissent le long du génome, par exemple un A remplacé par erreur par un T. Il faut donc essayer de repérer le long de la séquence où on eu lieu ces mutations. Dans le cas A->T, on parle de substitutions, c’est une mutation ponctuelle, par opposition aux insertions-délétions (appelés « indels« ) de plusieurs nucléotides comme AATGTCA->AACA où TGT a disparu.

Et bien figurez-vous que ce type de question, les informaticiens se les étaient déjà posées, avant même que l’on sache que l’ADN était le support de l’information génétique. En 1950, Hamming proposait une distance pour quantifier les substitutions entre deux chaines de symboles de même longueur: d_H=2 entre AATGC et ATTGA (A_2->T_2 et C_5->A_5). En 1965, Levenshtein étendait la mesure aux cas de séquences de longueur différente pour prendre en compte les insertions-délétions. Le problème revient alors à trouver le nombre de changements qu’il faut faire subir à une séquence pour la transformer en l’autre, par exemple AATGC -> ATTAGCA. Ainsi, on cherche à calculer la distance d’édition, le nombre d’opérations parmi « substitution », « insertion » et « délétion », qui fait passer d’une chaine de caractères (string en anglais) à une autre. En informatique, ce problème revient à trouver la plus longue séquence en commun entre deux chaines de caractères (LCS pour longest common substring).

Le problème LCS a deux caractéristiques majeures. Il a une sous-structure optimale: le problème peut être décomposé en sous-problèmes plus simples, qui eux-mêmes peuvent être décomposés en sous-sous-problèmes également plus simples, et ainsi de suite jusqu’à ce que la solution soit triviale à trouver. De plus, il a ce qu’on appelle des sous-problèmes chevauchants: la solution d’un sous-problème dépend des solutions aux sous-problèmes de niveau inférieur. De tels problèmes sont résolus grâce à une technique appelée la « programmation dynamique« : une solution optimale du problème peut être construite efficacement à partir de solutions optimales à des sous-problèmes. Pour cela on a besoin d’une procédure appelée « mémoisation » (et non pas « mémorisation »!): sauveguarder dans une table les solutions à un niveau de sous-problèmes de telle sorte qu’elles soient disponibles pour le niveau suivant de sous-problèmes. La programmation dynamique a été inventée par Richard Bellman dans les années 1940-1950. Cette méthode est plus rapide que les méthodes naïves de comparaison de deux séquences, et elle est exacte: contrairement aux méthodes heuristiques qui font des approximations pour aller plus vite, la programmation dynamique trouve la solution optimale du problème et rapidement. Et pour ces raisons, ça déboîte !

Levenshtein, quand il a généralisé la distance de Hamming, a également proposé un algorithme de programmation dynamique pour la calculer. Enfin, on y arrive, c’est en 1970 que Needleman et Wunsch vont comparer deux séquences de protéines à l’aide des méthodes décrites ci-dessus (ça marche pareil pour des séquences d’ADN). Et le reste du billet va justement montrer comment on fait, en pratique.

On suppose que les deux séquences sont homologues, c’est-à-dire qu’elles ont évoluées à partir d’un ancêtre commun, et que les différences entre elles ne sont dues qu’à deux types d’événements, des indels et des substitutions. Pour les aligner l’une avec l’autre, on applique le principe de parcimonie: l’alignement nécessitant le moins d’événements est dit « optimal », c’est la solution. Comme les deux types d’événements arrivent à des fréquences différentes, on utilise une parcimonie pondérée. Par exemple, à une certaine position, un match (deux lettres identiques) vaut +1, un mismatch (deux lettres différentes) vaut 0, et un gap vaut -1: c’est le scoring scheme. Needleman et Wunsch ont donc cherché l’alignement qui maximisait un score (la similarité entre les deux séquences), celui-ci étant la somme des matches, mismatches et pénalités de gap le long de l’alignement.

Quelques notations: on cherche à aligner la séquence s_1 de longueur l_1 avec la séquence s_2  de longueur l_2. Le i-ème élément de s_1 est noté s_1[i] et s_1[3:7] correspond à la sous-séquence de s_1 allant de l’élément 3 à l’élément 7. Le score est ici la similarité entre les deux séquences, notée S_{ij} pour l’alignement de s_1[1:i] avec s_2[1:j].

Comment marche la programmation dynamique? On commence par créer une matrice F dont l’entrée F_{i,j} correspond au score  de l’alignement entre s_1[1:i] et s_2[1:j]. Ensuite, on remplit la matrice récursivement en commençant en haut à gauche: remplir F(i,j) ne dépend que de 3 cases remplies précédemment F(i-1,j-1), F(i,j-1) et F(i-1,j). Bien sûr, il y a des conditions au bornes: dans la première ligne et la première colonne, on ne met que des gaps.

matrixf

Pour calculer F(i,j), on commence par  déterminer le score de s_1[i] aligné avec s_2[j] (match ou mismatch), puis on applique la règle de décision suivante:

equationfillf

Une fois que la matrice est remplie, on fait ce qu’on appelle le traceback. On part d’en bas à droite de la matrice, case avec le score le plus élevé (similarité de l’alignement optimal), et on remonte en suivant le meilleur chemin: ce chemin dans la matrice correspond à l’alignement optimal. En fin de compte, on dispose bien d’une méthode pour comparer quantitativement deux séquences d’ADN.

Prenons un exemple simple: AAT et ACAT.

1. on définit le scoring scheme:

match=2 ; mismatch=0 ; gap=-1

2. on remplit la première ligne et la première colonne avec des gaps:

fillf_row1col1

3. on remplit le reste de la matrice avec la formule donnée précédemment:

fillf1

4. on fait le traceback:

traceback

5. on montre l’alignement obtenu, avec un score de 5:

alignJe précise que cette méthode marche aussi pour les protéines, simplement, au lieu que l’alphabet soit de 4 lettres {A;T;G;C}, il est de 20 lettres, une par acides aminés. De plus, la comparaison de deux séquences peut avoir plusieurs alignements optimaux. Et pour finir d’expliquer le titre, on parle ici d’alignement global parce qu’on essaie d’aligner les deux séquences sur toute leur longueur. Vous vous en doutez, il existe aussi l’alignement local…

Maintenant qu’on dispose d’un algorithme pour comparer deux séquences, on peut s’amuser un peu à changer les paramètres (oui, on s’amuse de peu en ce bas monde). Avec les paramètres match=2, mismatch=0 et gap=-1, on obtient l’alignement suivant, de score 67:

aligngap-1

En passant la pénalité de gap à -2, on obtient un alignement tout différent, de score 58:

aligngap-2

Et en mettant le gap à -3, on obtient encore un autre alignement, de score 52 celui-là:

aligngap-3

Alors, lequel est le meilleur ? Chacun est optimal, mais cette optimalité est fonction du scoring scheme choisi, donc garder celui de score le plus élevé n’est pas forcément le plus judicieux. On peut aussi garder celui qui est le plus plaisant à l’oeil (vous rigolez mais au début les chercheurs faisaient comme ça !). Une autre possibilité est se fixer un scoring scheme qui semble relativement correct et de garder le même tout le temps. Enfin bon, tout cela n’est pas très satisfaisant, vous en conviendrez. Le vrai problème derrière tout ça consiste à écrire un modèle d’évolution des séquences et de l’utiliser lors de l’alignement: on garderait dans ce cas l’alignement ayant la vraisemblance la plus élevée (au sens statistique du terme: l’alignement qui maximise la vraisemblance d’observer ces séquences sous un modèle d’évolution donné). Mais ce sera pour d’autres billets…


Invasion d’ETs

8 mars 2009

Chaque cellule vivante possède un génome, la molécule d’ADN, qui contient l’information génétique nécessaire à sa survie (voir ce billet). On a vu précédemment que, selon les organismes, la taille du génome pouvait beaucoup varier (voir ce billet), et que ces variations entre génomes sont principalement dues aux éléments transposables, les ETs (voir ce billet).

En biologie, lorsque l’on s’intéresse à quelque chose sous l’angle de l’évolution, se posent généralement deux questions: l’origine (comment c’est apparu) et la maintenance (pourquoi ça s’est maintenu au cours du temps). Prenons l’exemple d’une famille d’ETs dans un génome. On parle de « famille » parce qu’un ET est capable de transposer (de changer de place) au sein du génome-hôte et, par extension, de se multiplier: on dit alors que toutes les copies ayant pour ancêtre commun cet ET-là forment une famille. Mais à l’origine, cet ET, il vient d’où ?

Mettons de côté la question de l’origine du tout premier ET (que j’ai brièvement abordée ici) et considérons simplement un organisme, par exemple le poisson rouge qui tourne dans son bocal. Si l’on séquence son génome, on y trouvera des ETs. Ils viennent  soit de son père et/ou de sa mère, comme tout autre constituant de son génome (transmission verticale), soit d’ailleurs (transmission horizontale), et bien sûr, c’est plus amusant de s’intéresser au deuxième cas… Mais avant de savoir si les ETs de notre poisson rouge lui viennent de l’espace, on peut déjà se demander si ça arrive fréquemment que des ETs soient transmis horizontalement.

En 2008, des chercheurs ont répondu à cette question dans le cas des animaux en montrant qu’une famille d’ETs (appelée SPIN pour SPace INvaders…) avait envahie par transferts horizontaux plusieurs génomes de mammifères au cours de centaines de millions d’années (voir cet article). Mais au juste, comment fait-on ça ?

Depuis quelques années, on a les moyens techniques pour séquencer des génomes appartenant à beaucoup d’espèces différentes. Dans notre cas, les auteurs de l’étude recherchaient les ETs présents dans le génome d’un lémurien, Otolemur garnettii (la jolie petite bête photographiée ci-dessous). Ils ont utilisé une approche bioinformatique d’alignement de séquences et ont trouvé un ET qu’ils ne connaissaient pas très bien, appartenant a priori à la famille des hAT. Pour en savoir plus, ils ont regardé si cet ET n’était pas présent chez d’autres espèces comme l’homme, la souris, l’éléphant, le chien, la chauve-souris… et ils l’ont trouvé chez certains d’entre eux mais pas tous ! Ça met la puce à l’oreille, vous en conviendrez… En effet, la phylogénie (l’arbre généalogique) ci-dessous montre que les SPINs sont présents dans 6 génomes de tétrapodes mais pas dans les autres (les barres verticales indiquent la distribution des copies en fonction de leur âge).

transferts horizontaux d'ETS chez les mammifères

On peut imaginer tout d’abord que l’ancêtre de tous les génomes analysés (à la racine de l’arbre, donc il a vécu il y a environ 350 millions d’années) possédait déjà cet ET et qu’au cours du temps, de nouvelles espèces sont apparues, certaines perdant cet ET pour diverses raisons, tandis que d’autres le gardaient. Dans ce cas-là, les copies de l’ET devraient être très anciennes et toutes avoir à peu près le même âge. L’autre possibilité est d’imaginer que plusieurs transferts horizontaux sont arrivés dans les différentes branches de l’arbre et donc que les ETs sont beaucoup plus jeunes et que les copies d’un génome peuvent avoir un âge différent de celles dans un autre génome (c’est-à-dire le long d’une autre branche de l’arbre). Tenez vous bien: on trouve justement que les copies d’ETs sont jeunes et que, par exemple, les copies chez la chauve-souris sont plus jeunes que celles chez le rat !

Tout ça veut donc dire qu’il y a eu des transferts horizontaux, et ce plusieurs fois, au cours de l’évolution menant à ces espèces. On ne sait pas très bien comment de telles choses arrivent, certains supposent qu’un parasite d’une espèce peut en parasiter une autre et qu’il peut faire la navette de l’une à l’autre en transférant du matériel génétique de temps en temps, ce qui pourait être le cas des poux ou bien des achariens, mais rien n’a encore été observé. Comme quoi, les morceaux d’ADN aussi ça aime se balader… !


Alignement de séquences

11 février 2009

Maintenant que l’on a vu ce qu’est une séquence d’ADN ainsi que le lien entre darwinisme et mendélisme, et après avoir introduit la génétique des populations par l’équilibre de Hardy-Weinberg, il est temps de parler de bioinformatique. Depuis que l’on sait que l’ADN est le support de l’hérédité, on s’est dépêché de trouver le moyen de séquencer des fragments d’ADN. On s’est donc retrouvé avec, d’un côté des séquences d’ADN, et de l’autre des questions du type: quelle est la fonction de cette séquence ? ces deux séquences ont-elles un ancêtre commun ?  à quoi ressemble l’action de la sélection naturelle sur une séquence d’ADN ? quels types de séquences sont trouvés dans un génome ? …

Ces questions ont ouvert de vastes champs de recherche tels que l’évolution moléculaire, la génomique fonctionnelle et comparative, la bioinformatique… Notamment, on s’est vite rendu compte que l’alignement de séquences était un point d’entrée quasi obligé. Par exemple, on se retrouve avec le problème suivant: étant donné une séquence X de longueur m et une séquence Y de longueur n, quel est le meilleur alignement de X avec Y ? Et par « meilleur alignement », on entend « celui qui reflète le mieux l’évolution des deux séquences à partir de leur ancêtre commun ».

Le cas le plus simple correspond à l’alignement deux-à-deux (pairwise en anglais). Mais on peut aussi chercher à comparer plus que deux séquences, faisant ainsi un alignement multiple. Par ailleurs, la notion de « meilleur alignement » requiert de connaître le processus évolutif auquel les séquences ont été soumises, on a donc besoin de modéliser les mutations au cours du temps. Cela permet ensuite de tracer la phylogénie d’un ensemble de séquences, cad leur généalogie. En fin de compte, les prédictions de la génétique des populations, et donc de la théorie de l’évolution, vont pouvoir être testées, et elles sont maintenant amplement prouvées.


%d blogueurs aiment cette page :