Traceback nell'implementazione della programmazione dinamica dell'algoritmo di Needleman-Wunsch

Question

Ho quasi fatto funzionare la mia implementazione di needleman-wunsch ma sono confuso su come gestire il traceback in un caso specifico.

L'idea è che per ricostruire la sequenza (il percorso più lungo) ricalcoli per determinare la matrice da cui proviene il punteggio. Il caso limite con cui sto avendo un problema è quando il punteggio in basso a destra non è nella matrice di corrispondenza, ma è nella matrice di inserimento colonna (significa che la risultante sequenza indietro dovrebbe avere un inserto.

Queste sequenze sono in corso registrato nel formato a2m, in cui gli inserimenti nella sequenza sono registrati come caratteri minuscoli Quindi nell'output finale l'allineamento di ZZ a AAAC deve essere AAac. Quando rintraccio a mano finisco con AAAc perché visito solo il Ic matrix una volta. Here è una foto della mia lavagna.Poiché puoi vedere, ho tre frecce nere e una freccia verde, motivo per cui il mio traceback mi dà AAAc. Dovrei contare la prima cella, quindi fermarmi alla posizione 1, 1? Non sono sicuro di come cambierei il modo in cui implego mentito a farlo.

Si noti che la matrice di sostituzione utilizzata qui è BLOSUM62. Le relazioni di ricorrenza sono

M(i,j) = max(Ic(i-1,j-1)+subst, M(i-1,j-1)+subst, Ir(i-1,j-1)+subst) 
Ic(i,j) = max(Ic(i,j-1)-extend, M(i,j-1)-open, Ir(i,j-1)-double) 
Ir(i,j) = max(Ic(i-1,j)-double, M(i-1,j)-open, Ir(i-1,j)-extend) 

MODIFICA: questa è la funzione traceback_col_seq riscritta per essere più pulita. Nota che score_cell ora restituisce thisM, thisC, thisR invece del massimo di questi. Questa versione segna l'allineamento come AaAc, avendo ancora lo stesso problema, e ora con un altro problema del perché torna in Ic di nuovo a 1,2. Questa versione è molto più leggibile comunque.

def traceback_col_seq(self): 
    i, j = self.maxI-1, self.maxJ-1 
    self.traceback = list() 
    matrixDict = {0:'M',1:'Ic',2:'Ir',3:'M',4:'Ic',5:'Ir',6:'M',7:'Ic',8:'Ir'} 
    while i > 0 or j > 0: 
     chars = self.col_seq[j-1] + self.row_seq[i-1] 
     thisM, thisC, thisR = self.score_cell(i, j, chars) 
     cell = thisM + thisC + thisR 
     prevMatrix = matrixDict[cell.index(max(cell))] 
     print(cell, prevMatrix,i,j) 
     if prevMatrix == 'M': 
      i -= 1; j -= 1 
      self.traceback.append(self.col_seq[j]) 
     elif prevMatrix == 'Ic': 
      j -= 1 
      self.traceback.append(self.col_seq[j].lower()) 
     elif prevMatrix == 'Ir': 
      i -= 1 
      self.traceback.append('-') 
    return ''.join(self.traceback[::-1]) 

Qui è la classe pitone che genera la matrice di programmazione dinamica e risale l'allineamento. C'è anche una funzione di punteggio utilizzata per verificare se l'allineamento prodotto è corretto.

