Recherche de mots dans une séquence d’ADN

L’objet de cette première partie est de calculer, pour une chaîne d’ADN donnée, son spectre. Expérimentez ici pour voir à quoi il ressemble :

Longueur :

Pour la séquence , et la longueur 3, le spectre est :

{AGG,GGT,GTC,TCA,CAG}

1. Pour cette question et les suivantes, on va utiliser le fait qu’une chaîne de caractères est une liste (de caractères). On peut donc appliquer à la chaîne S de longueur l la notation du tranchage. Les mots de longueur l dans S sont donc

S[0:l]
S[1:l+1]
S[2:l+2]
.
.
.
S[N-1:l+N-1]

en notant N le nombre cherché. Or l’expression l+N-1 désigne la longueur n de S, donc N est solution de l’équation l+N-1==n soit N==n+1-l

2. La liste des mots de longueur l peut donc être obtenue par cet algorithme (ici rédigé en Python), consistant à boucler sur les indices i allant de 0 (inclus) à n-l+1 (exclu) :

def liste_des_mots(S,l):
    return [S[i:i+l] for i in range(len(S)-l+1)]

def mots_longueur_l(sn, l):
    l_mots = []
    for indice in range(len(sn)-l+1):
        mot = ''
        for delta in range(l):
            mot += sn[indice + delta]
        l_mots.append(mot)
    return l_mots

3. On boucle n-l+1 fois, et à chaque passage dans la boucle, on effectue l fois le couple (append,affectation). On effectue donc en tout 2×l×(n-l+1) opérations : l’algorithme est linéaire par rapport à n et par rapport à l.

4. Dans une séquence de longueur n il y a, on l’a vu au 1, exactement n-l+1 mots de longueur l. Mais l’exemple du 3 montre qu’ils ne sont pas nécessairement distincts. En fait, comme il y a 4 lettres dans l’alphabet $\{A,C,G,T\}$, il ne peut y avoir plus de 4^l mots de longueur 4 dans une chaîne écrite sur cet alphabet. Une borne supérieure pour le nombre de mots distincts de longueur l est donc max(n-l+1,4^l).

5. Pour éviter ces répétitions, la théorie des ensembles vient à la rescousse : une structure de données particulièrement adaptée est donc l’ensemble (set en Python). Ce qui donne cette amélioration du script précédent, en remplaçant les crochets des listes, par les accolades des ensembles :

def spectre(S,l):
    return {S[i:i+l] for i in range(len(S)-l+1)}

6. Et il suffit de convertir l’ensemble en liste pour obtenir la liste des mots distincts :

def liste_mots(S,l):
    return list({S[i:i+l] for i in range(len(S)-l+1)})

Remarque : L’énoncé demandait d’ajouter comme paramètre, le nombre n. On ne l’a pas fait ici, parce que la fonction len de Python permet de le récupérer à partir du premier paramètre qui est, justement, de longueur n.

def liste_mots(sn, l):
    l_mots = []
    for indice in range(len(sn)-l+1):
        mot = ''
        for delta in range(l):
            mot += sn[indice + delta]
        if mot not in l_mots:
            l_mots.append(mot)
    return l_mots

Reconstruction d’une séquence ADN

Le problème de la reconstruction de la séquence d’ADN consiste, à partir de la connaissance d’un spectre [...], à déterminer une séquence d’ADN compatible.

7. Copier-coller ci-dessous les chaînes suivantes pour afficher les lrmax :

Le lrmax de et de est 3 car la chaîne "CGT" est à la fois suffixe de la première chaîne et préfixe de la seconde chaîne.

def lrmax(m1, m2):
    lg_max = min(len(m1), len(m2))
    for i in range(lg_max-1, -1, -1):
        prefixe = m2[:i]
        lg = len(prefixe)
        suffixe = m1[len(m1)-lg:]
        if prefixe == suffixe:
            return lg
    return 0

Première modélisation

La relation « le lrmax est égal à l » (relation entre mots du spectre) peut être représentée par un graphe orienté.

8. On obtient les graphes suivants (réalisés avec le module graphviz de Python) :

Pour le spectre {GTGA,ATGA,GACG,CGTG,ACGT,TGAC} (avec l=4) on obtient ce graphe :

Pour le spectre {TAC,ACC,ACG,CAC,CCG,CGT,CGC,GCA,GTA} (donc avec l=3) on obtient ce graphe :

import graphviz

def model_1(sp, l):
    g = graphviz.Digraph()
    for m1 in sp:
        g.node(m1)
        for m2 in sp:
            if lrmax(m1, m2) == l - 1:
                g.edge(m1, m2)
    return g

9. Une solution au problème de reconstruction peut être obtenue si on trouve un chemin (graphe) (ou un circuit (graphe)) qui passe par tous les sommets du graphe. Ce problème se rapproche de celui de la recherche d’un graphe hamiltonien qui est un problème classique (datant des travaux de William Rowan Hamilton au milieu du XIX^e siècle).

