Action ATGC

Précédent : Logiciels Remonter : Action ATGC, Action Transversale Génome Suivant : Actions industrielles

Résultats nouveaux

Un algorithme de comparaison par programmation dynamique prenant en compte les motifs.

Participants : Jean-Paul Comet et Jacques Henry

Mots-clés : alignement, motif

Résumé : On montre que l'on peut modifier l'algorithme de Needleman & Wunsch pour favoriser l'alignement de motifs.

Lorsque le biologiste aligne deux séquences, il cherche à mettre en correspondance les motifs qu'il sait être pertinents pour les séquences considérées. La méthode proposée vise alors à prendre en compte des informations autres que la seule structure primaire de la séquences. Typiquement, on pense à des séquences protéiques et à l'utilisation de motifs décrits dans la base PROSITE. Il ne s'agit pas d'imposer l'alignement de ces motifs car dans le cas où ces motifs apparaîtraient dans un ordre différent sur les deux séquences, on aboutirait à une impossibilité. La méthode consiste donc à réaliser un alignement par programmation dynamique lettre à lettre au sens de Needleman et Wunsch, en attribuant un bonus B supplémentaire à tout chemin mettant en correspondance le même motif sur les deux séquences.

La première phase consiste donc à repérer les occurrences des motifs de la base sur les deux séquences à l'aide d'un algorithme de ``pattern matching''.

On suppose ici pour simplifier que le motif apparaît avec la même longueur sur les deux séquences. Si $S_{i,j}$ désigne le score optimal à la position (i,j), le coût de la substitution de par , r le coût d'insertion, on a la récurrence:

$\begin{displaymath}S_{i,j}=max(S_{i-1,j-1}+d(a_i,b_j), S_{i-1,j}+r, S_{i,j-1}+r)\end{displaymath}$

en tout point (i,j) sauf ceux correspondant à la fin d'un motif détecté sur les deux séquences. En un tel point, on a le quadruple choix suivant, en supposant le motif de longueur :

$\begin{displaymath}\begin{array}{ll}S_{i,j}= &max(S_{i-1,j-1}+d(a_i,b_j), S_{... ...-l+1,j-l+1}+\sum_{k=0}^{k=l-1}d(a_{i-k},b_{j-k})+B)\end{array} \end{displaymath}$

Dans le cas d'un motif à longueur variable, la difficulté vient de l'alignement de deux occurrences du même motif. Il ne suffit pas de connaître l'occurrence du motif et sa longueur, mais il faut aussi connaître la structure du motif, pour pouvoir positionner les insertions aux bons endroits.

La méthode a été testée pour un petit nombre de cas et semble fournir des résultats prometteurs. Cette méthode est en cours d'intégration dans le logiciel LASSAP. Cependant, la difficulté principale vient de la détermination des bonus à attribuer à chaque motif. Pour la lever, il faut envisager une étude à grande échelle. Chaque motif se verrait affecter d'une valeur différente prenant en considération, sa longueur mais aussi sa signification biologique.

Amélioration de l'algorithme de Needleman & Wunsch

Participants : Jean-Paul Comet et Jacques Henry

Mots-clés : alignement, algorithme de Needleman & Wunsch

Résumé : On montre qu'une formulation différente de l'algorithme du plongement invariant permet de trouver une forme équivalente à l'algorithme de Needleman & Wunsch, mais nécessitant moins d'additions.

La programmation dynamique permet de donner en un temps quadratique une solution optimale au problème de l'alignement global de deux séquences.

Lorsqu'on reprend l'équation centrale de l'algorithme de programmation dynamique, on est tenté de modifier la fonction de coût total, pour exprimer le fait, que le score final se lit à la dernière itération de calcul, correspondant aux indices (n,m) des longueurs des deux séquences considérées. En effet le résultat sera lu sur dans cette dernière case; on peut alors prévoir d'aller par des verticales ou des horizontales jusqu'à la dernière case. La valeur du tableau modifié en (i,j) correspond alors au coût du chemin de à (n,m) passant par la case (i,j) et finissant uniquement par des insertions/délétions.

Qu'on se place dans le cas simple d'une pénalité des insertions/délétions linéaire, ou affine, cette modification de la programmation dynamique amène un gain de additions pour chaque couple d'indices (i,j). L'algorithme étant quadratique, ce gain est appréciable.

