Choix du coup minimal

Bonjour.
Lors de la phase 2, on remonte les nœuds en prenant le coût minimal parmi les 3 coûts amenant au nœud.
Mais comment choisir si il y a 2 minimums parmi les 3 coûts?
Dans le programme python fourni par la suite on choisit apparemment le coût de la substitution si ce dernier est un coût minimum. Est-ce exact? Est-ce toujours le cas? d’où des erreurs ensuite lors du dernier exemple ou on a inséré deux séquences. Est-ce toujours exact?

Il peut exister plusieurs alignements optimaux pour une même paire de séquences.

L’algorithme détaillé dans le cours “se contente” d’en trouver un. De ce fait, quand la situation que vous décrivez se présente, l’algorithme retient un noeud parmi les possibles et poursuit son parcours.

Ce bout de code non récursif devrait corriger le problème et fournir tous les chemins:

    # un exemple de phase pour tous les parcours minimaux avec une pile

def phase2bis(dna1, dna2, costs):
    pile=[([],[],len(dna1),len(dna2))]
    while len(pile)>0:
        a=pile.pop()
        r1,r2,i,j=a
        if i==0 and j==0:
            # à ce stade il nous reste à retourner les listes, les transformer en chaines et afficher
            s1 = "".join(reversed(r1))
            s2 = "".join(reversed(r2))
            print("-------------------------------")
            print(s1)
            print(s2)
        else:
            # la valeur courante
            c = costs[i,j]
            # si on est au bord, les formules en i-1 ou j-1 
            # ne vont pas faire ce qu'on veut, il faut traiter 
            # ces cas à part
            if i == 0:                  # bord = insertion
                rr1=r1+["-"]
                rr2 =r2 +[dna2[j-1]]
                pile.append((rr1,rr2,i,j-1))
            elif j == 0:                # bord = insertion
                rr1 = r1+[dna1[i]]
                rr2 = r2+ ["-"]
                pile.append((rr1,rr2,i-1,j))
            # dans le milieu du tableau on regarde de quelle direction nous vient le minimum
            else :
                if c == costs[i-1,j-1] + substitution_cost(dna1[i-1], dna2[j-1]):  # substitution
                    # s'agit-t-il d'une vraie substitution ? 
                    same = dna1[i-1] == dna2[j-1]
                    rr1= r1+ [outline(dna1[i-1], same)]
                    rr2= r2 + [outline(dna2[j-1], same)]
                    pile.append((rr1,rr2,i-1,j-1))
                if c == costs[i,j-1] + insertion_cost(dna2[j-1]):    # insertion
                    rr1 = r1 + ['-']
                    rr2 = r2 + [dna2[j-1]]
                    pile.append((rr1,rr2,i,j-1)) 
                if c == costs[i-1,j] + insertion_cost(dna1[i-1]):    # insertion
                    rr1=r1+[dna1[i-1]]
                    rr2 = r2+ ['-']
                    pile.append((rr1,rr2,i-1,j))
    return None

Il faut corriger bien sûr phase1 car je suis passé avec un dictionnaire:

def phase1(dna1, dna2):
    """
    Première phase de Needleman et Wunsch itératif
    Élabore itérativement le dictionnaire des coûts 
      par un parcours en diagonale
    On obtient un dictionnaire de taille 
      [len(dna1) + 1] x [len(dna2) + 1]
    Renvoie le dictionnaire en valeur
    """
    # initialisations
    len1 = len(dna1)
    len2 = len(dna2)
    # le tableau est initialisé à zéro
    costs = init_costs(len1, len2)

    # le parcours en diagonale - cf ci-dessus
    for c in range(len1 + len2 + 1):
        for i in range(c + 1):
            # on déduit j de c et i
            j = c - i
            # on ne considère que ceux qui tombent dans le rectangle 
            if 0 <= i <= len1 and 0 <= j <= len2:
                if i == 0 and j == 0:  # le coin en haut a gauche
                    costs[i,j] = 0
                elif j == 0:           # sur un bord : insertion 
                    costs[i,j] = costs[i-1,j] + insertion_cost(dna1[i-1])
                elif i == 0:           # l'autre bord : insertion
                    costs[i,j] = costs[i,j-1] + insertion_cost(dna2[j-1])
                else:                  # au milieu
                    costs[i,j] = min(
                        # substitution
                        costs[i-1,j-1] + substitution_cost(dna1[i-1], dna2[j-1]),
                        # insertion
                        costs[i,j-1] + insertion_cost(dna2[j-1]),
                       # insertion
                        costs[i-1,j] + insertion_cost(dna1[i-1]))
    # on renvoie le résultat
    return costs

Avec bien sûr celui-ci qui résume

# on peut maintenant écrire une fonction de commodité
def needleman_wunschbis(dna1, dna2):
    # on calcule les coûts avec la phase1
    costs = phase1(dna1, dna2)
    # on calcule la distance
    d = costs[len(dna1),len(dna2)]
    # on passe à la phase 2
    print("distance = ", d)
    phase2bis(dna1, dna2, costs)

Ainsi si nous mettons un coût de 2 pour une substitution et de 1 pour une insertion nous obtenons ceci:

# et sur un exemple cela donne
sample1 = "AC"
sample2 = "AT"

needleman_wunschbis(sample1, sample2)

Ce qui renvoie :

distance =  2
-------------------------------
A-C
AT-
-------------------------------
AC-
A-T
-------------------------------
A*C*
A*T*

Ou encore toujours avec le coût de 2 sur la substitution:

sample1 = "ACCTCTGTATCTATTCGGCATCGATCAT"
sample2 = "ACCTCGTGTATCTCTTCGGCATCATCAT"

needleman_wunschbis(sample1, sample2)

Qui nous donne :

distance =  4
-------------------------------
ACCTC-TGTATCT-ATTCGGCATCGATCAT
ACCTCGTGTATCTC-TTCGGCATC-ATCAT
-------------------------------
ACCTC-TGTATCTA-TTCGGCATCGATCAT
ACCTCGTGTATCT-CTTCGGCATC-ATCAT
-------------------------------
ACCTC-TGTATCT*A*TTCGGCATCGATCAT
ACCTCGTGTATCT*C*TTCGGCATC-ATCAT

Avec un coût de 1 sur la substitution nous n’obtenons qu’un seul chemin dans les 2 exemples précédents.
Cela semble fonctionner, non? Et cela a l’air assez rapide.
Merci encore pour ce Mooc très interessant!!

Merci @DraganP j’avais relevé la même interrogation.