Évitez la navigation.
AFUL · Parinux · FFII France · APRIL · ADULLACT · MongueursPerl · Wikipédia · OFSET · Scidéralle · LéaLinux · LinuxFrench · LinuxFr · FirstJeudi · AgendaLibre

sylvain(à)lhullier(.)org
À force d'aller au fond des choses, on y reste.
Jean Cocteau (J'en suis resté ...)

Description, recherche et prédiction de structures d'ARN

... ou deux ans de thèse de doctorat

 Institut Pasteur  UMLV

Durant les années 1999-2000, j'étais en thèse de doctorat en informatique fondamentale (en algorithmique génomique) à l'Université de Marne-la-Vallée. Le sujet de ma thèse était : Description, recherche et prédiction de structures d'ARN.

Mon directeur de thèse était Maxime Crochemore, directeur du Laboratoire d'Informatique. La personne qui m'encadrait était Marie-France Sagot, à l'époque chargée de recherches en informatique à l'Institut Pasteur.


Voici mon rapport de DEA, c'est une bonne introduction à la problématique de ma thèse et aux travaux réalisés durant ces années.

1. Motivation biologique

1.1. Introduction

Le patrimoine génétique de tout être vivant est contenu dans une ou plusieurs molécules d'ADN constitutives de son génome. Le fait que l'ADN contient toutes les informations concernant les fonctions vitales d'un organisme est un des fondements de la biologie moléculaire.

Les gènes sont les parties de la molécule d'ADN qui contiennent la majeure partie de l'information génétique. La grande majorité des gènes contient des instructions pour la synthèse des protéines ; les protéines sont des macromolécules qui jouent un rôle essentiel dans la vie de l'organisme.

Le modèle général de décodage du message génétique est constitué de deux mécanismes fondamentaux :

La figure suivante résume l'ensemble du mécanisme :

Schéma de principe du décodage du message génétique
Schéma de principe du décodage du message génétique

Comprendre le message génétique contenu dans l'ADN, et en particulier identifier les différents gènes à l'intérieur, est un des plus grands défis de la biologie moléculaire actuelle.

Ces dernières années, la vitesse de séquençage de l'ADN s'est accrue et un séquençage massif des séquences nucléiques s'est opéré, augmentant très fortement la taille des bases de données. Cette tendance ne va sans doute pas changer avant de nombreuses années.

Pour ces raisons, le décodage de l'information génétique et sa compréhension ne peuvent plus s'en remettre aux seules méthodes expérimentales et doivent faire appel à l'utilisation d'algorithmes de prédiction. Avec des moyens informatiques, on peut, entre autres, tenter de prédire les gènes, les séquences régulatrices, les structures d'ARN à partir de leur séquence ...

1.2. L'ARN

ADN (acide désoxyribonucléique) et ARN (acide ribonucléique) sont des macromolécules, plus précisément de longs polymères linéaires, formés de quatre types de monomères (petites molécules) appelés nucléotides. Un nucléotide est formé d'une base, d'un sucre et d'un ou plusieurs groupements phosphate.

Pour l'ARN, il y a quatre bases : l'adénine (A), la guanine (G), la cytosine (C) et l'uracile (U). Les quatre nucléotides déterminés par ces bases, sont attachés à la suite les uns des autres, un peu comme des perles enfilées sur un collier.

L'une des propriétés chimiques fondamentales des nucléotides est leur capacité à s'apparier en formant des liaisons hydrogène. L'adénine s'apparie avec l'uracile et la guanine s'apparie avec la cytosine (appariements Watson-Crick) ; de plus, il peut y avoir des appariements moins stables tels que l'appariement guanine-uracile.

Contrairement à l'ADN qui est composé de deux chaînes de nucléotides complémentaires reliées entre elles par des appariements, l'ARN est une molécule simple brin. Cette structure en simple brin n'empêche pas les appariements ; en effet les nucléotides d'un même ARN s'apparient souvent entre eux, certaines parties complémentaires de la chaîne de l'ARN se plaçant les unes en face des autres ; voyez l'exemple de la figure suivante :

Exemple d'appariements dans l'ARN
Exemple d'appariements dans l'ARN

Les différents éléments de structure que forment les ARNs par leurs repliements seront détaillés dans la section suivante.

2. Problème général : description de structures d'ARN

L'axe principal du travail entrepris est la description des structures d'ARN et leur recherche dans les bases de données génomiques.

Pour permettre une description simple et efficace de telles structures à des fins de recherche, nous avons identifié un certain nombre d'éléments de base constitutifs des structures d'ARN ainsi que leur caractéristiques propres et les relations qu'ils peuvent avoir entre eux.

2.1. Éléments structuraux de base

Nous avons dégagé les éléments structuraux suivants :

La figure suivante présente ces différents élément structuraux :

Éléments de structure
Éléments de structure

2.2. Propriétés des éléments structuraux

Ces éléments structuraux ont des propriétés qui permettent de les caractériser lors d'une recherche par exemple.

On peut vouloir donner à un élément une taille minimale et une taille maximale. Par exemple, on peut chercher une boucle dont la longueur de la séquence est comprise entre 3 et 5 bases.

À certaines positions, on peut vouloir retrouver les mêmes bases. Par exemple, sur une tige on veut trouver la base A à une position précise. On parle alors de base conservée.

Dans le cas où l'on cherche un élément comportant plusieurs bases conservées, on peut définir un taux d'erreur sur ces bases ; par exemple on peut dire qu'il faut qu'au moins 75% des bases soient conservées pour que la structure soit retenue.

On peut aussi établir un taux d'erreur en ce qui concerne les appariements dans une tige. Par exemple, sur une tige de longueur de 4 à 6, on peut vouloir tolérer 1 paire de bases qui ne seraient pas appariées.

Certaines séquences génomiques sont connues pour leur composition en bases A/C/G/U ; par exemple, 60% de A-U et 40% de C-G. Il est donc envisageable de vouloir sélectionner une structure en fonction de la composition en bases de la séquence qui la forme.

Enfin, on peut choisir de chercher un élément structural en fonction de son énergie libre.

2.3. Relations entre éléments structuraux

Pour rechercher une structure complète, il faut définir dans quel ordre et selon quel schéma s'arrangent les éléments structuraux qui la composent. La description d'une structure d'ARN passe donc par celle de ses éléments mais aussi par celle de la manière dont sont disposés les éléments.

On peut par exemple chercher une structure contenant une boucle multiple comportant (en telle et telle position) deux tiges fermées par une boucle, dont une des deux tiges comporte une excroissance.

Tous les objets et les propriétés définis précédemment ne concernent que la structure secondaire de l'ARN ; il est aussi possible de décrire certaines interactions tertiaires. On peut en effet envisager que certaines bases appartenant à des boucles (ou toutes autres bases non appariées) puissent s'apparier entre elles, par exemple pour former un pseudo-noeud.

Un pseudo-noeud est un élément de structure correspondant à deux sites qui s'apparient entre eux, entrelacés le long de la séquence primaire de l'ARN avec deux autres sites qui s'apparient entre eux, voir figure suivante :

Structure d'un pseudo-noeud
Structure d'un pseudo-noeud

2.4. Propriétés globales de la structure

Comme pour chaque élément structural, il est possible d'établir des propriétés globalement définies pour l'ensemble de la structure. Par exemple les longueurs minimales et maximales, le taux d'erreur global sur les bases conservées, le taux d'erreur global sur les appariements, la composition en bases A/C/G/U, l'énergie libre.

De plus il est possible de définir d'autres conditions sur la structure.

On peut parler d'un taux de présence des éléments structuraux. Cela permet par exemple de retenir une structure à laquelle il ne manque que peu d'éléments structuraux par rapport à la description mais qui peut tout de même nous intéresser.

On peut aussi envisager d'établir des propriétés conditionnelles de toute nature. Par exemple : « si l'élément 1 est présent, alors l'élément 2 n'est pas nécessaire ».

2.5. Classification des caractéristiques

Dans l'ensemble de ces propriétés, on peut distinguer deux groupes selon si les propriétés en question nous aident dans la recherche de la structure, ou si elles servent uniquement lors de la phase finale de vérification.

L'information sans doute la plus déterminante dans la phase de recherche est l'agencement des hélices. C'est pour cela que les critères de distance entre éléments, de même que les taux d'erreur sur les appariements sont utiles à la recherche. Une autre information importante est représentée par les bases conservées et donc par le taux d'erreur autorisé sur celles-ci.

Des critères tels que la composition en A/C/G/U, l'énergie libre ou les propriétés conditionnelles n'entreront sans doute en compte que dans la phase de vérification.

2.6. Bases dégénérées

Lorsque nous parlons de séquence conservée dans les éléments structuraux, nous voulons autoriser le fait qu'à une position donnée, plusieurs bases sont envisageables. On peut par exemple vouloir autoriser à une position soit une adénine soit une guanine : on dit alors qu'à cette position, on est en présence d'une base dégénérée.

Nous avons utilisé comme notation pour les bases dégénérées, le code IUPAC résumé dans la figure suivante :

Origine du choix des lettres IUPAC
Code Adénine Cytosine Guanine Uracile Origine du choix de la lettre
AA   Adenine
B CGUpas A, B suit A dans l'alphabet
C C  Cytosine
DA GUpas C, D suit C dans l'alphabet
G  G Guanine
HAC Upas G, H suit G dans l'alphabet
K  GUKeto
MAC  aMino
NACGUaNy (quelconque en anglais)
RA G puRine
S CG interaction forte (3 ponts H) Strong en anglais
T   UThymine
U   UUracile
VACG pas U, V suit U dans l'alphabet
WA  Uinteraction faible (2 ponts H) Weak en anglais
XACGUcomme une inconnue en mathématiques
Y C UpYrimidine

Ce tableau se lit ainsi : un R correspond à une adénine ou à une guanine.

Si la présence de bases dégénérées dans les motifs que nous recherchons a une sens clair (cela correspond à une séquence consensus : plusieurs bases sont autorisées à cette position), de telles bases dans les bases de données correspondent à des indécisions de lecture lors du séquençage. Nous avons entamé une réflexion sur la manière de les gérer lorsqu'elles apparaissent ainsi dans les bases de données.

3. Problème particulier

Au début de ce travail, le problème entier n'a pas été entamé de front, nous nous sommes concentrés sur un problème particulier.

En effet, nous nous sommes plus particulièrement intéressés au problème de décrire et de trouver dans les bases de données un objet de type tige-boucle avec éventuellement une ou des excroissances et/ou une ou des boucles internes.

3.1. Motivation biologique

Cette structure simple en tige-boucle trouve son intérêt notamment dans les introns de groupe II, objet qui intéresse beaucoup les biologistes. Voir la figure suivante :

Structure secondaire d'un intron de groupe II
Structure secondaire d'un intron de groupe II
(cliquez sur l'image pour zoomer)

Rechercher cette structure en entier pose plusieurs problèmes, dont un des principaux est dû à la taille de la base de données dans laquelle on pourrait souhaiter effectuer une telle recherche ; EMBL comportait 4640872 séquences pour 3456074245 bases le 10/09/1999. La recherche pourrait ainsi être trop complexe. L'idée est alors de rechercher une sous-partie caractéristique pour ensuite vérifier si le reste de la structure correspond à ce que l'on recherche. La recherche de cet élément caractéristique serait alors optimisé ; cela permettrait d'obtenir un temps de recherche raisonnable.

Un motif en tige-boucle correspondant au domaine V des introns de groupe II se trouve décrit dans la littérature comme étant caractéristique des introns de ce groupe.

Voir la figure suivante :

Tige-boucle caractéristique des introns de groupe II
Tige-boucle caractéristique des introns de groupe II

Ce motif est assez bien conservé à la fois en termes de structure et de séquence. Il va nous permettre de tester nos méthodes et nos algorithmes sur un cas réel et pertinent.

Nous ignorons à quel point cet élément de structure est effectivement caractéristique des introns de groupe II. Nous ne savons pas si d'autres objets biologiques comportent cette même structure, ni à partir de quel niveau de ressemblance à cet objet on peut espérer sélectionner tous les introns de groupe II de la base.

3.2. Représentations interne et externe

Pour utiliser cet objet biologique dans un algorithme, nous avons dû le représenter informatiquement de deux façons :

L'objet manipulé est de deux types : séquentiel et structural. Il nous faut manipuler ces deux informations. Pour les algorithmes implémentés, nous avons opté pour l'instant pour des représentations simples.

Une réflexion plus profonde a été entamée pour imaginer des représentations qui conviendraient pour des structures complexes, mais qui restent simples pour l'utilisateur (externe) et qui permettent d'être le plus efficace possible algorithmiquement (interne).

4. Travaux

4.1. Algorithme de Wu et Manber

L'idée principale de leur méthode est l'utilisation d'un champ de bits, chaque bit signifiant qu'une lettre du motif correspond à une lettre du texte.

Chaque lettre pouvant intervenir dans le texte est représentée par un mot machine de la longueur du motif, dont les bits indiquent si la lettre en question fait partie du motif à la position correspondante du bit. Il faut ensuite confronter ces champs de bit représentant le motif au texte : pour chaque lettre du texte, le bit courant est calculé grâce à des opérateurs logiques bit à bit, et par décalage du mot en mémoire d'un bit vers la droite, le programme passe à la lettre suivante du texte. La présence d'un bit à 1 à une certaine position permet de savoir si le motif est présent dans le texte.

Il est possible d'adapter cet algorithme afin de permettre les erreurs telles que des insertions, des suppressions et des substitutions, choses qui nous intéressent fortement.

Il a fallu adapter cet algorithme, initialement destiné à effectuer la recherche de motifs dans un texte classique, à un aspect propre à notre problème : les bases dégénérées. Dans leur article, Wu et Manber considèrent que seule la lettre A peut correspondre à la lettre A ; ce n'est pas le cas pour nous (voir section "Bases dégénérées").

Le principal avantage de cette méthode est sa rapidité ; du fait que l'on ne manipule principalement que des champs de bit avec des opérateurs bit à bit, les calculs sont très efficaces car les instructions effectuées par le programme sont très proches du processeur.

Par contre, les inconvénients sont nombreux.

Tout d'abord, la taille des motifs est a priori limitée à la taille d'un mot machine, même s'il est envisageable d'en gérer plusieurs pour les longs motifs. Dans ce cas il faudrait alourdir le programme pour la gestion des ``retenues'' et ce dernier perdrait en efficacité.

Ensuite, cette méthode est très peu souple, du fait que l'on ne peut pas affecter de note à un ``match''. La seule chose que l'on peut dire c'est que telle lettre du texte correspond ou non à telle lettre du motif ; il ne nous est pas permis d'affecter un score qui entrerait dans le calcul d'un score pour la séquence candidate.

Enfin, on ne tient pas compte, lors de cette recherche, de la structure secondaire du motif : seule la séquence primaire est utilisée comme information. Il est évident que nous ne pouvons pas nous contenter de chercher un motif uni-dimensionnel.

Nous avons appliqué notre programme à la banque de données EMBL et avons obtenu un certain nombre de séquences. Nous avons consulté les banques annotées pour voir à quoi correspondait certaines de ces séquences ; devant le côté fastidieux de la chose, nous ne l'avons pas fait pour toutes les séquences ; mais pour celles pour lesquelles nous avons cherché à quoi elles correspondaient, nous avons effectivement relevé qu'elles étaient à l'intérieur d'introns de groupe II.

Ce fut le premier élément tendant à confirmer l'idée que le motif en tige-boucle que nous utilisons était effectivement caractéristique de ces introns.

Devant la constatation de l'aspect fastidieux de la recherche dans les banques annotées, nous avons pour projet de mettre au point dans un avenir plus ou moins lointain un programme qui irait analyser automatiquement ces banques et nous renseignerait sur la nature des séquences détectées.

4.2. Programmation dynamique

Devant le manque de souplesse de l'algorithme précédant, nous nous sommes penchés sur la programmation dynamique pour effectuer un alignement global du motif et local des séquences.

Cette méthode demande l'établissement d'un score entre lettres. Le choix de ce score est difficile car il conditionne très fortement les résultats, et son choix est souvent un peu arbitraire. Nous avons utilisé le score suivant, qui est très simple, ou x est une lettre du texte et y une lettre du motif :

s(x,y) =

Puis nous utilisons une matrice rectangulaire de taille (taille du motif+1)*(taille de la séquence+1). Cette taille peut être réduite (en se basant sur une observation toute simple, il est facile de ramener l'occupation mémoire en O(taille du motif)).

La première ligne est initialisée à 0 ; la première colonne est initialisée de façon à prendre en compte les gaps successifs dus à ces positions. Voir figure suivante :

Matrice de programmation dynamique
Matrice de programmation dynamique

La matrice sera ensuite remplie au moyen de la relation suivante :

M(l,c) = max

m[l] est la lettre d'indice l dans le motif et t[c] celle d'indice c dans le texte ; '-' représente un gap.

La première ligne correspond à l'alignement d'une lettre du texte sur une lettre du motif. La seconde ligne correspond à l'insertion d'un gap dans l'appariement, c'est-à-dire qu'une lettre du motif ne sera alignée avec aucune lettre du texte. La troisième ligne correspond à l'insertion d'une lettre du texte dans le motif, c'est-à-dire que cette lettre ne sera alignée avec aucune lettre du motif.

La figure suivante résume les dépendances entre cellules de la matrice :

Dépendances entre cellules de la matrice
Dépendances entre cellules de la matrice

Pour savoir où sont situés les ``matches'' dans le texte, il suffit de considérer la dernière ligne de la matrice de programmation dynamique et d'y lire le score. C'est en effet sur cette ligne que l'on a pris en compte l'ensemble du motif.

Il faut ensuite établir un score seuil à partir duquel on considère avoir trouvé un fragment de séquence ressemblant au motif. Cette tâche est assez difficile car, pour l'instant, aucun autre critère que notre intuition ne nous permet de le choisir. Si on utilise le score présenté en section "Programmation dynamique", un score de 64 correspond à un alignement sans gap, c'est-à-dire que chacune des lettres du motif est mise en correspondance avec une lettre du texte sans insertion de trous.

À titre d'exemple, nous avons travaillé avec des scores proches de 60.

Ensuite, pour retrouver l'alignement à partir de sa position de fin, c'est-à-dire retrouver la suite d'événements qui ont donné lieu à un ``match'' concluant, il faut partir de cette case de score supérieur au seuil et reparcourir à l'envers le chemin dans la matrice ; cette méthode porte le nom anglophone de ``backtraking''.

Cet algorithme a donc été implémenté et plusieurs améliorations ont été réalisées :

Voici les résultats de l'algorithme classique de programmation dynamique, c'est-à-dire sans tester les appariements, sur le motif de la figure "Tige-boucle caractéristique des introns de groupe II" appliqué à une sous-banque d'EMBL. La courbe de la figure suivante présente le nombre de séquence ayant un score d'alignement local donné avec le motif :

Nombre de séquences en fonction du score d'alignement
Nombre de séquences en fonction du score d'alignement

Cette figure se lit ainsi : 57000 séquences ont un score de 59. Voici les valeurs pour les scores supérieurs à 59 : 60 : 5742, 61 : 1545, 62 : 481, 63 : 24 et 64 : 4.

Ces résultats ne nous permettent pas de conclure à quoi que ce soit d'intéressant aujourd'hui ; il est par contre certain que la structure secondaire doit être prise en compte dès la première étape de la recherche.

L'intérêt de cette méthode est sa grande souplesse, le score peut prendre en compte chaque cas d'alignement. Le problème majeur de l'algorithme est sa lenteur d'exécution.

4.3. Programmation dynamique éparse

L'idée principale est de réduire le temps de calcul en ne parcourant pas toute la matrice de programmation dynamique. En effet pour la grande majorité des zones des séquences, la différence avec le motif est si grande qu'il suffirait de calculer quelques cellules de la matrice pour se rendre compte que cette sous-séquence ne nous intéresse pas.

L'idée est alors d'introduire la notion de distance d'édition qui comptabilise le nombre d'opérations atomiques (insertion, suppression et substitution de lettres) nécessaires pour passer du motif recherché au texte. Il faudra aussi établir un seuil au delà duquel nous pouvons éliminer la séquence candidate.

Nous souhaiterions aussi vérifier les appariements de la structure au fur et à mesure de la recherche. En effet certains textes peuvent très bien convenir en ce qui concerne la séquence et être inacceptables en ce qui concerne les appariements : il est très probable, qu'en testant la complémentarité des bases censées s'apparier et en établissant un seuil de tolérance aux non-appariements, que nous pourions réduire l'espace de recherche et gagner ainsi du temps.

Une solution algorithmique existe : la programmation dynamique éparse (en anglais : ``sparse dynamic programming'').

Il s'agit de considérer la liste des positions envisageables du motif dans le texte et pour chaque position, le nombre d'opérations atomiques nécessaires pour y parvenir. Initialement cette liste comporte toutes les positions du texte ; selon que la lettre du texte correspond ou non à la première lettre du motif, la distance sera 1 ou 0. À chaque étape (c'est-à-dire, à chaque nouvelle lettre du motif que l'on considère), la distance de chaque élément de la liste est mise à jour en fonction de la présence ou non de la nouvelle lettre après cette position ; les éléments de la liste qui ont une distance d'édition trop grande sont oubliés. L'algorithme se termine lorsque la liste est vide, auquel cas le motif n'est pas présent dans le texte, ou bien lorsque l'on a lu la dernière lettre du motif et les positions restantes sont les positions de fin des occurrences du motif dans le texte.

Cet algorithme a besoin d'un espace en O(taille du texte) pour stocker la liste et aura une complexité inférieure à O(taille du texte*taille du motif) puisqu'au pire toutes les séquences correspondent au motif jusqu'à la fin.

La méthode que nous venons de décrire ne permet pas encore de prendre en compte la structure secondaire de l'objet de type tige-boucle. Pour cela, il faudrait effectuer la recherche de chacun des éléments des hélices en même temps (chacun par la méthode précédente dans des listes distinctes), et à chaque étape vérifier les appariements.

Cette vérification des appariements consiste surtout à tester si un brin d'une hélice a une séquence candidate pour être le brin qui lui correspond à une distance acceptable (les distances entre éléments sont définies par l'utilisateur, cf section "Problème général").

La figure suivante montre un exemple d'exécution. Nous sommes en présence de la recherche d'une tige-boucle avec boucle interne ; les éléments A et A' ainsi que B et B' doivent s'apparier, C et E forment la boucle interne et D la boucle en épingle.

Schéma de principe de la programmation dynamique éparse
Schéma de principe de la programmation dynamique éparse

Chaque motif A, B, B' et A' possède sa propre liste de positions représentées par des rectangles grisés ; A a, à cet instant, un préfixe en positions 1, 2 et 3 ; les carrés gris clair représentent les positions précédemment testées ; les gris foncés, la position actuelle. A et B sont étendus par la droite ; A' et B' par la gauche.

Par exemple, les éléments 4, 7, 2 et 10 seraient conservés lors du passage à l'étape suivante, si les distances a, b et c sont comprises dans leurs fourchettes respectives. Au contraire, l'élément 6 n'aura jamais d'élément pour s'apparier, il n'est donc pas nécessaire de le conserver dans la liste ; il en est de même pour le 3.

Il nous faut tester les distances entre chaque occurrence de chaque motif dont les distances sont connues. Pour cela nous allons utiliser le fait que les éléments de la liste restent ordonnés pour ne pas vérifier naïvement si chaque élément a son élément correspondant.

La vérification de ces distances va coûter, à chaque étape, un temps proportionnel au nombre d'éléments des listes et à la ``largeur'' des fourchettes des distances acceptables.

5. Conclusion

Cette solution de la programmation dynamique éparse mérite d'être approfondie. Il reste de nombreux questions et problèmes en suspens.

Le premier d'entre eux est l'établissement des fonctions de distance qui, comme à chaque fois pour les scores, est un problème épineux. La première distance concerne la correspondance entre séquences primaires, la seconde les appariements entre brins d'une même tige.

La question des seuils se pose une nouvelle fois ici aussi ; est-il envisageable de n'avoir qu'un seul seuil que la somme des deux distances ne devrait pas dépasser, ou bien aura-t'on deux seuils distincts ? Dans le premier cas, se posera le problème du réglage des deux fonctions de distance l'une par rapport à l'autre ; dans le second, il nous sera difficile de faire des réglages des seuils par essais successifs en raison de la dépendance des résultats à ces deux seuils.

Il faut ne pas oublier enfin de gérer les éléments de la structure qui ne s'apparient pas : les boucles et les excroissances doivent aussi être prises en compte pour l'établissement de la distance d'édition de l'ensemble de la structure.

Ce site respecte les standards de l'internet :
XHTML 1.1   ·   CSS v2   ·   Accessibilité
Plan du site  ·  Signature  ·  Imprimer la page
© 1999-2008 Sylvain Lhullier
http://sylvain.lhullier.org/cv/these.html
Creative Commons Attribution-ShareAlike