De plus, lorsqu'on se penche sur le flux des données pour le calcul d'une étape de l'algorithme, cette nouvelle formulation sépare certaines parties, ce qui pourrait permettre une parallélisation fine du calcul.

Utilisation de l'algèbre Max-Plus pour les algorithmes de programmation dynamique d'alignements de séquences.

Participants : Jean-Paul Comet et Jacques Henry

Mots-clés : Algèbre Max-Plus

Résumé : Les algorithmes de programmation dynamique de comparaisons de séquences s'expriment de façon linéaire dans l'algèbre Max-Plus.

L'algorithme de programmation dynamique pour l'alignement de séquences s'exprime de manière très simple grâce au formalisme de l'algèbre Max-Plus. Lorsque la fonction de coût d'une insertion est affine, lors de la comparaison de deux séquences de longueur n et m, l'algorithme revient à itérer n fois la multiplication (dans l'algèbre Max-Plus) d'un vecteur ligne $X_{l}$ de taille m par une matrice triangulaire supérieure $C_{l}$ :

$X_{l}$ = $X_{l-1}$ $\otimes$ $C_{l}$ .

Pour une matrice donnée, l'ensemble des matrices $C_{l}$ est très réduit: il y a autant de matrices que de lettres dans l'alphabet utilisé (4 dans le cas de l'ADN, et 20 dans le cas des protéines). Le semi-groupe des matrices engendré par l'ensemble des matrices $\{C_{a}, a \in {\cal A}\}$ est fini. Il s'en suit que le calcul des produits consécutifs des matrices peut être précalculé. Un automate de taille malheureusement très importante peut alors être construit pour répondre au problème de l'alignement de séquences par programmation dynamique.

Il a été implémenté dans le cas dégénéré de la recherche d'un motif avec erreurs. Les résultats pour des motifs de petites longueurs sont encourageants: le temps nécessaire pour la construction de l'automate, qui correspond au pré-calcul d'une partie du semi-groupe des matrices génératrices, n'est pas rédhibitoire, et est compensé par le gain de temps de l'algorithme de parcours de la banque de compléxite alors linéaire.

Ces travaux ont été réalisés en liaison avec S. Gaubert, du projet META2.

Applications de la programmation dynamique

Participant : Jacques Henry .

Mots-clés : programmation dynamique, problèmes aux limites elliptiques

Résumé : La méthode du plongement invariant appliquée spatialement permet d'obtenir une factorisation des problèmes aux limites elliptiques.

L'algorithme de comparaison de séquences de Needleman et Wunsch est fondé sur le principe de la programmation dynamique. Une réflexion a donc été menée sur les méthodes de programmation dynamique qui ont des applications dans des domaines très variés. Ceci conduit à reformuler le principe de la programmation dynamique comme un principe plus général d'élimination par substitution. Dans le cas de la minimisation d'un problème quadratique ce principe conduit à une méthode de directions conjuguées. Pour la résolution de systèmes linéaires il se traduit par la factorisation de Gauss par plongement dans une famille de systèmes linéaires de taille inférieure. L'algorithme QR résulte du plongement dans une famille de problèmes aux moindres carrés. L'application de la méthode à la factorisation de problèmes elliptiques dans des ouverts cylindriques s'est prolongée cette année en collaboration avec A. Ramos de l'universidade Complutense de Madrid. Une nouvelle collaboration débute sur le même sujet avec le CMAF de l'université de Lisbonne. La question étudiée est celle de la famille de problèmes conduisant après dicrétisation à une factorisation complète et non plus à une factorisation par blocs. On est alors conduit à définir une famillle de problèmes à n paramètres pour une équation dans un espace à n dimensions. Cela revient à faire un plongement invariant ``spatial'' à n-1 paramétres pour l'équation de Riccati trouvée par la méthode de plongement invariant à un paramètre. Pour la comparaison de séquences biologiques on a défini une analogie avec les problèmes précédents, qui est présentée au point .

Z-score

Participants : Jean-Christophe Aude , Jean-Paul Comet , Eric Glémet et Jean-Jacques Codani

Mots-clés : Statistique, loi des valeurs extrêmes, signification d'un alignement

Résumé : On montre que la distribution des Z-scores calculés sur des génomes complets suit une loi de Pareto.

Ce travail, mené conjointement avec le professeur Piotr Slonismki, a pour résultat la conception et l'implémentation de méthodes d'alignement de protéines (Z-score) conduisant à un modèle probabiliste fiable. Elles sont appliquées de manière intensive sur des génomes totalement séquencés, et permettent en particulier la construction rigoureuse et systématique de familles de protéines intra- et inter- organismes. Une étude, menée sur huit génomes microbiens complètement séquencés (S. Cervisiæ(levure), Synechocystis, Escherichia coli, Haemophilus Influenzæ, Methanococcus, Mycoplasma Genitalium et Pneumoniæ, B. Subtilis) montre que la distribution de ces familles et du nombre de séquences dans les familles ont une forte valeur prédictive, fonction de la taille du génome. Ces valeurs peuvent être considérées comme des lois de la génomique.

L'algorithme de Smith & Waterman ne fournit pas de probabilités sur les alignements et ne permet donc pas de générer un indice de dissimilarité nécessaire pour réaliser une classification des séquences. Nous avons donc adapté cet algorithme pour approcher la probabilité qu'un tel alignement soit dû ou non au hasard.

L'évaluation de la significativité d'un score est obtenu en calculant un Z-score. Pour deux séquences i et j ayant un alignement de score S(i,j), on compare ce score avec les scores obtenus à partir de i et d'un ensemble de séquences aléatoires. Cet ensemble respecte les critères suivants:

Les séquences sont de longueurs identiques à celle de j.
Les compositions en acides aminés sont identiques à celle de j.

Le calcul exact de la moyenne et de l'écart-type de cet échantillon n'étant pas possible, on utilise la méthode de Monte-Carlo pour estimer ces deux valeurs.

Cependant plusieurs problèmes subsistent: (i) existe-t-il un moyen de réduire le temps de calcul ? (ii) à partir de quelle valeur (cut-off) un Z-score est-il significatif ? (iii) existe-t-il une probabilité associée à ce Z-score?

Analyse de la variance du Z-score

L'étude de la variance du Z-score a permis d'établir une relation entre le nombre Nde séquences aléatoires à calculer, la variance estimée du Z-score et le Z-score:

$\begin{displaymath}\tilde{\sigma}^{2}(Z(i,j)) = (A \cdot N^{B}) \times \tilde{Z}(i,j) + \frac{1}{A' \cdot N + B'},\end{displaymath}$

où $\tilde{Z}(i,j)$ est le Z-score non-approximé (i.e. sans utiliser Monte-Carlo) entre les séquences i et j, A et les paramètres de la régression puissance, et ceux de la régression polynomiale inverse. L'utilisation de cette équation, en calculant itérativement le Z-score jusqu'à obtenir une variance satisfaisante, permet une réduction importante du temps de calcul et permet d'améliorer la qualité des résultats produits.

Détermination d'une valeur de cut-off

Ce que l'on appelle communément le ``cut-off'' est la valeur du Z-score au-delà de laquelle les alignements sont susceptibles de présenter un intérêt biologique. Jusqu'à présent, cette valeur était déterminée de manière expérimentale par les biologistes. Nous avons établi un modèle permettant de déterminer statistiquement ce ``cut-off''.

Pour ce faire, on a constitué un modèle de Bernouilli basé sur un échantillon de séquences réelles et un second échantillon de séquences aléatoires ayant une distribution en acides aminés réelle (afin de ne pas générer des séquences absurdes au sens biologique). On a pu ainsi déterminer une valeur du Z-score qui sépare l'échantillon aléatoire de l'échantillon réel; cette valeur (au voisinage de 8) est le ``cut-off'' recherché.

Détermination du modèle probabiliste sous-jacent

On cherche à associer au Z-score une loi de probabilité mesurant l'évènement: ``cet alignement est-il dû au hasard ?'' Le modèle pressenti dans le cas du Z-score est basé sur la loi de valeurs extrèmes type I (loi de Gumble). Nous avons montré que ce modèle n'était pas satisfaisant car il surestimait trop les Z-score supérieurs au ``cut-off''. Nous avons déterminé de manière empirique qu'une bonne approximation de la loi de distribution des Z-scores est la loi de Pareto.

Nous disposons désormais d'une probabilité. En utilisant la mesure de l'information de Wiener, nous pouvons en déduire un indice de similarité et a fortiori un indice de dissimilarité.

Classification

Participant : Jean-Christophe Aude

Résumé : Nous utilisons diverses méthodes de classification pour construire des familles de protéines à partir du Z-score.

La méthode que nous préconisons pour construire l'ensemble des familles recherchées est le ``composite clustering''. Ainsi, nous avons du combiner plusieurs méthodes pour trouver les familles de séquences. Les méthodes utilisées sont le clustering au sens de la théorie des graphes, la classification hiérarchique des clusters et la classification pyramidale pour chaque cluster.

Composite clustering

Le nombre important de séquences ne nous permet pas d'appliquer les méthodes standards de classification. En combinant plusieurs algorithmes de classification, nous sommes en mesure de résoudre ce problème. Les classes sont construites en utilisant la théorie des graphes et plus particulièrement la recherche de graphes fermés (ou composantes connexes). Les relations entre ces classes sont obtenues en appliquant les méthodes de classification hiérarchique sur les familles. La structure multi-domaines des familles composées de nombreuses séquences est explicitée en utilisant la classification pyramidale.

La première phase du composite clustering consiste à créer le graphe G(N,A), ou N l'ensemble des noeuds est composé des séquences et A l'ensemble des arcs des alignements. Les arcs sont valués suivant la valeur du Z-score de l'alignement. La méthode des composantes connexes à seuil variable permet de construire les familles. Cette méthode consiste à construire le graphe G'(N,A'), ou A' est le sous-ensemble de A tel que le Z-score de chaque arc est supérieur au seuil choisi. Les composantes connexes sont alors déterminées sur ce graphe.
L'analyse hiérarchique est basée sur la distance entre les composantes connexes. On procède alors par agrégation successive des individus pour créer un arbre de proximité de composantes. Le critère choisi est le min, c.à.d deux composantes sont agrégées ensemble en fonction de la proximité de ses deux plus proches individus.
La visualisation des familles de séquences est finalement réalisée en utilisant la classification pyramidale. La classification pyramidale est une généralisation du modèle hiérarchique. Une hiérarchie est un ensemble de partitions incluses les unes dans les autres. Une pyramide est un ensemble de recouvrements (deux ensembles ayant une intersection non nulle, et non inclus l'un dans l'autre). Cette caractéristique permet à une séquence d'appartenir à deux classes ou en d'autres termes d'être composée d'au moins deux domaines.

PMGD : Pyramidal Microbial Genome Databank

L'ensemble des ces méthodes est regroupé au sein de la librairie PLATOO (Post Lassap TOOlkit). Cette libraire a permis la réalisation d'un ensemble de programmes pour la construction de PMGD, une banque de famille de séquences homologues de l'ensemble des génomes microbiens séquencés. L'objet de cette banque est faciliter la consultation des résultats.

Optimisations de l'algorithme de Smith/Waterman

Participant : Andrzej Wozniak (action DYADE)

Résumé : On tire parti du parallélisme SIMD induit par le jeu d'instruction video du processeur UltraSparc (VIS) pour doubler la vitesse d'exécution de l'algorithme de Smith/Waterman.

Nous avons conçu et implanté une version optimisée de l'algorithme de Smith/Waterman (SW). Le principe est de vectoriser les opérations élémentaires nécessaires à la construction des cellules de la matrice calculée. En effet, les cellules composant une diagonale de la matrice peuvent être calculées au même instant. Notre implantation tire profit d'un parallélisme au niveau du jeu d'instruction des processeurs. L'algorithme de SW a donc été modifié pour en produire une version ``vectorisable''. Cette version a ensuite été instanciée sur le processeur UltraSparc (Sun), en utilisant le jeu d'instruction VIS (Video Instruction Set). VIS permet de décomposer un mot de 64 bits en plusieurs sous-mots de 8, 16 ou 32 bits, et d'appliquer un opérateur sur chaque sous-mot dans le même cycle d'horloge. Il s'agit d'un micro-parallélisme SIMD. Nous avons considéré des sous-mots de 16 bits, ce qui permet de calculer des scores qui couvrent la quasi-totalité des cas rencontrés. Plus important, l'implémentation assure l'intégrité des calculs en cas de saturation, c'est à dire qu'elle garantit qu'aucun dépassement de capacité ne perturbe le résultat. Cette propriété a été prouvée par J.P. Comet.

Un facteur deux a été gagné sur un UltraSparc à 167 MHZ. La performance atteinte est de 18 millions de cellules de matrices calculées par seconde (MMCS), ce qui est, à notre connaissance, l'implantation la plus rapide de SW sur station de travail. De plus, elle a la propriété d'être basée sur la comparaison de deux séquences. Elle n'est donc pas liée au cas particulier de ``database scanning'', mais garde au contraire son efficacité dans tous les cas de figures (banque contre banque, Z-score, ...). Elle a donc été intégrée dans LASSAP. Elle tire ainsi parti d'un parallélisme macroscopique. Des mesures de performances réalisées chez Sun France montrent que 200 MMCS sont atteints sur un serveur Enterprise6000 à 12 processeurs/167 MHZ. Une séquence de 300 acides aminés est comparée à SWISSPROT (release 33) en 30 secondes. Ces performances sont du même ordre de grandeur que celles du matériel spécialisé.

Parallélisme dans LASSAP

Participant : Eric Glémet

Mots-clés : Parallélisme, répartition de charge

Résumé : On montre montre que, pour assurer une répartion de charge efficace dans LASSAP, un découpage statique associé à une répartition dynamique est la solution qui permet la variation des paramètres tout en garantissant une efficacité de calcul.

Le module de flot de contrôle de LASSAP permet une implantation répartie des comparaisons ainsi que le pilotage de matériels spécialisés.

En d'autres termes, le module de flot de contrôle découpe et répartit automatiquement l'ensemble des comparaisons sur un réseau de processeurs, et ce quel que soit l'algorithme. La conception d'un tel parallélisme se heurte au triple problème: (i) les complexités des algorithmes sont très variables, (ii) la taille des problèmes est très variable, (iii) les moyens de calcul parallèle de la communauté sont très divers.

Une stratégie de découpage et de répartition qui minimise le temps de calcul, quels que soient les paramètres est théoriquement réalisable moyennant certaines approximations, en particulier en faisant abstraction de la localité des données et de la charge des processeurs.

En effet, pour estimer le temps de calcul parallèle, avec un nombre de processeurs donné, on est confronté aux problèmes d'hétérogénéité. L'hétérogénéité s'exprime au niveau:

des processeurs. On peut en effet disposer de plusieurs processeurs n'ayant pas la même rapidité.
des programmes. Des implantations particulières peuvent être faites pour certains types de processeurs. C'est l'algorithme ici qui est adapté à la machine (par exemple: cas de l'implantation VIS de l'algorithme SW sur UltraSparc).
des réseaux. C'est le cas lorsque l'on utilise plusieurs machines distantes.

Une étude des différents cas de figures entre découpage (statique ou dynamique) et répartition (statique ou dynamique) montre qu'un découpage statique associé à une répartition dynamique est la solution qui permet la variation des paramètres tout en garantissant une efficacité de calcul. Ceci nous a donc conduit à implanter des stratégies de découpage et de répartition de la charge en fonction de l'algorithme, du nombre de processeurs, et de la taille des données.

Cette stragégie a été implémentée dans LASSAP.

Recherche de similarités entre protéines

Participant : Pierre Nicodème (projet ALGO) .

Mots-clés : Blast, séries génératrices, multialignements

Résumé : Des résultats de marche aléatoire permettent de démontrer les formules de Karlin, bases de l'algorithme Blast.

P. Nicodème a soutenu sa thèse en Septembre 1997 sur la recherche de similarités entre protéines. Ce travail a été effectué en collaboration avec l'INRA-Toulouse qui développe la base ProDom (Protein Domains), où les familles de protéines sont regroupées suivant leurs domaines fonctionnels. Une similarité entre deux protéines s'exprime au moyen d'une fonction de score. On peut de la même manière exprimer une similarité entre une protéine et une famille de protéines. P. Nicodème utilise les formules de Karlin-Iglehart pour définir la pertinence probabiliste des scores obtenus par les algorithmes de recherche de similarité du type Blast. Il a regroupé dans sa thèse les différents développements mathématiques permettant d'obtenir ces formules. P. Nicodème utilise les séries génératrices dans le cas d'un modèle de scores discret, qui est celui utilisé par les biologistes, ce qui simplifie beaucoup la démonstration des formules de Karlin à partir de l'identité de Spitzer, qui est un résultat classique de marche aléatoire. Il avait par ailleurs développé précedemment le logiciel BlastMultAl qui permet des recherches de similarités entre une protéine et une famille de protéines. Il a calibré ce logiciel pour les domaines de ProDom33 en utilisant la méthode de clumping-declumping et d'approximations Poissonniennes de Waterman et Vingron.

Précédent : Logiciels Remonter : Action ATGC, Action Transversale Génome Suivant : Actions industrielles