Capes NSI 2020 Epreuve 1, Problème 1

mardi 28 juillet 2020
par  Alain BUSSER , Sébastien HOARAU

Dans le séquençage de l’ADN, on cherche une chaîne de caractères écrite dans l’alphabet $\{A,C,G,T\}$ à partir de ses mots de longueur l (un entier typiquement beaucoup plus petit que la longueur de la chaîne à reconstituer).

Le problème n°2 fera l’objet d’un article ultérieur.

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

  1. S[0:l]
  2. S[1:l+1]
  3. S[2:l+2]
  4. .
  5. .
  6. .
  7. S[N-1:l+N-1]

Télécharger

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) :

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

Télécharger

Version détaillée

Dans la question suivante, on se basera sur cette variante plus détaillée :

  1. def mots_longueur_l(sn, l):
  2.     l_mots = []
  3.     for indice in range(len(sn)-l+1):
  4.         mot = ''
  5.         for delta in range(l):
  6.             mot += sn[indice + delta]
  7.         l_mots.append(mot)
  8.     return l_mots

Télécharger

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 4l 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,4l).

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 :

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

Télécharger

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

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

Télécharger

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.

Version détaillée

On peut modifier la version détaillée de la question 2 pour avoir le même résultat sans passer par les ensembles (donc sans utiliser la question précédente) :

  1. def liste_mots(sn, l):
  2.     l_mots = []
  3.     for indice in range(len(sn)-l+1):
  4.         mot = ''
  5.         for delta in range(l):
  6.             mot += sn[indice + delta]
  7.         if mot not in l_mots:
  8.             l_mots.append(mot)
  9.     return l_mots

Télécharger

C’est ce qui a été fait pour la version interactive qui ouvre cet article (mais en JavaScript, pas en Python).

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 :

ATGC GGTA
TGGCGT CGTAAATG
GCTAGGCTAA AGGCTAAGTCGAT
TCTAGCCAGCTAGC TAGCCAGCTAGCACT

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.

Calcul du lrmax

Ce script Python permet de calculer le lrmax, sous forme d’une fonction :

  1. def lrmax(m1, m2):
  2.     lg_max = min(len(m1), len(m2))
  3.     for i in range(lg_max-1, -1, -1):
  4.         prefixe = m2[:i]
  5.         lg = len(prefixe)
  6.         suffixe = m1[len(m1)-lg:]
  7.         if prefixe == suffixe:
  8.             return lg
  9.     return 0

Télécharger

La boucle doit se faire par valeurs décroissantes de la longueur, car sinon elle renverrait une valeur non maximale.

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 :

%3 GACG GACG ACGT ACGT GACG->ACGT CGTG CGTG ACGT->CGTG TGAC TGAC TGAC->GACG GTGA GTGA CGTG->GTGA ATGA ATGA ATGA->TGAC GTGA->TGAC

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

%3 ACC ACC CCG CCG ACC->CCG CGT CGT CCG->CGT CGC CGC CCG->CGC GTA GTA CGT->GTA TAC TAC GTA->TAC TAC->ACC ACG ACG TAC->ACG GCA GCA CAC CAC GCA->CAC CAC->ACC CAC->ACG CGC->GCA ACG->CGT ACG->CGC

Graphes en Python

La fonction suivante construit, à partir du spectre, le graphe orienté (« digraph ») au format graphviz :

  1. import graphviz
  2.  
  3. def model_1(sp, l):
  4.     g = graphviz.Digraph()
  5.     for m1 in sp:
  6.         g.node(m1)
  7.         for m2 in sp:
  8.             if lrmax(m1, m2) == l - 1:
  9.                 g.edge(m1, m2)
  10.     return g

Télécharger

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 XIXe 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 :

%3 GAC GAC ACG ACG GAC->ACG CGT CGT ACG->CGT TGA TGA TGA->GAC GTG GTG CGT->GTG ATG ATG ATG->TGA GTG->TGA

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 :

%3 AC AC CC CC AC->CC CG CG AC->CG CC->CG GT GT CG->GT GC GC CG->GC TA TA GT->TA TA->AC CA CA GC->CA CA->AC

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