10. Pour le premier des graphes ci-dessus, le seul moyen de passer par le sommet ATGA est de commencer par lui (son degré entrant est nul). Ensuite comme le lrmax de chaque sommet au suivant est 3, il suffit, en suivant les flèches, de rajouter la dernière lettre de chaque sommet parcouru, ce qui donne ATGACGTGAC. L’outil qui ouvre cet article donne le spectre de ATGACGTGAC qui est bien {ATGA,TGAC,GACG,ACGT,CGTG,GTGA}. On note que le chemin est bien hamiltonien (mais ce n’est pas un circuit, il ne revient pas au point de départ).

Pour le second graphe la situation est plus complexe, parce qu’on ne sait pas d’où partir. Comme graphviz a eu la gentillesse de proposer un ordre (de haut en bas) sur les sommets, on peut commencer par le sommet ACC qui est tout en haut puis suivre les sommets de haut en bas (chacun est relié au précédent) pour obtenir ACCGTACGCAC dont le spectre est {ACC,CCG,CGT,GTA,TAC,ACG,CGC,GCA,CAC}.

Deuxième modélisation

Cette fois-ci les sommets du graphe ne sont plus les mots du spectre mais leurs préfixes ou suffixes de longueur l-1 (obtenus en enlevant la première ou dernière lettre du mot).

11. Pour le premier spectre {GTGA,ATGA,GACG,CGTG,ACGT,TGAC} (l==4) on obtient ce graphe :

C’est le même que la première modélisation, au renommage des sommets près.

Pour le second spectre {TAC,ACC,ACG,CAC,CCG,CGT,CGC,GCA,GTA} (donc avec l==3) on obtient ce graphe :

On constate qu’il y a moins de sommets qu’avec le modèle précédent (6 au lieu de 9).

import graphviz

def model_2(sp, l):
    g = graphviz.Digraph()
    for m in sp:
        g.edge(m[:-1], m[1:])
    return g

On voit apparaître là un graphe de de Bruijn. Ces graphes ont été décrits en 1946 par de Bruijn, dans cet article.

12. Le parcours du premier graphe que l’on a vu à la question 10 n’est pas seulement hamiltonien, il est aussi eulérien. Ce chemin eulérien n’est pas un circuit. Le graphe ne possède pas de circuit eulérien parce qu’un tel circuit passerait par le sommet ATG et aucun arc ne permet d’arriver à ce sommet.

Le second graphe par contre possède un circuit eulérien : on fait un tour à gauche (qui revient en haut) puis un tour à droite avec AC→CC→CG→GT→TA→AC→CG→GC→CA→AC qui redonne la solution ACCGTACGCAC au problème de la reconstruction.

Remarque : tout graphe de de Bruijn est à la fois hamiltonien et eulérien, c’est ce que l’on retrouve ici.

13. Un chemin eulérien fournit une solution au problème de reconstruction, tout simplement en le parcourant et en recollant les mots au fur et à mesure du parcours (en évitant de répéter le préfixe du prochain mot puisqu’il est déjà présent comme suffixe du mot précédent). Pour un circuit il faut penser à répéter à la fin le sommet de départ.

14. On doit prouver deux implications :

• Si le graphe (orienté) possède un circuit eulérien alors tous ses sommets ont un degré entrant égal au degré sortant.

• Si chaque sommet a des degrés entrant et sortant égaux entre eux alors il existe un circuit eulérien sur ce graphe.

Histoire de la preuve

On peut, pour montrer cela, « désorienter » le graphe en remplaçant chaque arc par une arête. On obtient alors un multigraphe non orienté et connexe. De plus chaque sommet de ce multigraphe a pour degré la somme des degrés entrant et sortant du graphe orienté initial, c’est-à-dire, par hypothèse, le double du degré sortant. On peut alors appliquer à ce multigraphe le théorème d’Euler (1741) qui donne l’existence d’un circuit eulérien sur ce graphe. Ce circuit est également un circuit eulérien sur le graphe de départ.

Voici l’article de 1741 :

On constate qu’Euler n’a pas réellement prouvé l’existence du circuit eulérien. La première preuve publiée de cette existence est celle donnée par Carl Hierholzer dans cet article de 1873 :

On remarque que la preuve de Hierholzer s’inscrit dans le constructivisme (mathématiques) puisque pour prouver l’existence du circuit eulérien, il fournit un algorithme permettant de construire ce circuit. D’où l’idée d’une preuve rapide :

Une preuve constructiviste est l’algorithme qui sera donné dans la question 20.

Le graphe possédant un circuit eulérien et qui servira d’exemple par la suite, peut se coder ainsi en Python (dictionnaire dont les valeurs sont des listes d’adjacence) :

