Site WWW de Laurent Bloch

Plan du cours d'algorithmes pour la biologie (BNF 103)

2010-03-07T17:19:00Z

Le cours suivra le plan suivant :

Définitions, notations, pseudocode

Exemple, Pseudo-code, pseudo-code->scheme, recherche dans une liste, complexité, notations asymptotiques.

Automates et machines de Turing

À l'attaque des vrais fondements.

Automates II : non-déterminisme

Après les automates à états finis déterministes, les automates non-déterministes, qui sont d'ailleurs équivalents.

Expressions régulières

La suite logique des automates, auxquels elles sont équivalentes.

Accélération d'un algorithme

Recherche du maximum d'un ensemble de valeurs, memoization de Fibonacci, recherche de segments maximaux O(n³), O(n²), O(n). On attaque l'algorithmique proprement dite.

Accélération de l'algorithme de Fibonacci

Un exemple.

Tris

Algorithmes de tri par insertion, recherche de maximum dans une liste, tri fusion, tri rapide, tri linéaire.

Structures de données élémentaires

Les pile, les files, le tas et tris par tas, les arbres, arbres binaires de recherche, tables de hashage.

Recherche de mots dans un texte

Force brute, Knuth-Morris-Pratt (KMP), Karp Rabin.

Le précédent en style récursif

Instrumentation d'un programme pour l'observer : KMP

Une variation autour du précédent.

Recherche de mots dans un texte (suite)

Structures de données Les arbres de mot-clés (keyword trees), les arbres de suffixes (suffix trees).

Recherche d'un motif dans un texte

Prosite, les expressions régulières, grep, les automates à état fini, expression régulières -> automates à état finis indéterministes, déterminisation d'un automate à état fini, minimisation d'automates à état finis.

Algorithmes sur les graphes

Le métabolisme, représentation de graphe, recherche de plus court chemin, recherche d'arbre couvrant minimal.

Algorithme de Needleman et Wunsch

Alignement de séquences.

Programmation dynamique

Recherche de séquences hydrophobes, comparaison de séquence par l'algorithme de Needleman et Wunsch, par l'algorithme de Smith et Waterman, optimisation de Gotoh, multiplication de matrices : minimisation des calculs.

Recherche dans des banques de données

Blast, Fasta, les matrices PAM.

Chaines de Markov

Prédiction de modèles de gènes, prédiction de structure, représentation de motifs fonctionnels.

Des problèmes difficiles : NP-complétude

Problème de la clique, Sat, voyageur de commerce.

Solution presque optimales

Simulated anealing, algorithmes génétiques...

Algorithmes de tri

Laurent Bloch, William Saurin — 2009-06-22T14:14:00Z

Auteurs : William Saurin, Laurent Bloch

#Sommaire-

Recherche dichotomique

Aurait-on-pu aller plus vite dans la recherche d'un élément dans une liste ?

Comment cherche-t-on un mot dans un dictionnaire ? on regarde les mots de la page qui est au milieu du dictionnaire et on se demande si celui que l'on cherche est dans la page, s'il plus grand que le dernier ou s'il est plus petit que le premier. S'il est dans la page : c'est trouvé, s'il est plus petit on recommence sur la moitié inférieure du dictionnaire et s'il est plus grand on recommence sur la moitié supérieure.

Algorithme dichotomique

Nom de l'algorithme : dichot Donnée d'entrée : vecteur trié t et valeur v Résultat : numéro dans la liste petit <- 0 grand <- longueur de t -1 pivot <- Entier(grand/ 2) si v = t[grand] alors retourner grand sinon si v =t[petit] alors retourner petit sinon tant que pivot != petit et t[pivot] != v faire si t[pivot] < v alors petit <- pivot sinon grand <- pivot finsi pivot <- Entier((grand - petit) / 2) fin tant que si pivot != petit alors retourner pivot sinon retourner -1 fin si

Exemple

Rechercher 176 dans le tableau suivant :


0	1	2	3	4	5	67	8	9	10	11
1	17	43	89	100	113	156	176	231	248	301

petit					pivot					grand
0	1	2	3	4	5	67	8	9	10	11
1	17	43	89	100	113	156	176	231	248	301

					petit		pivot			grand
0	1	2	3	4	5	67	8	9	10	11
1	17	43	89	100	113	156	176	231	248	301

Combien de fois fait on le test “grand != petit et t[pivot] != v”

En fait à chaque étape on coupe l'espace de recherche (le tableau en 2) . On cesse de faire des tests lorsque l'espace de recherche contient au plus un ou deux éléments.

étape	1	2	3	4	...	n
taille du tableau	t	t/2	t/4	t/8	...	t/2ⁿ⁻¹

A l'étape n (où on s'arrête) on a une taille 2<= t/2ⁿ⁻¹ et donc

2ⁿ<= t et donc n <= log₂t

Selon l'algorithme le nombre de test est :

En Scheme :

(define (dichot t v) ; t est un vecteur triée v une valeur (let* ((petit 0) (grand (- (vector-length t) 1)) (pivot (quotient (+ grand petit) 2))) (if (= v (vector-ref t grand)) grand (if (= v (vector-ref t petit)) petit (let loop () (if (and (not (= petit pivot)) (not (= v (vector-ref t pivot)))) (begin (if (< (vector-ref t pivot) v) (set ! petit pivot) (set ! grand pivot)) (set ! pivot (quotient (+ grand petit) 2)) (loop)) (if (not (= petit pivot)) pivot -1)))))))

Tri par insertion

Un tableau A :

0	1	2	3	4	5	6	7
100	10	8	56	78	4	67	34

Nous nous proposons de trier un tableau V représenté par un vecteur de longueur Vlen.

L'idée est simple : à un certain moment de l'algorithme la région qui va de 0 à j − 1 dans le tableau est déjà triée. Il suffit pour avoir une séquence triée de 0 à j d'insérer V[j] dans cette séquence. Insérer c'est trouver le premier élément dans la séquence de 0 à j−1 tel que cet élément soit plus grand que V[j], on décale alors tous les éléments y compris celui là d'un rang vers le haut et on insère l'élément V[j] dans le trou.

Il y a un piège : si on décale tous les éléments vers la droite, à la fin du décalage V[j] doit contenir V[j−1]. Il faut donc avoir sauvegardé V[j] avant le décalage : il nous faut une variable auxiliaire.

On peut observer que la séquence qui va de 0 à 0 dans un tableau est toujours triée.

Pour une étape donnée il faut faire :

trouver où insérer ;
sauvegarder V[j] ;
décaler ;
insérer V[j] dans le trou.

En somme :

0	1	2	3	4	5	6	7
8	10	56	100	78	4	67	34
				__

0	1	2	3	4	5	6	7
8	10	56	100	78	4	67	34
			^	__

puis

0	1	2	3	4	5	6	7
8	10	56	100		4	67	34
			^	__

0	1	2	3	4	5	6	7
8	10	56		100	4	67	34
			^	__

0	1	2	3	4	5	6	7
8	10	56	78	100	4	67	34

Pseudo-code

Nom de l'algorithme : tri-insertion Donnée : un tableau V Résultat : V est trié pour j allant de 1 à (longueur de V) -1 faire i <- 0 tant que i < j et V[i] <= V[j] faire i <- i + 1 fait rotation-droite (V, i , j - 1) fait retourner V Nom de l'algorithme : rotation-droite Données : un tableau V, i, j Soit Vlen <- longueur(V) tmp <- V[j+1] k <- j tant que k > i faire V[k+1] <- V[k] k <- k - 1 fait V[i] <- tmp

Traduisons en Scheme :