class global_aligner(): 
    def __init__(self, subst, open=12, extend=1, double=3): 
     self.extend, self.open, self.double, self.subst = extend, open, double, subst 
    def __call__(self, row_seq, col_seq): 
     #add alphabet error checking? 
     score_align(row_seq, col_seq) 
     return traceback_col_seq() 
    def init_array(self): 
     """initialize three numpy arrays, set values in 1st column and row""" 
     self.M = zeros((self.maxI, self.maxJ), float) 
     self.Ic = zeros((self.maxI, self.maxJ), float) 
     self.Ir = zeros((self.maxI, self.maxJ), float) 
     for i in xrange(1,self.maxI): 
      self.M[i][0], self.Ic[i][0], self.Ir[i][0] = \ 
        -float('inf'), -float('inf'), -(self.open+self.extend*(i-1)) 
     for j in xrange(1,self.maxJ): 
      self.M[0][j], self.Ir[0][j], self.Ic[0][j] = \ 
        -float('inf'), -float('inf'), -(self.open+self.extend*(j-1)) 
     self.Ic[0][0] = self.Ir[0][0] = -float('inf') 
    def score_cell(self, i, j, chars): 
     """score a matrix cell based on the 9 total neighbors (3 each direction)""" 
     thisM = [self.M[i-1][j-1]+self.subst[chars], self.Ic[i-1][j-1]+ \ 
       self.subst[chars], self.Ir[i-1][j-1]+self.subst[chars]] 
     thisC = [self.M[i][j-1]-self.open, self.Ic[i][j-1]-self.extend, \ 
         self.Ir[i][j-1]-self.double] 
     thisR = [self.M[i-1][j]-self.open, self.Ic[i-1][j]-self.double, \ 
       self.Ir[i-1][j]-self.extend] 
     return max(thisM), max(thisC), max(thisR) 
    def score_align(self, row_seq, col_seq): 
     """build dynamic programming matrices to align two sequences""" 
     self.row_seq, self.col_seq = list(row_seq), list(col_seq) 
     self.maxI, self.maxJ = len(self.row_seq)+1, len(self.col_seq)+1 
     self.init_array() #initialize arrays 
     for i in xrange(1, self.maxI): 
      row_char = self.row_seq[i-1] 
      for j in xrange(1, self.maxJ): 
       chars = row_char+self.col_seq[j-1] 
       self.M[i][j], self.Ic[i][j], self.Ir[i][j] = self.score_cell(i, j, chars) 
    def traceback_col_seq(self): 
     """trace back column sequence in matrices in a2m format""" 
     i, j = self.maxI-1, self.maxJ-1 
     self.traceback = list() 
     #find which matrix to start in 
     #THIS IS WHERE THE PROBLEM LIES I THINK 
     cell = (self.M[i][j], self.Ic[i][j], self.Ir[i][j]) 
     prevMatrix = cell.index(max(cell)) 
     while i > 1 and j > 1: 
      if prevMatrix == 0: #M matrix 
       i -= 1; j -= 1 #step up diagonally 
       prevChars = self.row_seq[i-1]+self.col_seq[j-1] 
       diag = self.score_cell(i, j, prevChars) #re-score diagonal cell 
       prevMatrix = diag.index(max(diag)) #determine which matrix that was 
       self.traceback.append(self.col_seq[j]) 
      elif prevMatrix == 1: #Ic matrix 
       j -= 1 
       prevChars = self.row_seq[i-1]+self.col_seq[j-1] 
       left = self.score_cell(i, j, prevChars) 
       prevMatrix = left.index(max(left)) 
       self.traceback.append(self.col_seq[j].lower()) 
      elif prevMatrix == 2: #Ir matrix 
       i -= 1 
       prevChars = self.row_seq[i-1]+self.col_seq[j-1] 
       up = self.score_cell(i, j, prevChars) 
       prevMatrix = up.index(max(up)) 
       self.traceback.append('-') 
     for j in xrange(j,0,-1): #hit top of matrix before ending, add chars 
      self.traceback.append(self.col_seq[j-1]) 
     for i in xrange(i,0,-1): #hit left of matrix before ending, add gaps 
      self.traceback.append('-') 
     print(''.join(self.row[::-1])) 
     return ''.join(self.traceback[::-1]) 
    def score_a2m(self, s1, s2): 
     """scores an a2m alignment of two sequences. I believe this function correctly 
     scores alignments and is used to test my alignments. The value produced by this 
     function should be the same as the largest value in the bottom right of the three 
     matrices""" 
     s1, s2 = list(s1.strip('.')), list(s2.strip('.')) 
     s1_pos, s2_pos = len(s1)-1, len(s2)-1 
     score, gap = 0, False 
     while s1_pos >= 0 and s2_pos >= 0: 
      if s1[s1_pos].islower() and gap is False: 
       score -= self.open; s1_pos -= 1; gap = True 
      elif s1[s1_pos].islower() and gap is True: 
       score -= self.extend; s1_pos -= 1 
      elif s2[s2_pos].islower() and gap is False: 
       score -= self.open; s2_pos -= 1; gap = True 
      elif s2[s2_pos].islower() and gap is True: 
       score -= self.extend; s2_pos -= 1 
      elif s1[s1_pos] == '-' and gap is False: 
       score -= self.open; s1_pos -= 1; s2_pos -= 1; gap = True 
      elif s1[s1_pos] == '-' and gap is True: 
       score -= self.extend; s1_pos -= 1; s2_pos -= 1 
      elif s2[s2_pos] == '-' and gap is False: 
       score -= self.open; s1_pos -= 1; s2_pos -= 1; gap = True 
      elif s2[s2_pos] == '-' and gap is True: 
       score -= self.extend; s1_pos -= 1; s2_pos -= 1 
      elif gap is True: 
       score += self.subst[s1[s1_pos].upper() + s2[s2_pos].upper()] 
       s1_pos -= 1; s2_pos -= 1; gap = False 
      else: 
       score += self.subst[s1[s1_pos].upper() + s2[s2_pos].upper()] 
       s1_pos -= 1; s2_pos -= 1 
     if s1_pos >= 0 and gap is True: 
      score -= self.extend*s1_pos 
     elif s1_pos >= 0 and gap is False: 
      score -= self.open+s1_pos*self.extend 
     if s2_pos >= 0 and gap is True: 
      score -= self.extend*s2_pos 
     elif s2_pos >= 0 and gap is False: 
      score -= self.open+s2_pos*self.extend 
     return score 