{'AC':['CC', 'CG'], 'CC':['CG'], 'CG':['GT', 'GC'], 
            'GT':['TA'], 'GC':['CA'], 'TA':['AC'], 'CA':['AC']}

15. Pour savoir si un graphe orienté G possède un circuit eulérien, on peut donc lui appliquer cette fonction :

def presence_circuit_eulerien(g):
    return all(degre_entrant(g, v) == degre_sortant(g, v) for v in g)

def degre_sortant(g, v):
    return len(g[v])

def degre_entrant(g, v):
    return sum(v in g[s] for s in g)

16. Pour construire un circuit partant d’un sommet v du graphe G, on peut utiliser cette fonction qui effectue un parcours en profondeur du graphe :

def construction_circuit(g, v):
    dep = v
    copie_g = {s:g[s].copy() for s in g}
    circuit = [v]
    while copie_g[v] and (len(circuit) <= 1 or circuit[-1] != dep):
        v = copie_g[v].pop()
        circuit.append(v)
    return circuit

Appliquée au graphe précédent avec le sommet CG, cette fonction donne le circuit ['CG', 'GC', 'CA', 'AC', 'CG'] (une liste Python). Ce cycle est de longueur 4 (arcs non pointillés) :

17. Pour supprimer (« remove ») les arcs d’un circuit C de G, on peut utiliser cette fonction :

def enleve_circuit(g, c):
    for i, v in enumerate(c[:-1]):
        g[v].remove(c[i+1])

Appliquée aux graphe et cycle précédents, elle donne ce nouveau graphe qui n’est plus connexe :

L’algorithme de Hierholzer (1873) consiste à enlever des circuits jusqu’à ce qu’il n’en reste plus puis à joindre ces circuits en un circuit eulérien. On voit ici que ce qui reste du graphe après avoir enlevé un premier circuit, est constitué de sommets isolés et d’un circuit.

18. Pour recoller des circuits, il faut d’abord trouver un sommet qui leur soit commun. Ce que fait cette fonction :

def sommet_commun(g, c):
    for v in c:
        if v in g and g[v]:
            return v

Appliquée au graphe ci-dessus et au circuit C enlevé précédemment, elle renvoie le sommet AC qui est bien commun au graphe et au circuit.

19. Pour recoller deux circuits dont un sommet commun est connu, on peut leur appliquer cette fonction :

def fusion_circuits(c1, c2, v):
    assert v in c1.keys() and v in c2.keys()
    c3 = []
    fusion = False
    for s in c1:
        c3.append(s)
        if s == v and not fusion:
            c3.extend(rotation(c2, v))
            fusion = True
    return c3

def rotation(c, v):
    ind = 1
    tmp = []
    while c[ind] != v:
        tmp.append(c[ind])
        ind += 1
    return c[ind+1:]+tmp+[v]

La fusion des deux circuits précédents sur le sommet CG produit ce circuit eulérien : ['AC', 'CG', 'GT', 'TA', 'AC', 'CC', 'CG', 'GC', 'CA', 'AC'].

20. Pour calculer un circuit eulérien selon l’algorithme de Hierholzer, on utilise les 4 fonctions précédentes :

def circuit_eulerien(g, v):
    assert presence_circuit_eulerien(g)
    def graphe_sans_arete(g):
        return all(g[v] == [] for v in g)
    
    copie_g = {s:g[s].copy() for s in g}
    circuit = construction_circuit_b(copie_g, v)
    while not graphe_sans_arete(copie_g):
        v = sommet_commun(copie_g, circuit)
        circuit = fusion_circuits(circuit, construction_circuit_b(copie_g, v), v)
    return circuit

def construction_circuit_b(g, v):    
    dep = v
    circuit = [v]
    while g[v] and (len(circuit) <= 1 or circuit[-1] != dep):
        v = g[v].pop()
        circuit.append(v)
    return circuit

Cet algorithme est de complexité linéaire.

21. Le circuit eulérien trouvé avec la fonction de la question 20 sur le graphe précédent donne CGTACCGCACG.

def obtenir_sequence(circuit):
    s = circuit[0]
    for v in circuit[1:]:
        s += v[-1]
    return s

22. Si, au lieu de partir du sommet CG, on part du sommet AC, on obtient la séquence ACCGTACGCAC qui n’est pas la même que CGTACCGCACG, même à rotation près. Le 4 circuits eulériens issus du sommet AC donnent d’ailleurs 4 séquences possibles à rotations près :

• ACCGTACGCAC
• ACGTACCGCAC
• ACCGCACGTAC
• ACGCACCGTAC

On peut le vérifier en suivant le graphe du regard, à cet effet le revoici :

L’idée de reconstituer une séquence ADN en cherchant un chemin eulérien dans un graphe de de Bruijn remonte à 2001 avec cet article de recherche :