(define (rotation-droite V i j) (let ((Vlen (vector-length V))) (if (>= j (- Vlen 1)) 'impossible (let ((tmp (vector-ref V (+ j 1)))) (do ((k j (- k 1))) ((< k i) (vector-set ! V i tmp)) (vector-set ! V (+ k 1) (vector-ref V k)))))))

(define (tri-insert V) (let ((Vlen (vector-length V))) (do ((j 1 (+ j 1))) ((> j (- Vlen 1)) V) (let loop ((i 0)) (if (and (<= (vector-ref V i) (vector-ref V j)) (< i j)) (loop (+ i 1)) (rotation-droite V i (- j 1)))))))

Complexité

Le nombre d'opérations de comparaisons d'éléments de tableau de cet algorithme est en O(n²) en effet il y a n−1 étape à chacune d'entre elles on fait au plus le nombre de comparaison correspondant au nombre d'éléments de la séquence trié de l'étape, pour l'étape 1 1, pour l'étape 2 2, et ainsi de suite soit 1+2+3+...+n−1 si n est la taille du tableau et donc au plus (n−1)n/2 comparaisons à chaque étape il y a au plus i décalage si i est le numéro de l'étape et donc de la même façon (n−1)n/2 décalage élémentaires.

Tri fusion

L'idée est la suivante. Supposons que nous ayons deux ensembles de valeurs triées (dans une liste par exemple) pour fabriquer une liste triée à partir des deux listes triées il suffit de se dire que si la première liste est vide le résultat est la seconde, si la seconde est vide, le résultat est la première si aucune des deux n'est vide il suffit de prendre la tête de la liste qui est la plus petite et de la concaténer à la liste résultant de la fusion du reste de cette liste et de l'autre liste. En Scheme :

(define (fusion l1 l2) (if (null ? l1) l2 (if (null ? l2) l1 (if (< (car l1) (car l2)) (cons (car l1) (fusion (cdr l1) l2)) (cons (car l2) (fusion l1 (cdr l2)))))))

Ainsi on voit que pour trier une liste il suffit de la séparer en deux sous listes, de trier chacune de ces listes et de les fusionner. Sauf si la liste est vide au départ, dans ce cas elle est déjà triée !

Séparer une liste en deux sous-listes et appliquer une fonction f à chacun de ces sous-listes et enfin appliquer une fonction g à ces deux sous-listes peut s'écrire :

(define (separe-et-applique l f g l1 l2) (if (null ? l) (g (f l1) (f l2)) (separe-et-applique (cdr l) f g l2 (cons (car l) l1))))

Pour trier par fusion la liste l il suffirait d'appeler separe-et-applique à a avec comme fonction g fusion et comme fonction f une fonction qui trie une liste !

(separe-et-applique l trier fusion '() '())

et en somme :

(define (fusion l1 l2) (if (null ? l1) l2 (if (null ? l2) l1 (if (< (car l1) (car l2)) (cons (car l1) (fusion (cdr l1) l2)) (cons (car l2) (fusion l1 (cdr l2)))))))(define (separe-et-applique l f g l1 l2) (if (null ? l) (g (f l1) (f l2)) (separe-et-applique (cdr l) f g l2 (cons (car l) l1))))(define (tri-fusion l) (if (or (null ? l) (null ? (cdr l))) l (separe-et-applique l tri-fusion fusion '() '())))

cons ?—>

Combien de `cons` ?

A la première étape on fera pas plus de cons qu'il y a d'éléments dans la liste de départ -1. Mais on aura aussi du faire le nombre de cons nécessaire à trier les deux sous-liste qui ont chacune pas plus que n/2 +1 éléments.

N(n) <= n−1 + 2 * (N(n/2 +1))

N(n) <= n−1 + 2 * n/2 * 4(N(n/2 +1 / 2 +1))

N(n) <= n−1 + n + 4

Faire un arbre !!!

(Quicksort)—>

Tri rapide

L'algorithme Quicksort a été publié par C.A.R. Hoare¹ en 1962.

Une autre façon de faire consiste à choisir un élément du tableau à trier, nous appellerons cet élément le pivot, puis à s'assurer que tous les éléments plus petits ou égaux au pivot sont au début du tableau et les éléments plus grand ou égaux à la fin. Il y a de la sorte deux région créées dont tous les éléments de la première sont plus petits ou égaux à ceux de la seconde. De la sorte le tableau est partiellement trié, il suffit alors de trier la partie basse et la partie haute.

On peut écrire qu'à chaque étape on fera :

choisir le pivot ;
créer deux régions une au début qui ne contienne aucun élément plus grand que pivot, une à la fin qui ne contienne aucun élément plus petit que pivot ; soit q la position fin de la première région :
- trier les éléments du début a q,
- trier les éléments de q +1 à dernier élément ;

Supposons qu'il s'agisse de trier les éléments du tableau entre les positions imin et imax (comprises), on pourrait dire que le pivot est imin puis on examinerait en descendant tous les éléments du tableau de imax vers imin jusqu'à en avoir trouvé un plus petit ou égal aupivot appelons jmax cette position. On peut ensuite parcourir le tableau de imin jusqu'à ce qu'on ait trouvé un élément plus grand ou égal que le pivot ; appelons jmin cette position. Si jmin est inférieur à jmax alors on échange les valeurs en jmin et jmax on incrémente jmin et on décrémente jmax puis on recommence à balayer vers le bas (pour jmax) et vers le haut (pour jmin). Si jmin n'est pas inférieur à jmax alors il est clair que jmax est le haut de la région de valeur inférieures au pivot.

On peut alors trier entre imin et jmax et entre jmax +1 et imax.

Nom de l'algorithme : tri-rapide Donnée d'entrée : un vecteur v, imin et imax les bornes dans lesquelles trier si imin < imax alors q <- partition(v, imin, imax) tri-rapide(v, imin, q) tri-rapide(v, q+1, imax) finsi

avec l'algorithme suivant pour partition :

Nom de l'algorithme : partition Données d'entrée : un vecteur v, imin et imax les bornes dans lesquelles partitionner x <- v[imin] tant que vrai faire tant que v[imax] > x faire imax <- imax -1 fait tant que v[imin] < x faire imin <- imin + 1 fait si imin < imax alors echanger v[imax] et v[ iman] imax <- imax - 1 imin <- imin + 1 sinon retourner imax finsi fait

En Scheme :

(define (tri-rapide V)
 (define (tri-aux V imin imax)
 (if (< imin imax)
 (let ((q (partition V imin imax)))
 (tri-aux V imin q)
 (tri-aux V (+ q 1) imax)))
 V)
 (tri-aux V 0 (- (vector-length V) 1)))
(define (partition V imin imax)
 (let ((x (vector-ref V imin)))
 (let boucle ()
 (do ((i imax (- i 1)))
 ((<= (vector-ref V imax) x)) (set! imax i))
 (do ((i imin (+ i 1)))
 ((>= (vector-ref V imin) x)) (set! imin i))
 (if (< imin imax)
 (begin
 (permute V imax imin)
 (set! imax (- imax 1))
 (set! imin (+ imin 1))
 (boucle))
 imax))))
(define (permute V i j)
 (let ((tmp (vector-ref V i)))
 (vector-set! V i (vector-ref V j))
 (vector-set! V j tmp)))

Télécharger

En voici une version compilée et générique, qui reçoit en paramètre le prédicat qui définit la relation d'ordre à appliquer pour le tri, et qui en infère aussi le type des éléments à trier, nombres ou chaînes de caractères :

(module tri-rapide-2
 (main def-tri))
;; ./bin/make-tri-rapide "string<=?" 5 8 2 0 78 0 1 12 42 2 4
;; #(0 0 1 12 2 2 4 42 5 78 8)
;; ./bin/make-tri-rapide ">=" 5 8 2 0 78 0 1 12 42 2 4
;; #(78 42 12 8 5 4 2 2 1 0 0)
(define (def-tri Args)
 (let* ((nom-predicat (cadr Args))
 (le-tri
 (make-tri-rapide
 (eval (string->symbol nom-predicat))))
 (le-vecteur
 (list->vector
 (cond
 ((string=? (substring nom-predicat 0 1) "s")
 (cddr Args)) ;; on trie des chaînes de caractères
 (else ;; on trie des nombres
 (map string->number (cddr Args)))))))
 (print (le-tri le-vecteur))))
;; (define mon-tri (make-tri-rapide (eval (string->symbol "string<=?"))))
;; (mon-tri '#("Alice" "Julot" "Bernard" "Élodie"))
;; #(Alice Bernard Julot Élodie)
(define (make-tri-rapide predicat)
 (define (tri-aux V imin imax)
 (if (< imin imax)
 (let ((q (partition V imin imax)))
 (tri-aux V imin q)
 (tri-aux V (+ q 1) imax)))
 V)
 (define (tri-rapide V)
 (tri-aux V 0 (- (vector-length V) 1) ))
 (define (permute V i j)
 (let ((tmp (vector-ref V i)))
 (vector-set! V i (vector-ref V j))
 (vector-set! V j tmp)))
 (define (partition V imin imax)
 (let ((x (vector-ref V imin)))
 (let boucle ()
 (do ((i imax (- i 1)))
 ((predicat (vector-ref V i) x)
 (set! imax i)))
 (do ((i imin (+ i 1)))
 ((predicat x (vector-ref V i) )
 (set! imin i)))
 (if (< imin imax)
 (begin
 (permute V imax imin)
 (set! imax (- imax 1))
 (set! imin (+ imin 1))
 (boucle))
 imax))))
 tri-rapide) ;; make-tri-rapide renvoie une procédure, 
 ;; qui est le programme de tri voulu

Télécharger

Tri linéaire

Il s'agit de trier un tableau V mais nous savons dans quel intervalle sont présentes les valeurs de V.

On peut alors créer un tableau C dont les indices vont de valeur la plus petite de V (vmin) à valeur la plus grande de V (vmax). C est initialisé à 0.

On parcourt V et pour chaque valeur de V[i] on fait C[V[i]] ← C[V[i]] +1. A la fin dans c[v[i]] on a le nombre de fois que la valeur C[V[i]] est présente dans V. On peut alors calculer les sommes cumulées de C[vmin] ... C[i] on a maintenant en C[i] le nombre de valeurs de V qui sont inférieures ou égales à i. On parcourt une fois V et pour chaque V[i] rencontré on écrit sa valeur en C[V[i]] dans un tableau B et on fait C[V[i]] ← C[V[i]]−1. À la fin B est une version triée de A.

Nom de l'algorithme : tri-linéaire Données d'entrée : v un vecteur Donnée en sortie : b un vecteur qui contient toutes les valeurs de v triées Créer un vecteur c de longueur valeur maximum de v - valeur minimum de v + 1 initialisé à 0 pour toutes ses valeurs. pour i allant de 0 à longueur de v -1 faire c[v[i] - vmin] <- c[v[i] - vmin] + 1 fait pour i allant de 1 à longueur de c -1 faire c[i] <- c[i] + c[i-1] fait créer un vecteur b de longueur longueur de c pour i allant de longueur de v -1 à 0 faire b[c[v[i] -vmin] - 1] <- v[i] c[v[i]-vmin] <- c[v[i]-vmin] - 1 fait retourner b

Combien d'opérations ?

En Scheme :

(define (rmax v) (let loop ((rmax 0) (i 1)) (if (>= i (vector-length v)) rmax (loop (if (> (vector-ref v i) (vector-ref v rmax)) i rmax) (+ i 1)))))(define (rmin v) (let loop ((rmin 0) (i 1)) (if (>= i (vector-length v)) rmin (loop (if (< (vector-ref v i) (vector-ref v rmin)) i rmin) (+ i 1)))))

(define (tri-lineaire v) (let ((vmin (vector-ref v (rmin v))) (vmax (vector-ref v (rmax v)))) (let ((b (make-vector (vector-length v))) (c (make-vector (+ 1 (- vmax vmin)) 0))) (do ((i 0 (+ i 1))) ((>= i (vector-length v)) 'done) (vector-set ! c (- (vector-ref v i) vmin) (+ 1 (vector-ref c (- (vector-ref v i) vmin))))) (do ((i 1 (+ i 1))) ((>= i (vector-length c)) 'done) (vector-set ! c i (+ (vector-ref c (- i 1)) (vector-ref c i)))) (do ((i 0 (+ i 1))) ((>= i (vector-length v)) b) (vector-set ! b (- (vector-ref c (- (vector-ref v i) vmin)) 1) (vector-ref v i)) (vector-set ! c (- (vector-ref v i) vmin) (- (vector-ref c (- (vector-ref v i) vmin)) 1))))))

Tas et tri par tas

Un tas est un objet qui peut être vu comme un arbre binaire complet : à part peut-être un seul de ses sous arbres tous ont soit un ou deux fils qui sont des arbres-vides soit deux fils qui ne sont pas des arbres vides. Tous les sous-arbres non vides se présentent sur deux niveaux au plus. Enfin lorsqu'on examine deux sous arbres vides dont l'un est à gauche de l'aure, le sous arbre de gauche est toujours à une profondeur supérieure ou égale à celui de droite. La valeur d'un sous arbre est toujours plus élevée que les valeurs contenues dans ses sous arbres droit et gauche.

Représentation d'un tas : on peut le représenter par un tableau (de longueur l) et par une valeur qui représente le nombre de sous-arbres non vides du tas. Les valeurs du tas seront stockées dans les cellules 1 à taille du tas (attention à la taille maximum du tas qui doit etre plus petite d'un que la taille du tableau). De la sorte la valeur du fils gauche de l'arbre dont la valeur sera en position i du tableau sera stockée en 2i et le fils droit en 2i+1, de la même manière le père d'un noeud en position i sera en Entier(i/2).

Nom de l'algo : père Donnée : i Résultat : la position du père de i retourner Entier(i/2)

De le même façon :

Nom de l'algo : fils-gauche Entrée : i Résultat : la position du fils gauche de i retourner 2 * i

et fils droit :

Nom de l'algo : fils-droit Donnée : i Résultat : la position du fils droit de i retourner (2 * i) + 1

(define (tas:creer taillemax) (vector 'tas 0 (make-vector (+ 1 taillemax))))(define (tas:taille T) (vector-ref T 1))(define (tas:tableau T) (vector-ref T 2))(define (tas:pere i) (quotient i 2))(define (tas:fils-gauche i) (* i 2))(define (tas:fils-droit i) (+ 1 (* 2 i)))(define (tas:permute ! T i j) (let ((tmp (vector-ref T i))) (vector-set ! T i (vector-ref T j)) (vector-set ! T j tmp)))

Problème : comment ajouter un élément à un tas. On augmente la taille du tas de 1, on positionne le nouvel élément à l'indice taille dans le tableau et on le laisse remonter (en échangeant sa valeur avec son père) tant que son père existe (on n'est pas à la racine) ou que celui-ci à une valeur inférieure à la valeur du noeud :

Nom de l'algo : augmenter-tas Données : T, un tas, et v, une valeur Résultat : le tas comportant la nouvelle valeur taille de T <- taille de T + 1 i <- taille de T T [i] <- v tant que i > 1 et T [i] > T[pere(i)] faire échanger T [i] et T[pere(i)] i <- pere (i) fait

La complexité de l'algorithme ? log₂(taille(t)).

On peut envisager d'insérer un nouvel élément dans le tas à la place du sommet de la façon suivante :

Nom de l'algo : insérer-tas Données : v, T une valeur et un tas. Résultat : T le tas T[1]<-v i<-1 tant que (fils-gauche(i) <= taille de T et T[fils-gauche(i)] > T[i]) ou (fils-droit[i] <= taille de T et T[fils-droit(i)]> T[i]) faire si fils-droit(i) > taille de T ou T[fils-droit(i)] <= T[fils-gauche(i)] alors echanger T[fils-gauche(i)] et T[i] sinon echanger T[fils-droit(i)] et T[i] finsi retourner T

Exemple :

Figure 1 : Tas = 10 sommets

(Cette figure et la suivante sont empruntées à H. W. Lang de l'université de Flensburg en Allemagne, qui propose un excellent site sur le sujet : http://www.inf.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm)

Peut-on extraire l'élément maximum du tas et conserver au tas sa structure ? Il suffit de prendre l'élément en position 1 (que l'on extrait), de placer en position 1 l'élément qui est en position taille du tas puis de le faire descendre tant qu'il n'est pas en bas et que l'un de ses fils plus grand que lui après avoir réduit d'un la taille du tas. Il descend de la sorte en étant échangé avec son fil de valeur la plus élevée.

Nom de l'algo : extraire-tas Données : T un tas Résultat : v le plus grand élément de T v<- tableau de T[1] taille de T <- taille de T -1 inserer-tas(T[taille de T + 1], T) retourner v

On peut alors utiliser ces deux algorithmes pour trier un tableau.

(a) (b) (c)

(d) (e)

Figure 2 : Extraire l'élément le plus grand et restaurer la structure de tas

Nom de l'algo : Tri-par-tas Données : V un tableau Résultat : V trié. t <- creer un tas contenant un tableau de taille longueur de V + 1 taille de t = 0 pour i allant de 0 à longueur(V) - 1 faire augmenter-tas(t, V[i]) fait pour i allant de 0 à longueur de V - 1 faire V[longueur(V) - 1 - i] <- extraire-tas(t) fait retourner V

Figure 3 : Représentation par un tableau du tas de la figure 1

(define (tas:ajoute ! T v) (vector-set ! T 1 (+ (tas:taille T) 1)) (let ((i (tas:taille T))) (vector-set ! (tas:tableau T) i v) (let boucle ((i i)) (if (and (> i 1) (> (vector-ref (tas:tableau T) i) (vector-ref (tas:tableau T) (tas:pere i)))) (begin (tas:permute ! (tas:tableau T) i (tas:pere i)) (boucle (tas:pere i))) T))))(define (tas:insere-sommet ! T v) (vector-set ! (tas:tableau T) 1 v) (let boucle ((i 1) (g 2) (d 3)) (if (or (and (<= g (tas:taille T)) (< (vector-ref (tas:tableau T) i) (vector-ref (tas:tableau T) g))) (and (<= d (tas:taille T)) (< (vector-ref (tas:tableau T) i) (vector-ref (tas:tableau T) d)))) (let ((m (if (> d (tas:taille T)) g (if (< (vector-ref (tas:tableau T) g) (vector-ref (tas:tableau T) d)) d g)))) (tas:permute ! (tas:tableau T) i m) (boucle m (tas:fils-gauche m) (tas:fils-droit m))))))(define (tas:extraire T) (let ((res (vector-ref (tas:tableau T) 1)) (nouveau-sommet (vector-ref (tas:tableau T) (tas:taille T)))) (vector-set ! T 1 (- (tas:taille T) 1)) (tas:insere-sommet ! T nouveau-sommet) res))

On peut utiliser ces algorithmes pour déterminer quels sont les k plus petits éléments d'un flux de données de données. On commencerait construire un tas avec les k premiers éléments du tas. A partir de l'élément k+1 du flux on pourrait pour chacun des éléments suivants tester s'il est plus petit que l'élément au sommet du tas et si c'est le cas l'insérer au sommet du tas.

(define (trier-liste L) (let ((T (tas:creer (length L)))) (let loop ((L L)) (if (not (null ? L)) (begin (tas:ajoute ! T (car L)) (loop (cdr L))))) (let loop ((L '())) (if (zero ? (tas:taille T)) L (loop (cons (tas:extraire T) L))))))(define (trier-vecteur V) (let ((T (tas:creer (vector-length V)))) (do ((i 0 (+ i 1))) ((>= i (vector-length V))) (tas:ajoute ! T (vector-ref V i))) (do ((i (- (vector-length V) 1) (- i 1))) ((< i 0) V) (vector-set ! V i (tas:extraire T)))))(define (n-min V n) (if (<= (vector-length V) n) (trier-vecteur V) (let ((T (tas:creer n))) (do ((i 0 (+ i 1))) ((>= i n) T) (tas:inserer-sommet ! T (vector-ref V i))) (do ((i n (+ i 1))) ((>= i (vector-length V)) T) (entasser T (vector-ref V i))) (let ((W (make-vector n))) (do ((i 0 (+ i 1))) ((>= i n) W) (vector-set ! W (- (- n 1) i) (tas:extraire T)))))))

Exercice : Ecrire l'algorithme décrit informellement ci-dessus en supposant qu'une fonction fct rend des valeurs successives ou bien qu'elle rend -1 comme dernière valeur.

1
Le Britannique C. Antony R. Hoare est à l'origine de quelques contributions de premier plan. Des études universitaires de tonalité plutôt littéraire lui ont donné l'occasion de partir en stage en 1960 dans le laboratoire d'Andreï Nikolaiévitch Kolmogorov à l'Université de Moscou pour travailler sur un projet de traduction automatique. Pour réaliser un dictionnaire électronique nécessaire à ce projet (par ailleurs vite abandonné) il inventa l'algorithme de tri « Quicksort » que tout étudiant en informatique se doit d'avoir étudié en détail et programmé. Ce stage lui avait donné le goût de la programmation, mais, comme il le raconte lui-même avec humour, il eut la chance que sa formation initiale lui fermât les portes du UK National Physical Laboratory, et il entra dans l'industrie. C'est là qu'il eut l'occasion de participer au développement de systèmes d'exploitation, domaine pour lequel il élabora la méthode des moniteurs de Hoare afin de contrôler les accès concurrents et exclusifs à des ressources partagées et la théorie des processus séquentiels communicants, encore aujourd'hui le seul modèle complet et cohérent qui excède réellement le modèle de von Neumann.

This document was translated from L^AT_EX by H^EV^EA.

Recherche de mots dans un texte, KMP

Laurent Bloch, William Saurin — 2009-05-08T19:33:00Z

L'algorithme de Knuth-Morris-Pratt

William Saurin, Laurent Bloch

Le 8 mai 2009

#Sommaire-

Nous nous intéressons à déterminer si un mot (ici un mot est une simple chaîne de caractères) figure dans un texte. Nous noterons Mot le mot et nous considérerons que les éléments successifs du mot recherché sont stockés en Mot[0], Mot[1], ... Mot[m-1] Le texte sera noté Texte et ses éléments successifs seront stockés en Texte[0], Texte[1], ... Texte[t-1].

L'indice j correspond à une progression dans le mot et l'indice i à une progression dans le texte. À un moment donné on compare le caractère d'indice i+j du texte avec le caractère d'indice j du mot.

Algorithme force brute.

On peut imaginer positionner successivement la sous chaine sous le texte de manière à ce que sa première position (celle indexée par 0) soit successivement sous les chacune des positions du texte. On vérifie alors la coincidence des différents caractères superposés.

Algorithme : brute-force Données : Texte, Mot ; le texte, et le mot recherché Résultat : la première position où le mot a été trouvé, ou bien -1. Soient Ltexte la longueur du Texte Lmot la longueur du Mot pour i allant de 0 à (Ltexte - Lmot) faire j <- 0 tant que j < Lmot et Texte[i+j] = Mot[j] faire j <- j+1 fait si j = Lmot alors retourner i finsi fait si j = Lmot alors retourner i sinon retourner -1

Complexité de l'algorithme : l × m

En Scheme :

(module force-brute
 (main get-args))
(define (get-args args)
 (print (force-brute (cadr args)
 (caddr args))))
;; les instruments
(define (force-brute Texte Motif)
 (let* ((L-Texte (string-length Texte))
 (L-Motif (string-length Motif))
 (DerPos-Texte (- L-Texte 1)))
 (let boucle1 ((i 0)
 (j 0))
 (let boucle2 ()
 (if (and (< j L-Motif)
 (char=? (string-ref Texte (+ i j))
 (string-ref Motif j)))
 (begin
 (set! j (+ j 1))
 (boucle2))))
 (cond ((= j L-Motif) i)
 ((< i DerPos-Texte)
 (boucle1 (+ i 1) 0))
 (else -1)))))

Télécharger

Exposé de l'algorithme KMP

L'idée de l'algorithme est simple : lorsque nous ne réussissons plus à faire progresser l'indice j vers la droite, en fait nous faisons progresser l'indice i mais dans ce cas nous réexaminons à nouveau des caractères du texte. Puisque nous les avons déjà vus ne pourrions nous pas nous en passer ?

La deuxième idée est que puisque les k+1 caractères du texte, de rang i à i+k, que nous venons d'examiner sont « bons » alors ce sont aussi les k+1 premiers caractères du mot : nous pouvons savoir ce qui nous arrivera d'avance uniquement en examinant d'abord les caractères du mot. Par exemple, si les n derniers caractères parmi ces k+1, soit les caractères du mot de rangs k−n+1 à k, se trouvent constituer un préfixe du mot, nous pouvons recommencer la comparaison dans le texte à la position i+k−n+1, et comme nous savons que les caractères i+k−n+1 à i+k sont bien égaux aux n premiers caractères du mot, en fait nous allons recommencer la comparaison avec le caractère de rang i+k+1 dans le texte, et k−n+1 dans le mot.

La troisième idée est celle-ci : la seconde idée consiste à exploiter le fait que si nous examinons les différents préfixes du mot, de longuurs croissantes, certains de ces préfixes vont avoir des suffixes identiques à un certain préfixe du mot lui-même, ce qui est la circonstance qui autorise le recours à la seconde idée exposée ci-dessus.

Reprenons : si les k+1 caractères de rangs i à i+k étaient « bons », c'est-à-dire égaux aux k+1 premiers caractères du mot, on peut en tirer des conséquences. Par exemple, si aucun des caractères de i+1 à i+k du texte n'est égal au premier caractère du mot, on en déduit qu'il est inutile de commencer une comparaison à ces emplacements, on peut sauter tout de suite à la position i+k+1 du texte. Par contre, si un suffixe de la chaîne Texte[i..i+k], disons Texte[i+k-m..i+k] est un préfixe du mot, alors il faut reprendre la comparaison à la position Texte[i+k+1] du texte et à partir de la position Mot[k-m], puisque l'on sait déjà que la comparaison sera positive pour les m+1 premiers caractères.

Si la comparaison était positive pour les k+1 premiers caractères du mot, soit Texte[i..i+k] = Mot[0..k], et si en outre Texte[i+k-m..i+k] = Mot[k-m..k], alors c'est que Mot[0..m]= Mot[k-m..k]. Or ce résultat concerne uniquement Mot, est peut être pré-calculé.

Dans le cas de notre exemple : i=2, k=4, m=2.

Imaginons donc que pour tous les mots allant de 0 à k dans la chaîne nous connaissions la positions la plus grande j telle que le mot allant de 0 à j soit identique à celle allant de k−j à k. Considérons que prec soit une fonction telle que prec(k) soit la valeur de j en question. Imaginons que nous ayons examiné tous les caractères du texte allant de i à i+k avec succès, c'est à dire tous ceux de la chaine de 0 à k. Nous cherchons donc à trouver quel caractère de la chaîne peut-être mis en correspondance avec le caractère i+k+1 du texte. Le caractère k+1 de la chaîne est un bon candidat mais s'il n'est pas identique à celui du texte alors le caractère prec(k)+1 de la chaîne en est autre bon candidat, si ce caractère n'est pas identique alors prec(prec (k))+1 est à nouveau un bon candidat, et ainsi de suite jusqu'à ce qu'on ait trouvé une identité ou jusqu'à ce que prec(...(prec(k))...)+1 vaille 0 et ait échoué. Dans le premier cas on examinera le caractère suivant du texte. Dans tous les cas on examine le caractère suivant du texte.

Figure 1 :

La fonction qui construit la table des préfixes est écrite selon la méthode exposée à la section ?? plus bas.

L'algorithme KMP proprement dit

Nous supposerons donc l'existence d'un tableau Tpref qui, pour chaque position k du mot recherché à laquelle une comparaison échoue (et qui correspond donc, si nous avons commencé dans le texte à la position i, à la position i+k du texte), nous donne la position où recommencer la comparaison en évitant les positions qu'il est inutile de vérifier parce que c'est déjà fait. Ou, en d'autres termes, de combien de positions dans le texte il faut revenir en arrière à partir de i+k. Soit :

si nous avons comparé fructueusement de Texte[i] à Texte[i+k-1], mais échoué en comparant Texte[i+k] à Mot[k], alors la comparaison doit pour avoir des chances de réussite reprendre en Texte[i+k-Tpref[k-1]].Notons que Tpref[-1] existe et vaut −1, parce que la comparaison peut échouer au premier caractère du texte, Texte[0], et donc du mot, Mot[0], et alors il faut recommencer en Texte[1].

Nous verrons plus bas la construction de Tpref, pour l'instant supposons ce tableau construit et donné par la procédure kmp:tableau. Voici le principe de l'algorithme :

Algorithme : kmp :KMP Données : Mot, Texte Résultat : la position à laquelle est Mot dans Texte Soient Tpref <- kmp :tableau(Mot) Lmot <- longueur(Mot) Ltexte <- longueur(Texte) i <- 0 j <- 0 tant que j < Lmot et j+i < Ltexte faire si Texte[j+i] = Mot[j] alors j <- j+1 sinon i <- j+i - Tpref[j] si j > 0 alors j <- Tpref[j] finsi finsi fait si j = Lmot alors retourner i sinon retourner -1 fin

En moyenne cet algorithme s'exécute en un temps moyen qui croît comme O(n).

Signalons l'excellent article de Wikipédia, en français.

En Scheme :

(module kmp
 (main appel)
 (import kmp-tableau))
(define (appel args)
 (print (kmp:KMP (cadr args) (caddr args))))
(define (kmp:KMP Mot Texte)
 (let ((Tpref (kmp:tableau Mot))
 (L-texte (string-length Texte))
 (L-mot (string-length Mot))
 (i 0)
 (j 0))
 (let boucle ()
 (if (char=? (string-ref Texte (+ i j))
 (string-ref Mot j))
 (set! j (+ j 1))
 (begin
 (set! i (- (+ i j) (vector-ref Tpref j)))
 (if (> j 0)
 (set! j (vector-ref Tpref j)))))
 (if (and (< (+ i j) L-texte)
 (< j L-mot))
 (boucle)))
 (if (= j L-mot)
 i
 (+ i j))))

Télécharger

Construction de la table des préfixes

Le but de la construction de la table des préfixes est de déterminer, en fonction du résultat d'une comparaison partielle entre une zone du texte et le début du mot cherché, début que nous appelons motif, quelle est la position suivante, dans le mot, à partir de laquelle il est utile de reprendre la comparaison.

Si la comparaison a été positive entre les k caractères de la zone Texte[i..i+k-1] et les k caractères du mot, soit Mot[0..k-1], c'est que ces caratères sont identiques, dirait Monsieur de la Palice. Nous allons « pré-chercher » dans le mot les sous-chaînes, que nous appellerons motifs, telles que Mot[0..j] = Mot[k-j..k]. Munis de la liste de ces motifs, nous saurons combien de caractères sauter lorsqu'une comparaison échouera au caractère de rang k du mot.

Pour chaque position dans le mot nous déterminerons la longueur du plus long motif qui se termine à la position courante.

Tpref vérifiera la propriété suivante : si Texte[i..i+k-1] = Mot[0..k-1] mais que Texte[i+k] != Mot[k], alors la première position à laquelle une occurrence du mot peut apparaître est Texte[i+k-Tpref[k-1]]. En particulier, Tpref[-1] existe et vaut −1, parce que la comparaison peut échouer au premier caractère du texte, et donc du mot.

Tpref[j] aura pour valeur la longueur du motif initial le plus long qui se termine en Mot[j]. Notamment, Tpref[-1] = -1, comme signalé ci-dessus.

Voici l'exposé de l'algorithme. Afin d'éviter un cas particulier, nous adopterons la convention suivante : Mot[-1] existe et sa valeur est différente de toute autre valeur possible pour les caractères de Mot. Nous choisirons ici le caractère nul ASCII, qui s'écrit en Scheme

. Pour des raisons de syntaxe de Scheme, les indices de notre tableau seront augmentés de 1 : Tpref[i] dans le pseudo-code est équivalent à Tpref[i+1] dans les développements ci-dessus.

Algorithme : kmp :tableau Donnée : Mot Soient Lmot <- longueur de Mot Tpref un tableau de longueur Lmot+1 i <- 0 j <- -1 Tpref[0] <- -1 c <- #a000 ;; le caractère nul tant que i < Lmot faire si c = Mot[i] alors Tpref[i+1] <- j+1 j <- j+1 i <- i+1 sinon si j > 0 alors j <- Tpref[j] sinon Tpref[i+1] <- 0 i <- i+1 j <- 0 finsi finsi c <- Mot[j] fait retourner Tpref fin

En Scheme (algorithme adapté de l'article de Wikipedia, en français) :

(module kmp-tableau
 (export (kmp:tableau Motif)))
(define (kmp:tableau Motif)
 (let* ((L-motif (string-length Motif))
 (Tpref (make-vector (+ L-motif 1) 0))
 (i 0)
 (j -1))
 (vector-set! Tpref 0 j)
 (let boucle ((c #a000)) ;; le caractère nul
 (print Tpref)
 (cond ((= i L-motif) 'fini)
 ((char=? c (string-ref Motif i))
 (vector-set! Tpref (+ i 1) (+ j 1))
 (set! j (+ j 1))
 (set! i (+ i 1)))
 ((> j 0)
 (set! j (vector-ref Tpref j)))
 (else
 (vector-set! Tpref (+ i 1) 0)
 (set! i (+ i 1))
 (set! j 0)))
 (if (< i L-motif)
 (boucle (string-ref Motif j))
 Tpref))))

Télécharger

Trace d'exécution

Voici un exemple de déroulement de l'algorithme de construction du tableau des préfixes :

Init. : Motif <- "ATATAGA" Tpref <- 00000000 i <- 0 j <- -1 c <- #a000 Tpref[0] <- -1 ---- c = #a000 != Motif[0] = A, donc : j = -1, donc : Tpref[1] <- 0 i <- 1 j <- 0 c <- A ---- c = A != Motif[1] = T, donc : j = 0, donc : Tpref[2] <- 0 i <- 2 j <- 0 c <- A ---- c = A = Motif[2], donc : Tpref[3] <- 1 ; (j + 1) j <- 1 i <- 3 c <- T ---- c = T = Motif[3], donc : Tpref[4] <- 2 ; (j + 1) j <- 2 i <- 4 c <- A ---- c = A = Motif[4], donc : Tpref[5] <- 3 ; (j + 1) j <- 3 i <- 5 c <- T ---- c = T != Motif[5] = G , donc : j <- Tpref[j = 3] = 1 c <- T ---- c = T != Motif[5] = G, donc : j <- Tpref[j = 1] = O c <- A ---- c = A != Motif[5] = G, donc : j = 0, donc : Tpref[i+1 = 6] <- 0 i <- 6 j <- 0 c <- A ---- c = A = Motif[6] Tpref[7] <- 1 j <- 1 i <- 7 c <- T ---- i = longueur(Motif) => fait Tpref = #(-1 0 0 1 2 3 0 1)

En moyenne cet algorithme s'exécute en un temps moyen qui croît comme O(n).

KMP en style récursif

Aux fins d'une meilleure conformité aux canons de l'élégance lispienne, voici les mêmes programmes réécrits dans un style récursif dépourvu d'affectations :

(module kmp
 (main main)
 (import kmp-table))
(define (main args)
 (print (kmp:KMP (cadr args) (caddr args))))
(define (kmp:KMP Word Text)
 (let ((Tpref (kmp:table Word))
 (Ltext (string-length Text))
 (LastCharPos (- (string-length Word) 1)))
 (let loop ((i 0)
 (j 0))
 (cond ((>= (+ i j) Ltext)
 -1)
 ((char=? (string-ref Text (+ i j))
 (string-ref Word j))
 (if (= j LastCharPos)
 i
 (loop i (+ j 1))))
 (else
 (loop
 (- (+ i j) (vector-ref Tpref j))
 (if (> j 0)
 (vector-ref Tpref j)
 j))) ) ) ))

Télécharger

(module kmp-table
 (export (kmp:table Word)))
(define (kmp:table Word)
 (let* ((WordLength (string-length Word))
 (Tpref (make-vector (+ WordLength 1) 0)) )
 (vector-set! Tpref 0 -1)
 (let loop ((i 0)
 (j -1)
 (c #a000)) ;; null character
 (if (>= i WordLength)
 Tpref
 (cond ((char=? c (string-ref Word i))
 (vector-set! Tpref (+ i 1) (+ j 1))
 (loop (+ i 1)
 (+ j 1)
 (string-ref Word (+ j 1))))
 ((> j 0)
 (let ((j2 (vector-ref Tpref j)))
 (loop i
 j2
 (string-ref Word j2))))
 (else
 (vector-set! Tpref (+ i 1) 0)
 (loop (+ i 1)
 0
 (string-ref Word 0)))) ) )))

Télécharger

This document was translated from L^AT_EX by H^EV^EA.

Algorithme de Needleman et Wunsch

2008-05-27T21:04:00Z

http://www.w3.org/TR/REC-html40/loo...">

L'algorithme de Needleman et Wunsch

Laurent Bloch

Le 10 mai 2006

Sommaire-

Principes de l'algorithme

Dans un autre article de ce site sont présentés des algorithmes de recherche d'un mot dans un texte, notamment celui de Knuth-Morris-Pratt (KMP). Ces algorithmes sont dévolus au problème de la recherche exacte : il s'agit de trouver, si elle existe, la première occurrence exacte de ce mot dans ce texte.

Nous allons maintenant étudier, parce que c'est un problème central en bioinformatique, une recherche approximative : il s'agit de savoir si deux mots se ressemblent, quel est leur degré de ressemblance, ou de trouver, dans un ensemble de mots, celui qui ressemble le plus à un mot-cible. Et nous allons voir que ce problème relève de solutions très différentes de celles qui valent pour la recherche exacte.

Notons d'abord que la ressemblance est une notion imprécise : la plupart des algorithmes utilisés proposent différents paramètres pour ajuster les facteurs de ressemblance aux caractéristiques du problème traité.

Les algorithmes utilisés fournissent en général deux résultats :

pour chaque comparaison de deux chaînes, un score de ressemblance, qui permet ensuite de trouver la meilleure ressemblance parmi un ensemble de comparaisons ;
un alignement des deux chaînes (qui n'ont pas forcément la même longueur) selon la configuration qui procure le meilleur score ; on dit bien un alignement, et non pas l'alignement, parce qu'en effet, comme nous le verrons plus loin, le problème peut admettre plusieurs solutions conduisant au même score.

Le plus caractéristique de cette famille d'algorithmes est peut-être celui de Needleman et Wunsch, que nous étudierons ici ; il réalise un alignement global de deux séquences (chaînes de caractères).

Calculer un alignement global peut être coûteux si les séquences à aligner sont longues, ou s'il y en a beaucoup. D'autres algorithmes, qui ressemblent à celui-ci, ont été conçus pour limiter la taille du problème en ne réalisant l'alignement que pour des régions « intéressantes ». La détermination des régions intéressantes est bien sûr en elle-même un problème intéressant. Citons l'algorithme de Smith et Waterman, qui réalise des alignements locaux, et le logiciel BLAST¹, qui mettent en oeuvre des méthodes similaires à celles de Needleman et Wunsch, après des optimisations éventuellement complexes.

Le problème de la comparaison de séquences est exponentiel, la solution est en O(kⁿ) ; ces algorithmes sont susceptibles d'une multiplicité de solutions ; une des techniques les plus généralement utilisées pour en réduire la complexité est la programmation dynamique, qui fait l'objet d'un autre article sur ce site.

La programmation dynamique résout des problèmes en combinant des solutions de sous-problèmes. (Thomas Cormen, Charles Leiserson, Ronald Rivest et Clifford Stein, Introduction à l'algorithmique)

L'idée de la programmation dynamique est de mémoriser les résultats de calculs intermédiaires qui seront probablement répétés. La programmation dynamique est par exemple souvent un bon choix lorsque l'on aura besoin, après les avoir calculées, des valeurs stockées dans tous les noeuds d'un arbre ou dans toutes les cases d'un tableau. Parfois aussi cette conservation des résultats intermédiaires est imposée par un problème tel que le calcul d'une valeur se fait en fonction de toutes les précédentes. L'art algorithmique consiste à chercher des solutions qui évitent ce type de contrainte mais c'est parfois impossible. Et puis il y a des problèmes intrinsèquement récursifs pour lesquels n'existe pas d'algorithme itératif.

Nous allons donc chercher des procédés pour associer un algorithme qui calcule des valeurs successives avec une structure de données qui les archive.

Cette section doit beaucoup au travail préalable de William Saurin pour cet enseignement, ainsi qu'à une page créée par Eric C. Rouchka, de l'université Washington à Saint-Louis, et reprise par Per Kraulis au Stockholm Bioinformatics Center :

http://www.sbc.su.se/~pjk/molbioinf...

Supposons que nous souhaitions calculer un alignement global de deux séquences :

séquence n° 1 : G A A T T C A G T T A

séquence n° 2 : G G A T C G A

La séquence n° 1 a m=11 résidus, la séquence n° 2 n=7 résidus.

Nous allons ici étudier l'algorithme avec des paramètres particulièrement simples, peut-être même simplistes : pénalité nulle pour les trous (gaps) et les discordances (substitutions), une pénalité négative, ou prime, égale à 1 pour les concordances (matches). Le but est d'acquérir une vue d'ensemble de l'architecture de la solution, qui permettra au lecteur d'envisager ensuite des exemples plus compliqués, avec des formules de calcul plus élaborées pour les scores et pour les pénalités de gap.

Le principe de pondération que nous adopterons sera le suivant :

la « prime de score » pour la comparaison du résidu de rang i de la première séquence avec le résidu de rang j de le seconde séquence sera S_i,j = 1 si les deux résidus sont identiques, sinon :
S_i,j = 0 (score de discordance) ;
w = 0 (pénalité de gap).

L'algorithme opère en trois étapes :

initialisation ;
calcul des scores et remplissage de la matrice ;
calcul de l'alignement en « remontant » dans la matrice.

Initialisation

Création d'une matrice M de m+2=13 colonnes et n+2=9 lignes : la ligne et la colonne de rangs 0 contiendront les textes des séquences, la seconde ligne (de rang 1, les M_1,j) et la première colonne (les M_i,1) de M sont remplies de 0 parce que nous avons posé qu'il n'y avait pas de pénalités pour des gaps initiaux ou finals.

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0
G	0
A	0
T	0
C	0
G	0
A	0

Remplissage de la matrice

À chaque position M_i,j de la matrice M (i est le numéro de ligne, j le numéro de colonne) le score se calcule ainsi :

M_i,j = Maximum de :

M_i−1,j−1 + S_i,j	(concordance ou discordance
	dans la diagonale)
M_i,j−1 + w	(gap dans la séquence n° 1)
M_i−1,j + w	(gap dans la séquence n° 2)

Nous voyons que pour calculer M_i,j il faut (et il suffit de) connaître M_i−1,j, M_i,j−1 et M_i−1,j−1 ; de ce point de vue le problème est assez analogue à ceux posés par Fibonacci ou par le triangle de Pascal.

Ainsi, comme chaque séquence commence par le résidu G (concordance), S_1,1=1. Nous avons posé par hypothèse w=0. Donc :

M_1,1	= Max[M_0,0+1, M_1,0+0, M_0,1+0]	(1)
	= Max[1,0,0]	(2)
	= 1	(3)

Nous pouvons donc inscrire un 1 en M_1,1 :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1
G	0
A	0
T	0
C	0
G	0
A	0

Ceci fait, toujours parce que w=0, nous pouvons facilement remplir la ligne 1 et la colonne 1 avec des 1 ; ainsi :

M_2,1	= Max[M_1,0+0, M_2,0+0, M_1,1+0]	(4)
	= Max[0,0,1]	(5)
	= 1	(6)

soit :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1
A	0	1
T	0	1
C	0	1
G	0	1
A	0	1

Finalement :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Nous avons signalé ci-dessus que le problème général de la comparaison de séquences était exponentiel (O(kⁿ)). L'utilisation de la programmation dynamique, avec le graphe représenté par ce tableau, permet de le réduire à un problème quadratique (O(m × n), m et n étant les longueurs respectives des séquences). En effet, il y a m × n valeurs dans la table, et le calcul de chacune s'effectue en temps constant.

Déterminer l'alignement optimal

L'étape précédente nous a déjà permis de savoir que le score d'alignement maximum pour nos deux séquences est 6. Souvent, cette information est suffisante, parce que l'on cherche en fait les meilleurs scores parmi une collection de séquences à comparer à la cible. Mais peut être aussi intéressant de connaître un alignement qui donne ce score.

Nous allons maintenant déterminer l'alignement effectif qui donne ce résultat.

Pour cela, on considère la case du tableau qui contient le score maximum, qui est M_m,n, et on la compare à ses voisines. Ici toutes les voisines contiennent la valeur 5. Comme la différence de scores est 1 dans tous les cas, et que le seul moyen d'avoir un accroissement de 1 est une concordance (match) (toutes les autres situations donnent un accroissement nul), c'est que la case précédente était la voisine en diagonale :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Ce qui nous donne un alignement :

Maintenant nous considérons la case courante et cherchons celle qui la précède : c'est la voisine avec le score maximum, soit celle de la même ligne.

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Cet alignement correspond à un gap dans la séquence n° 2 :

T	A
	\|
_	A

Encore une fois, le prédécesseur immédiat donne un gap dans la séquence n° 2 :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

T	T	A
		\|
_	_	A

Au bout du compte :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

G	A	A	T	T	C	A	G	T	T	A
\|		\|		\|	\|		\|			\|
G	G	A	_	T	C	_	G	_	_	A

Il y a une autre solution possible :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

qui donne l'alignement suivant :

G	_	A	A	T	T	C	A	G	T	T	A
\|			\|		\|	\|		\|			\|
G	G	_	A	_	T	C	_	G	_	_	A

L'algorithme

Algo : NW Données : s1, s2, S et gap ; des chaînes numérotées de ; 1 à longueur(s1) et de 1 à longueur(s2), ; S le score de concordance entre caractères ; et gap un coût de gap. Résultat : la matrice de calcul des scores Créer C une matrice à longueur(s1) + 2 colonnes et à longueur(s2) + 2 lignes pour j allant de 2 à longueur(s1)+2 faire C[1, j] <- gap * j-2 fait pour i allant de 2 à longueur(s2)+2 faire C[i, 1] <- gap * i-2 pour j allant de 2 à longueur(s2)+2 faire C[i, j] = max(C[i-1,j-1] + si s1[i] = s2[j] alors 1 sinon 0, C[i-1,j] + gap, C[i,j-1] + gap) fait fait retourner C

Nous verrons dans un prochain article l'algorithme de remontée dans le graphe (backtracking) pour trouver un alignement optimal.

Au sujet de ces algorithmes on consultera avec profit le livre de Maxime Crochemore, Christophe Hancart et Thierry Lecroq, Algorithmique du texte, chez Vuibert.

Le programme

(module nw:lb (main main) (import nw:matrices) (import nw:chains) (import nw:alignment)) (define (nw-2 s1 s2 match-bonus gap-penalty) (let ((ncol (+ (chain-length s1) 2)) (nlin (+ (chain-length s2) 2))) (let ((C (make-matrix nlin ncol 0))) (matrix:margins C s1 s2) ;; (do ((j 2 (+ j 1))) ((= j nlin) 'fait) (matrix:set! C j 1 (* j gap-penalty))) (do ((i 2 (+ i 1))) ((= i ncol) 'fait) (matrix:set! C 1 i (* gap-penalty i))) (do ((i 2 (+ i 1))) ((= i ncol) C) (do ((j 2 (+ j 1))) ((= j nlin) 'fait) (let ((val (max (+ (matrix:ref C (- j 1) (- i 1)) (if (char=? (matrix:ref C 0 i) (matrix:ref C j 0)) match-bonus 0)) (+ (matrix:ref C j (- i 1)) gap-penalty) (+ (matrix:ref C (- j 1) i) gap-penalty)))) (matrix:set! C j i val))))))) ;; là on suppose qu'une séquence est dans un fichier fasta ;; read-fasta lit un fichier fasta et rend la premiere séquence trouvée (define (read-fasta port) (let ((titre (read-line port))) (if (or (eof-object? titre) (zero? (string-length titre)) (not (char=? (string-ref titre 0) #\>))) (error 'read-fasta "not a fasta file" port) (let loop ((str "")) (let ((lu (read-line port))) (if (or (eof-object? lu) (char=? #\> (string-ref lu 0))) (cons titre str) (begin (print lu) (loop (string-append str lu))))))))) (define (f-concorde base1 base2 match-bonus) (if (char=? base1 base2) match-bonus 0)) ;; (define (usage) (print "nw fichier-1 fichier-2 match-bonus gap-penalty") (exit 1)) (define (main argv) (if (not (= (length argv) 5)) (usage) (let ((f1 (cadr argv)) (f2 (caddr argv)) (match-bonus (string->number (cadddr argv))) (gap-penalty (string->number (cadddr (cdr argv))))) (let* ((s1 (make-chain (cdr (let ((port (open-input-file f1))) (read-fasta port))))) (s2 (make-chain (cdr (let ((port (open-input-file f2))) (read-fasta port))))) (the-score-matrix (nw-2 s1 s2 match-bonus gap-penalty))) (matrix:print the-score-matrix) (matrix:print (alignment the-score-matrix match-bonus gap-penalty))))))

Le module de manipulation de matrices

Voici les deux séquences au format FASTA :

> La séquence 1GAATTCAGTTA

Chaque fichier ne peut comporter qu'une séquence, mais le texte de la séquence peut s'étendre sur plusieurs lignes. Vous pouvez visiter au NCBI le site de référence du format FASTA.

> La séquence 2GGATCGA

Le module de matrices :

(module nw:matrices (export (make-matrix n m . fill) (matrix? obj) (matrix:ref T i j) (matrix:set! T i j val) (matrix:nlines T) (matrix:ncols T) (matrix:margins M s1 s2) (matrix:print T)) (import nw:chains)) ;; il nous faut des matrices (define matrix:tag "*MATRIX*") (define (make-matrix lin col . fill) (let ((the-table (vector "*MATRIX*" (make-vector lin #f)))) (do ((i 0 (+ i 1))) ((= i lin)) (vector-set! (vector-ref the-table 1) i (if (null? fill) (make-vector col) (make-vector col (car fill))))) the-table)) ;; un prédicat d'appartenance, pour vérifier qu'un ;; objet appartient bien au type : (define (matrix? obj) (and (vector? obj) (string=? (vector-ref obj 0) "*MATRIX*") (vector? (vector-ref obj 1)))) ;; un mutateur, pour modifier un objet du type en ;; affectant une valeur à un élément du tableau : (define (matrix:set! T i j val) (if (matrix? T) (vector-set! (vector-ref (vector-ref T 1) i) j val))) ;; un sélecteur, pour accéder à un élément d'un ;; tableau du type : (define (matrix:ref T i j) (if (matrix? T) (vector-ref (vector-ref (vector-ref T 1) i) j))) (define (matrix:margins M s1 s2) (let ((nlin (matrix:nlines M)) (ncol (matrix:ncols M))) (do ((j 2 (+ j 1)) (c (chain-ref s1 1) (chain-ref s1 (min j (- ncol 2))))) ((= j ncol) 'fait) (matrix:set! M 0 j c)) (do ((i 2 (+ i 1)) (c (chain-ref s2 1) (chain-ref s2 (min i (- nlin 2))))) ((= i nlin) 'fait) (matrix:set! M i 0 c)) (matrix:set! M 0 0 #\space) (matrix:set! M 0 1 #\space) (matrix:set! M 1 0 #\space))) ;; diverses procédures utilitaires dont la fonction se ;; comprend d'elle-même : (define (matrix:nlines T) (if (matrix? T) (vector-length (vector-ref T 1)))) (define (matrix:ncols T) (if (matrix? T) (vector-length (vector-ref (vector-ref T 1) 0)))) (define (matrix:print T) (if (matrix? T) (let ((n (matrix:nlines T)) (m (matrix:ncols T))) (do ((i 0 (+ 1 i))) ((= i n)) (let ((this-line (vector-ref (vector-ref T 1) i))) (do ((j 0 (+ 1 j))) ((= j m)) (display (vector-ref this-line j)) (display " ")) (newline))))))

Le module de chaînes :

(module nw:chains (export (make-chain s) (chain-ref s i) (chain-set! s i c) (chain-length s))) ;; il nous faut des chaîne numérotées de 1 à longueur de chaîne: (define (make-chain s) ; prend une string et rend un chaîne à partir de 1 s) (define (chain-ref s i) (string-ref s (- i 1))) (define (chain-set! s i c) (string-set! s (- i 1) c)) (define (chain-length s) (string-length s))

Le Makefile :

BIGLOO = bigloo AFILE = NW-afile.scm BGL_FLAGS = -afile $(AFILE) -Obench -farithmetic -static-bigloo CIBLE = nw REPERT_CIBLE = . %.o: %.scm @ $(BIGLOO) $(BGL_FLAGS) -c $*.scm -o $*.o OBJECTS = NW-lb.o NW-matrices.o NW-chains.o NW-alignment.o SOURCES = NW-lb.scm NW-matrices.scm NW-chains.scm NW-alignment.scm all: $(REPERT_CIBLE)/$(CIBLE) $(REPERT_CIBLE)/$(CIBLE): $(OBJECTS) @ echo "Edition de liens..." @ $(BIGLOO) $(BGL_FLAGS) $(OBJECTS) \ -o $(REPERT_CIBLE)/$(CIBLE) @ echo "$(REPERT_CIBLE)/$(CIBLE) construit." @ echo "-------------------------------" $(REPERT_CIBLE): @ mkdir -p $(REPERT_CIBLE) install: cp $(REPERT_CIBLE)/$(CIBLE) $(REPERT_CGI) clean: -rm -f $(OBJECTS) $(SOURCES_C) -rm -f *~ Src/*~ Src/*.o Src/*.mco @ echo "nettoyage fait..." @ echo "-------------------------------" # destruction aussi des binaires : cleanall: clean -rm -f $(REPERT_CIBLE)/$(CIBLE)

1: BLAST (Basic Local Alignment Search Tool) est un logiciel écrit par S.F. Altschul, W. Gish, W. Miller, Gene Myers et D.J. Lipman. Il est l'outil de travail quotidien de tous ceux qui travaillent en biologie moléculaire (cf. http://en.wikipedia.org/wiki/BLAST).

This document was translated from L^AT_EX by H^EV^EA.

Needleman et Wunsch, once again

2008-05-27T16:09:00Z

http://www.w3.org/TR/REC-html40/loo...">

Needleman et Wunsch : le retour

Laurent Bloch

Le 18 mai 2006

Calculer un alignement d'après un graphe

Cet article est la suite de l'Algorithme de Needleman et Wunsch, où nous avions montré comment déterminer le score maximum, c'est-à-dire le coût minimum, d'alignement de deux séquences. Nous avions eu recours à une méthode de programmation dynamique, par la construction d'une matrice des scores. Ceci fait, restait à exhiber un alignement optimal : c'est l'objet du présent article. Nous obtiendrons ce résultat en remontant dans le graphe représenté par la matrice, selon un chemin qui va nous donner l'alignement. Reprenons la démarche de l'article.

Pour cela, on part de la case du tableau qui contient le score maximum, qui est M_m,n, et on considère ses voisines. Ici toutes les voisines contiennent la valeur 5. Comme la différence de scores est 1 dans tous les cas, et que le seul moyen d'avoir un accroissement de 1 est une concordance (toutes les autres situations donnent un accroissement nul), c'est que la case précédente était la voisine en diagonale :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Ce qui nous donne un alignement :

Maintenant nous considérons la case courante et cherchons celle qui la précède : c'est la voisine avec le score maximum, soit celle de la même ligne.

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Cet alignement correspond à un gap dans la séquence n° 2 :

T	A
	\|
_	A

Encore une fois, le prédécesseur immédiat donne un gap dans la séquence n° 2 :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

T	T	A
		\|
_	_	A

Au bout du compte :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

G	A	A	T	T	C	A	G	T	T	A
\|		\|		\|	\|		\|			\|
G	G	A	_	T	C	_	G	_	_	A

Il y a une autre solution possible :

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

qui donne l'alignement suivant :

G	_	A	A	T	T	C	A	G	T	T	A
\|			\|		\|	\|		\|			\|
G	G	_	A	_	T	C	_	G	_	_	A

Au sujet de cet algorithme on consultera avec profit le livre de Maxime Crochemore, Christophe Hancart et Thierry Lecroq, Algorithmique du texte, chez Vuibert.

Le programme

(module nw:alignment (export (alignment C match-bonus gap-penalty)) (import nw:matrices)) (define (alignment-col-set! align col . vals) (matrix:set! align 0 col (car vals)) (matrix:set! align 1 col (cadr vals)) (matrix:set! align 2 col (caddr vals))) (define (alignment C match-bonus gap-penalty) (let* ((ncol (matrix:ncols C)) (nlin (matrix:nlines C)) (line-length (+ nlin ncol (- 4))) (this-alignment (make-matrix 3 line-length #\space))) (let boucle ((i (- nlin 1)) (j (- ncol 1)) (p-align (- (+ nlin ncol) 5))) (if (and (= i 1) (= j 1)) 'fait (cond ((and (= j 1) (> i 1)) (alignment-col-set! this-alignment p-align #\_ #\space (matrix:ref C i 0)) (boucle (- i 1) j (- p-align 1))) ((and (= i 1) (> j 1)) (alignment-col-set! this-alignment p-align (matrix:ref C 0 j) #\space #\_) (boucle i (- j 1) (- p-align 1)) ) ((and (= (matrix:ref C i j) (+ (matrix:ref C (- i 1) (- j 1)) match-bonus)) (char=? (matrix:ref C i 0) (matrix:ref C 0 j))) (alignment-col-set! this-alignment p-align (matrix:ref C 0 j) #\| (matrix:ref C i 0)) (boucle (- i 1) (- j 1) (- p-align 1))) ((= (matrix:ref C i j) (matrix:ref C (- i 1) (- j 1))) (alignment-col-set! this-alignment p-align (matrix:ref C 0 j) #\x (matrix:ref C i 0)) (boucle (- i 1) (- j 1) (- p-align 1)) ) ((= (matrix:ref C i j) (+ (matrix:ref C (- i 1) j) gap-penalty)) (alignment-col-set! this-alignment p-align #\_ #\space (matrix:ref C i 0)) (boucle (- i 1) j (- p-align 1))) ((= (matrix:ref C i j) (+ (matrix:ref C i (- j 1)) gap-penalty)) (alignment-col-set! this-alignment p-align (matrix:ref C 0 j) #\space #\_) (boucle i (- j 1) (- p-align 1))) )) this-alignment)))

This document was translated from L^AT_EX by H^EV^EA.

Knuth-Morris-Pratt en style récursif

2008-05-13T20:54:41Z

Voici comme promis une version en style récursif du programme Knuth-Morris-Pratt. Ce type d'exercice est plein de pièges : si vous découvrez que je suis tombé dans l'un d'entre eux, je serais content que vous me le signaliez.

(module kmp-table
 (export (kmp:table Word)))
(define (kmp:table Word)
 (let* ((WordLength (string-length Word))
 (Tpref (make-vector (+ WordLength 1) 0)) )
 (vector-set! Tpref 0 -1)
 (let loop ((i 0)
 (j -1)
 (c #a000)) ;; null character
 (if (>= i WordLength)
 Tpref
 (cond ((char=? c (string-ref Word i))
 (vector-set! Tpref (+ i 1) (+ j 1))
 (loop (+ i 1)
 (+ j 1)
 (string-ref Word (+ j 1))))
 ((> j 0)
 (let ((j2 (vector-ref Tpref j)))
 (loop i
 j2
 (string-ref Word j2))))
 (else
 (vector-set! Tpref (+ i 1) 0)
 (loop (+ i 1)
 0
 (string-ref Word 0)))) ) )))

Télécharger

(module kmp
 (main main)
 (import kmp-table))
(define (main args)
 (print (kmp:KMP (cadr args) (caddr args))))
(define (kmp:KMP Word Text)
 (let ((Tpref (kmp:table Word))
 (L-texte (string-length Text))
 (LastCharPos (- (string-length Word) 1)))
 (let loop ((m 0) ;; match
 (i 0)) ;; index
 (cond ((>= (+ m i) L-texte)
 -1)
 ((char=? (string-ref Text (+ m i))
 (string-ref Word i))
 (if (= i LastCharPos)
 m
 (loop m (+ i 1))))
 (else
 (loop
 (- (+ m i) (vector-ref Tpref i))
 (if (> i 0)
 (vector-ref Tpref i)
 i))) ) ) ))

Télécharger

Corrigé de l'examen BNF-103 de juin 2006

2007-06-05T14:06:41Z

http://www.w3.org/TR/REC-html40/loo...">

correction-exam-BNF103-200606

Conservatoire National des Arts et Métiers

292 rue Saint-Martin – 75141 PARIS Cedex 03

–0–

Chaire de Bioinformatique

Corrigé de l'examen du :
lundi 26 juin 2006 Heure : 18h15 – 20h15

TITRE DE L'ENSEIGNEMENT : Algorithmique de la Bioinformatique
(BNF103)

1.

Programme min

L'algorithme suivant permet de calculer l'indice de la valeur minimale d'un tableau.

1 Nom de l'algorithme : Min 2 Entrée : A; un tableau indice de 0 à longueur(A)-1 3 Sortie : l'indice de la plus petite valeur de A 4 imin <- 0; 5 pour i allant de 1 à longueur(A) -1 faire 6 si A[imin] > A[i] alors 7 imin <- i 8 finsi 9 fait 10 retourner(imin)

Ecrire une fonction Scheme qui implémente cet algorithme.
Combien de fois le test de la ligne 6 sera-t-il exécuté ?
Combien de fois l'affectation de la ligne 7 sera-t-elle exécutée dans le pire des cas, dans le meilleur des cas, en moyenne ?
Quand la ligne 7 est exécutée que peut-on dire de la relation entre A[i] et les valeurs des cases précédentes du tableau ?

Réponse :

Fonction Scheme :

(define (minr A) (let ((imin 0)) (do ((i 1 (+ i 1))) ((>= i (vector-length A)) imin) (if (> (vector-ref A imin) (vector-ref A i)) (set! imin i)))))
Soit L la longueur du tableau A. Le test de la ligne 6 sera effectué L-1 fois.
Si le tableau est trié en ordre décroissant, le nombre contenu dans chaque case sera inférieur au nombre contenu dans la case précédente, et ainsi l'affectation de la ligne 7 sera effectuée L-1 fois (le pire des cas).
Si le tableau est trié en ordre croissant, l'élément de la case 0 sera le plus petit, et l'affectation de la ligne 7 ne sera jamais effectuée (le meilleur des cas).

En moyenne, l'affectation de la ligne 7 sera effectuée integer(L−1/2) fois.
La ligne 7 est exécutée si A[i] est inférieur à A[imin], ce qui signifie que A[i] est inférieur aux valeurs de toutes les cases précédentes du tableau.

2.

Tri par sélection

Nous voulons définir un algorithme qui permette de trier un tableau A dont les indices vont de 0 à longueur(A)-1. Pour cela nous allons procéder par étapes. La première sera numérotée 0, la seconde 1, etc...

Au cours d'une étape donnée i nous considérerons que les éléments du tableau jusqu'à l'élément numéro i−1 sont triés et nous chercherons dans les éléments de numéro i, i+1, ... longueur(A)−1 le plus petit, nous échangerons alors cet élément avec l'élément de numéro i puis nous recommencerons.

On considère maintenant la fonction Scheme minr telle qu'un appel à (minr A i) permet de trouver l'indice du plus petit élément de A entre i et longueur(A)-1.

Cette fonction s'écrit :

Cette fonction implémente un algorithme que nous appelons MinR. C'est une variante de l'algorithme Min de la question 1. Ecrire le pseudocode de l'algorithme MinR.

La fonction Scheme suivante réalise l'algorithme du tri par sélection:

Écrire maintenant le pseudocode de l'algorithme de tri par sélection.

Combien de fois l'algorithme de tri par sélection fait-il appel à l'algorithme MinR ?

Réponse :

Pseudo-code de la fonction qui détermine le plus petit élément du tableau A de longeur L entre les indices i et L-1 :

1 Nom de l'algorithme : MinR 2 Entrée : A; un tableau indicé de 0 à L-1 3 i; un indice compris entre 0 et L-1 4 Sortie : l'indice de la plus petite valeur de A 5 entre 0 et L-1 6 imin <- i; 7 pour j allant de i+1 à L-1 faire 8 si A[imin] > A[j] alors 9 imin <- j 10 finsi 11 fait 12 retourner(imin)

Pseudo-code du tri par sélection :
1 Nom de l'algorithme : Tri-select 2 Entrée : A; un tableau indicé de 0 à L-1 3 pour i allant de 0 à L-1 4 j <- MinR(A, i) 5 tmp <- A[j] 6 A[j] <- A[i] 7 A[i] <- tmp 8 fait 9 retourner A
L'algorithme de tri par sélection fait appel L fois à la fonction MinR. L'algorithme de MinR traite à chaque appel un sous-tableau de longueur K variant entre L-1 et 1, soit en moyenne L−2/2, et, comme nous l'avons vu, effectue K−1/2 comparaison, soit L−3/4. L'algorithme de tri effectue donc en moyenne L × L−3/4 opérations élémentaires, il est en O(n²).

3.

Une chaîne d'ADN peut être représentée comme une chaîne de caractères contenant uniquement les caractères A,T,G,C.

On peut aussi représenter une chaîne d'ADN par un nombre en considérant que chaque lettre représente un chiffre et que le nombre représentant la chaîne d'ADN est écrit en base 4. Par exemple en admettant que A représente le chiffre 0, T, le chiffre 1, G, le chiffre 2 et C le chiffre 3, la chaîne d'ADN ATG serait représentée par le nombre 6 = 0 * 16 + 1 * 4 + 2 * 1. Par quel nombre serait représentée la chaîne GTG? la chaîne TAA ?

Réponse :

GTG → 2 × 16 + 1 × 4 + 2 × 1 = 38

TAA → 1 × 16 + 0 × 4 + 0 × 1 = 16

Si une chaîne d'ADN a une longueur de 2 nucléotides (deux lettres), quelle est la chaîne représentée par le nombre 15 ? par le nombre 0 ?

Réponse :

15 = 3 × 4 + 3 × 1 = CC

00 = 0 × 4 + 0 × 1 = AA

Quelle est le plus grand nombre qui représente une chaîne d'ADN de longueur 10 ? Quelle chaîne représente-t-il ? Le nombre représentant une chaîne d'ADN dépend de l'association que nous définissons entre les lettres A,T,G et C et les chiffre 0, 1, 2, 3. Si nous décidons d'associer T à 0, G à 1, C à 2 et A à 3, par quel nombre serait représentée la chaîne ATG ? la chaîne GTG ? Quelle serait la chaîne d'ADN représentée par le nombre 15 ?

Réponse :

Dans notre système à base 4, le « chiffre » qui a la plus grande valeur est C, qui vaut 4−1=3. Le plus grand nombre représenté par dix chiffres en base 4 sera constitué de dix « chiffres » C, et vaudra :

(4−1) × 4⁹ + (4−1) × 4⁸ + ... + (4−1) × 4² + (4−1) × 4¹ + (4−1) × 4⁰

soit en développant :

4¹⁰ − 4⁹ + 4⁹ − 4⁸ + 4⁸ − ... + 4³ − 4² + 4² − 4¹ + 4¹ −1

et grâce à ces combinaisons d'additions et de soustractions :

CCCCCCCCCC → 4¹⁰−1 = 1048575

Combien y a-t-il de façons d'associer les lettres A, T, G et C chacune avec un chiffre différent pris entre 0 et 3 ?

Réponse :

S'il y a 4 choix possibles pour A, il en reste 3 pour T, 2 pour G et 1 pour C : donc 4!=24 solutions possibles, qui sont appelées les permutations de l'ensemble { A,T,G,C } (cf. http://fr.wikipedia.org/wiki/Permutation).

Une chaîne de protéine est une chaîne d'acides aminés, elle peut être représentée comme une chaîne de caractères ne contenant que les caractères suivants:

A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y.

En effet il y a 20 acides aminés. On peut représenter une chaîne de protéine par un nombre en considérant que chacun des caractères pouvant apparaître dans une chaîne de protéine est associé à un nombre pris entre 0 et 19. Le calcul pour une protéine se fera donc en base 20.

On peut prendre comme correspondance entre acide aminé et chiffre la correspondance suivante par exemple:

A: 0; C: 1; D: 2; E: 3; F: 4; G: 5; H: 6; I: 7; K: 8; L: 9; M:10; N:11; P:12; Q:13; R:14; S:15; T:16; V:17; W:18; Y:19;

Dans ce cas quelle est la valeur du nombre représentant le tripeptide CAY ?

Réponse :

CAY → 1 × 20² + 0 × 20¹ + 19 × 20⁰

soit : 1 × 400 + 0 × 20 + 19 × 1 = 419

Quel est le plus grand nombre représentant un tripeptide (une protéine formée de trois acides aminés) ? Quel est le plus grand nombre représentant une protéine de 300 acides aminés ?

Réponse :

Le raisonnement est le même qu'à la question :

YYY → 20³−1 = 7999

YYYYYYYYYYYYYYYYYYYY → 20³⁰⁰−1

Soit la fonction Scheme:

Cette fonction rend la valeur d'un nombre associé à une chaîne d'ADN. Quelle correspondance a-t-on choisi entre les nucléotides (lettres) de la chaîne d'ADN et les chiffres 0,1,2,3 ?

Réponse :

T → 0

G → 1

C → 2

A → 3

Donner le pseudocode de l'algorithme décrit par la fonction adn->number.

Réponse :

1 Nom de l'algorithme : ADN->number 2 Entrée : A ; une séquence d'ADN représentée par une 3 chaîne de caractères 4 Résultat : la valeur numérique qui représente la séquence 5 valeur <- 0 6 L <- longueur de A 7 pour i allant de 0 à L-1 faire 8 si A[i] = T alors valeur <- valeur + 0 9 sinon si A[i] = G alors valeur <- valeur + 1 10 sinon si A[i] = C alors valeur <- valeur + 2 11 sinon si A[i] = A alors valeur <- valeur + 3 12 finsi 13 fait 14 retourner valeur

On peut se demander si deux chaînes sont identiques de la façon suivante: ont-elle la même longueur et dans ce cas sont-elles représentées par le même nombre ? Ecrire la fonction.

Réponse :

Ce document a été traduit de L^AT_EX par H^EV^EA

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

Site WWW de Laurent Bloch

Plan du cours d'algorithmes pour la biologie (BNF 103)

Recherche de mots dans un texte (suite)

Recherche d'un motif dans un texte

Algorithmes sur les graphes

Recherche dans des banques de données

Chaines de Markov

Des problèmes difficiles : NP-complétude

Solution presque optimales

Algorithmes de tri

#Sommaire-

Recherche dichotomique

Aurait-on-pu aller plus vite dans la recherche d'un élément dans une liste ?

Algorithme dichotomique

Exemple

Combien de fois fait on le test “grand != petit et t[pivot] != v”

En Scheme :

Tri par insertion

Pseudo-code

Traduisons en Scheme :

Complexité

Tri fusion

Combien de cons ?

Tri rapide

En Scheme :

Tri linéaire

En Scheme :

Tas et tri par tas

Recherche de mots dans un texte, KMP

L'algorithme de Knuth-Morris-Pratt

William Saurin, Laurent Bloch

Le 8 mai 2009

#Sommaire-

Algorithme force brute.

Exposé de l'algorithme KMP

L'algorithme KMP proprement dit

Construction de la table des préfixes

Trace d'exécution

KMP en style récursif

Algorithme de Needleman et Wunsch

L'algorithme de Needleman et Wunsch

Laurent Bloch

Le 10 mai 2006

Sommaire-

Principes de l'algorithme

Initialisation

Remplissage de la matrice

Déterminer l'alignement optimal

L'algorithme

Le programme

Le module de manipulation de matrices

Needleman et Wunsch, once again

Needleman et Wunsch : le retour

Laurent Bloch

Le 18 mai 2006

Calculer un alignement d'après un graphe

Le programme

Knuth-Morris-Pratt en style récursif

Corrigé de l'examen BNF-103 de juin 2006

1.

Réponse :

2.

Réponse :

3.

Réponse :

Réponse :

Réponse :

Réponse :

Réponse :

Réponse :

Réponse :

Réponse :

Réponse :

Combien de `cons` ?

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6

		G	A	A	T	T	C	A	G	T	T	A
	0	0	0	0	0	0	0	0	0	0	0	0
G	0	1	1	1	1	1	1	1	1	1	1	1
G	0	1	1	1	1	1	1	1	2	2	2	2
A	0	1	2	2	2	2	2	2	2	2	2	3
T	0	1	2	2	3	3	3	3	3	3	3	3
C	0	1	2	2	3	3	4	4	4	4	4	4
G	0	1	2	2	3	3	4	4	5	5	5	5
A	0	1	2	3	3	3	4	5	5	5	5	6