test = global_aligner(blosumMatrix) 
s1,s2 = 'ZZ','AAAC' 
test.score_align(s1, s2) 
align = test.traceback_col_seq() 
print('This score: ', test.score_a2m(s1,align) 
print('Correct score: ', test.score_a2m(s1,'AAac')) 

BLOSUM parser

def parse_blosum(blosumFile): 
    blosumMatrix, commentFlag = dict(), False 
    for line in blosumFile: 
     if not line.startswith('#') and not commentFlag: 
      alphabet = line.rstrip().split() 
      commentFlag = True 
     elif commentFlag: 
      line = line.rstrip().split() 
      thisChar, line = line[0], line[1:] 
      for x in xrange(len(line)): 
       alphaChar, thisValue = alphabet[x], line[x] 
       blosumMatrix[thisChar+alphaChar] = int(thisValue) 
    return blosumMatrix 

Answer 1

def traceback_col_seq(self): 
    """ 
    Traces back the column sequence to determine final global alignment. 
    Recalculates the score using score_cell. 
    """ 
    i, j = self.maxI-1, self.maxJ-1 
    self.traceback = list() #stores final sequence 
    matrixDict = {0:'M',1:'Ic',2:'Ir'} #mapping between numeric value and matrix 
    chars = self.col_seq[j-1] + self.row_seq[i-1] #store first characters 
    thisM, thisC, thisR = self.score_cell(i,j, chars) 
    cell = max(thisM), max(thisC), max(thisR) #determine where to start 
    prevMatrix = matrixDict[cell.index(max(cell))] #store where to go first 
    while i > 0 or j > 0: 
     #loop until the top left corner of the matrix is reached 
     if prevMatrix == 'M': 
      self.traceback.append(self.col_seq[j-1]) 
      i -= 1; j -= 1 
      prevMatrix = matrixDict[thisM.index(max(thisM))] 
     elif prevMatrix == 'Ic': 
      self.traceback.append(self.col_seq[j-1].lower()) 
      j -= 1 
      prevMatrix = matrixDict[thisC.index(max(thisC))] 
     elif prevMatrix == 'Ir': 
      self.traceback.append('-') 
      i -= 1 
      prevMatrix = matrixDict[thisR.index(max(thisR))] 
     chars = self.col_seq[j-1] + self.row_seq[i-1] #store next characters 
     thisM, thisC, thisR = self.score_cell(i,j, chars) #rescore next cell 
    return ''.join(self.traceback[::-1]) 

Answer 2

Come accennato con i suoi commenti su freccia di colore, il problema fondamentale è che la seconda freccia orizzontale è sempre etichettato come un movimento M (freccia nera) quando è davvero un movimento in Ic (freccia verde). Questo sembra accadere perché la variabile prevMatrix indica la migliore matrice in (i-1, j-1), (i-1, j) o (i, j-1). Questa è la matrice nella cella precedente dalla prospettiva in avanti. Poiché stai facendo trace * back * a questo punto, la cella che stai "provenendo da" nella traccia è in realtà (i, j) - quella in cui ti trovi attualmente. Ciò determina la direzione in cui ti stai muovendo e quindi se si sta allineando a un divario o no. L'opzione migliore è probabilmente quella di avere una variabile che tenga traccia della matrice corrente dall'iterazione del ciclo all'iterazione del ciclo e utilizzarla per determinare quale carattere aggiungere alla stringa di allineamento. Ad ogni iterazione puoi aggiornarlo dopo aver determinato quale cella stai per passare.

Penso che apportare queste modifiche chiarirà il problema che si sta verificando con la funzione traceback_col_sequence riscritta, ma in via preliminare sembra che non si tenga conto delle informazioni che si sono acquisite tracciando indietro. Stai riutilizzando score_cell(), che sembra essere stato progettato per calcolare i punteggi in futuro. Quando stai risalendo, conosci già il finale, quindi non vuoi scegliere il punteggio massimo per la cella precedente. Vuoi scegliere il punteggio massimo che potrebbe portare alla posizione nella matrice in cui ti trovi attualmente.Ad esempio, se sai che la matrice in cui ti trovi nel tuo traceback è M, allora sai che non ci sei arrivato estendendo un gap, quindi non dovresti considerare queste opzioni.