¿Existe alguna función que pueda calcular una puntuación para secuencias alineadas dados los parámetros de alineación?

Question

Intento marcar las secuencias ya alineadas. Vamos a decir

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE' 
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE' 

con parámetros dados

substitution matrix : blosum62 
gap open penalty : -5 
gap extension penalty : -1 

me veía a través del libro de cocina Biopython pero todo lo que puedo conseguir es la sustitución blogsum62 matriz pero siento que debe tener alguien ya implementado este tipo de bibliotecas .

¿Alguien puede sugerir alguna biblioteca o el código más corto que pueda resolver mi problema?

Thx de antemano

Answer 1

Jessada,

la matriz BLOSUM62 (nótese la ortografía;) está en Bio.SubsMat.MatrixInfo y es un diccionario con tuplas resolver a la puntuación (de modo ('A', 'A') vale 4 pts) . No tiene las lagunas, y es solo un triángulo de la matriz (por lo que podría tener ('T', 'A') pero no ('A', 'T'). Hay algunas funciones auxiliares en Biopython, incluyendo algunos en Bio.Pairwise, pero esto es lo que se me ocurrió como respuesta:.

from Bio.SubsMat import MatrixInfo 

def score_match(pair, matrix): 
    if pair not in matrix: 
     return matrix[(tuple(reversed(pair)))] 
    else: 
     return matrix[pair] 

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e): 
    score = 0 
    gap = False 
    for i in range(len(seq1)): 
     pair = (seq1[i], seq2[i]) 
     if not gap: 
      if '-' in pair: 
       gap = True 
       score += gap_s 
      else: 
       score += score_match(pair, matrix) 
     else: 
      if '-' not in pair: 
       gap = False 
       score += score_match(pair, matrix) 
      else: 
       score += gap_e 
    return score 

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE' 
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE' 

blosum = MatrixInfo.blosum62 

score_pairwise(seq1, seq2, blosum, -5, -1) 

que devuelve 82 para su alineación Hay casi certianly maneras más bonitas que hacer todo esto, pero que debe ser un buen empezar.

Answer 2

BLOSUM62 es una Dictonary de 276 artículos.

yo prefería para completar con los elementos que faltan, porque representa una iteración de sólo 276 vueltas, mientras que las secuencias a analizar ar Es probable que tenga más de 276 elementos. En consecuencia, si encuentra el puntaje de cada par con la ayuda de la función score_match(), esta función deberá realizar la prueba if pair not in matrix para cada uno de los elementos de las secuencias, es decir, mucho más de 276 veces.

Otra cosa que lleva mucho tiempo: cada score += something crea un nuevo entero y vincula el nombre puntuación a este nuevo objeto. Cada enlace toma una cantidad de tiempo que no existe con un flujo de enteros por un generador que se agrega inmediatamente a la cantidad actual.

from Bio.SubsMat.MatrixInfo import blosum62 as blosum 
from itertools import izip 

blosum.update(((b,a),val) for (a,b),val in blosum.items()) 

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e, gap = True): 
    for A,B in izip(seq1, seq2): 
     diag = ('-'==A) or ('-'==B) 
     yield (gap_e if gap else gap_s) if diag else matrix[(A,B)] 
     gap = diag 

seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE' 
seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE' 

print sum(score_pairwise(seq1, seq2, blosum, -5, -1)) 

Este score_pairwise() es una función de generador porque hay rendimiento en lugar de retorno.