En Python

Ces graphes ont été construits avec cette fonction Python :

  1. import graphviz
  2.  
  3. def model_2(sp, l):
  4.     g = graphviz.Digraph()
  5.     for m in sp:
  6.         g.edge(m[:-1], m[1:])
  7.     return g

Télécharger

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.

Une idée de preuve

Pour montrer cela on peut imaginer que les arcs du graphe sont des ponts fragiles qui sont détruits à chaque passage par l’arc (de toute façon le circuit étant eulérien, on ne passe qu’une seule fois sur chaque pont, comme à Königsberg). Alors chaque passage par un sommet du graphe consomme (détruit)

    • un arc entrant en venant vers le sommet,
    • un arc sortant en repartant du sommet.

Ceci a pour effet de diminuer d’exactement une unité les degrés entrant et sortant à chaque passage par le sommet (l’égalité entre les degrés entrant et sortant est donc un invariant). Le circuit parcouru étant eulérien, à la fin de son parcours il ne reste plus aucun arc. Ce qui prouve qu’initialement les degrés entrant et sortant étaient déjà égaux (au nombre de passages par le sommet).

  • 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) :

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

Télécharger

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

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

Télécharger

degrés

la fonction précédente est basée sur ces calculs de degrés :

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

Télécharger

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

Télécharger

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 :

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

Télécharger

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) :

image/svg+xml %3 AC AC CC CC AC->CC CG CG AC->CG CC->CG GT GT CG->GT GC GC CG->GC TA TA GT->TA TA->AC CA CA GC->CA CA->AC

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

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

Télécharger

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

image/svg+xml %3 AC AC CC CC AC->CC CG CG CC->CG GT GT CG->GT GC GC TA TA GT->TA TA->AC CA CA

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 :

  1. def sommet_commun(g, c):
  2.     for v in c:
  3.         if v in g and g[v]:
  4.             return v

Télécharger

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 :

  1. def fusion_circuits(c1, c2, v):
  2.     assert v in c1.keys() and v in c2.keys()
  3.     c3 = []
  4.     fusion = False
  5.     for s in c1:
  6.         c3.append(s)
  7.         if s == v and not fusion:
  8.             c3.extend(rotation(c2, v))
  9.             fusion = True
  10.     return c3

Télécharger

Où recoller ?

Pour recoller sur le sommet commun, il est nécessaire d’effectuer une rotation adéquate, des sommets du second circuit (pour qu’il commence par le sommet commun). Pour cela, on utilise cette fonction qui implémente une file (structure de données) :

  1. def rotation(c, v):
  2.     ind = 1
  3.     tmp = []
  4.     while c[ind] != v:
  5.         tmp.append(c[ind])
  6.         ind += 1
  7.     return c[ind+1:]+tmp+[v]

Télécharger

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 :

  1. def circuit_eulerien(g, v):
  2.     assert presence_circuit_eulerien(g)
  3.     def graphe_sans_arete(g):
  4.         return all(g[v] == [] for v in g)
  5.    
  6.     copie_g = {s:g[s].copy() for s in g}
  7.     circuit = construction_circuit_b(copie_g, v)
  8.     while not graphe_sans_arete(copie_g):
  9.         v = sommet_commun(copie_g, circuit)
  10.         circuit = fusion_circuits(circuit, construction_circuit_b(copie_g, v), v)
  11.     return circuit

Télécharger

fonction auxiliaire

La fonction construction_circuit a été remplacée par celle-ci :

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

Télécharger

qui évite d’avoir à enlever des circuits.

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.

Du circuit eulérien à la séquence

Cette fonction passe du circuit eulérien de la question 20 à la chaîne de caractères :

  1. def obtenir_sequence(circuit):
  2.     s = circuit[0]
  3.     for v in circuit[1:]:
  4.         s += v[-1]
  5.     return s

Télécharger

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 :

%3 AC AC CC CC AC->CC CG CG AC->CG CC->CG GT GT CG->GT GC GC CG->GC TA TA GT->TA TA->AC CA CA GC->CA CA->AC

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 :


Commentaires