El alineamiento de secuencias es un método computacional que permite inferir información biológica sólo a partir de la información que tienen las secuencias.
Introducción
La forma más fiable de determinar la estructura o función de una molécula biológica es mediante la experimentación directa. Sin embargo, es mucho más fácil obtener la secuencia de ADN del gen correspondiente a un ARN o una proteína, que determinar experimentalmente su función o estructura.
Por tanto parte del reto es simplemente organizar, clasificar y analizar la inmensa riqueza de datos de secuencias, y desarrollar métodos computacionales que puedan inferir información biológica sólo a partir de la secuencia.
En este enlace tiene el proyecto de soporte: https://gitlab.com/xtec/bio/sequence-alignment
Secuencias
Para entender la naturaleza no debes pensar en ella como un inventor que diseña los organismos, sino en una “entidad” que se dedica a arreglar y experimentar con organismos y sus componentes.
Las nuevas secuencias se adaptan a partir de secuencias preexistentes en lugar de inventarse de cero, lo que nos permite analizar las secuencias de forma computacional para reconocer una similitud significativa entre una secuencia nueva y una secuencia sobre la que ya se sabe algo.
Cuando hacemos esto podemos transferir información sobre la estructura y/o la función a la nueva secuencia. Decimos que las dos secuencias relacionadas son homólogas y que estamos transfiriendo información por homología.
Dos secuencias son homologas si se derivan de una secuencia común, es decir, tienen un ancestro común.
La homología no es cuantificable, las secuencias son o no homologas.
Por tanto, una de las preguntas más básicas sobre un gen o proteína es si está relacionado con algún otro gen o proteína:
- La relación de dos proteínas a nivel de secuencia sugiere que son homólogas.
- La relación también sugiere que pueden tener funciones comunes.
Mediante el análisis de muchas secuencias de ADN y proteínas, es posible identificar dominios o motivos que se comparten entre un grupo de moléculas.
Estos análisis de la relación de proteínas y genes se logran alineamiento secuencias.
La mayoría de los problemas en el análisis de secuencias computacionales son esencialmente estadísticos. Las fuerzas evolutivas estocásticas actúan sobre los genomas. Discernir similitudes significativas entre secuencias antiguamente divergentes en medio de un caos de mutación aleatoria, selección natural y deriva genética presenta serios problemas de señal al ruido. Muchos de los métodos de análisis más potentes disponibles utilizan la teoría de la probabilidad.
Alineamiento de proteínas
A la hora de comparar proteínas, podemos comparar la secuencia de aminoácidos que la componen, o comparar los genes que codifican la secuencia de aminoácidos.
Las secuencias de aminoácidos son más informativas porque:
-
Muchos cambios en una secuencia de ADN (especialmente en la tercera posición de un codón) no cambian el aminoácido que se especifica.
-
Muchos aminoácidos comparten propiedades biofísicas relacionadas (por ejemplo, la lisina y la arginina son ambos aminoácidos básicos).
El resultado es que las comparaciones de secuencias de proteínas pueden identificar secuencias homólogas mientras que las comparaciones de secuencias de ADN correspondientes no pueden.
Por tanto, cuando se analiza una secuencia de codificación de nucleótidos a menudo es preferible estudiar su proteína traducida.
En Blast veremos que podemos movernos fácilmente entre los mundos del ADN y las proteínas.
Por ejemplo, la herramienta TBLASTN del sitio web NCBI BLAST permite buscar proteínas relacionadas derivadas de una base de datos de ADN con una secuencia de proteínas.
Esta opción de consulta se logra traduciendo cada secuencia de ADN a las seis proteínas que potencialmente codifica.
Sin embargo todo lo que hemos explicado hasta ahora, en muchos casos es adecuado comparar secuencias de ADN.
Esta comparación puede ser importante para confirmar la identidad de una secuencia de ADN en una búsqueda de bases de datos, para buscar polimorfismos, para analizar la identidad de un fragmento de ADNc clonado o para comparar regiones reguladoras, por poner algunos ejemplos.
Homología
Dos secuencias son homólogas si comparten un ascendente evolutivo común.
No existen grados de homología; las secuencias son homólogas o no.
Las proteínas homólogas casi siempre comparten una estructura tridimensional significativamente relacionada. Un ejemplo de homología es el de la mioglobina humana (NP_005359.1) y la beta globina (NP_000509.1), dos proteínas que están relacionadas lejanamente, pero significativamente.
Se cree que la mioglobina y las cadenas de hemoglobina (alfa, beta y otros) divergieron hace unos 450 millones de años, cerca del momento en que los linajes de peces humanos y cartílagos divergieron.
La mioglobina y la beta globina tienen estructuras muy parecidas, tal y como determina la cristalografía de rayos X tal y como puedes ver en esta imagen:
| Human myoglobin (3RGK) | Human beta globin (subunit of 2H35) |
|---|---|
Estas proteínas son homólogas (descienden de un antepasado común) y comparten estructuras tridimensionales muy similares.
Sin embargo, el alineamiento por parejas de las secuencias de aminoácidos de estas proteínas revela que las proteínas comparten una identidad de aminoácidos muy limitada.
Identidad
Cuando dos secuencias son homólogas, sus secuencias de aminoácidos o de nucleótidos suelen compartir una identidad significativa.
Si bien la homología es una inferencia cualitativa (las secuencias son o no homólogas), la identidad y la semejanza son cantidades que describen la relación de las secuencias.
En particular, dos moléculas pueden ser homólogas sin compartir una identidad estadísticamente significativa de aminoácidos (o nucleótidos).
En la familia de la globina, todos los miembros son homólogos, pero algunos tienen secuencias que tanto han divergido que no comparten ninguna identidad de secuencia reconocible
Por ejemplo:
- La beta globina humana y la neuroglobina humana sólo comparten un 22% de identidad de aminoácidos.
- Las cadenas de globina individuales comparten la misma forma general que la mioglobina, aunque las proteínas de mioglobina y alfa globina sólo comparten un 26% de identidad de aminoácidos.
Por lo general, las estructuras tridimensionales divergen mucho más lentamente que la identidad de la secuencia de aminoácidos entre dos proteínas.
¡Y es la estructura la que define la función de la proteína!
Reconocer este tipo de homología es un problema bioinformático especialmente complejo.
Ortología
Las proteínas que son homólogas pueden ser ortologues o paralogas.
Los ortólogos son secuencias homólogas en distintas especies que tienen un gen antepasado común y que han divergido por especiación.
La historia del gen refleje la historia de la especie, por lo que los genes se llaman ortólogos (orto = exacto). Se supone que los ortólogos tienen funciones biológicas similares, pudiendo deducir la función de una secuencia a partir de otra secuencia homóloga.
Por ejemplo, los humanos y roedores divergieron hace unos 90 millones de años (MYA), momento en que un solo gen de la mioglobina ancestral divergió por especiación, y en ambos casos la función del gen es transportar oxígeno a las células musculares
A continuación tienes un dibujo con el árbol de los ortólogos de la mioglobina:
Paralogía
Los parálogos son secuencias homólogas surgidas por un mecanismo como la duplicación de genes.
En este caso la homología es el resultado de la duplicación génica por lo que ambas copias han bajado una al lado de la otra durante la historia de un organismo.
Los genes se llaman paralógicos (para = en paralelo) porque evolucionan dentro de la misma especie.
Por ejemplo, la globina alfa 1 humana (NP_000549.1) es paraloga en la globina alfa 2 (NP_000508.1). De hecho, estas dos proteínas comparten una identidad de aminoácidos del 100%.
La globina alfa 1 y la globina beta humanas también son paralogos, como todas las proteínas humanas que se muestran a continuación y que son miembros de la familia globina:
Todas las globinas tienen propiedades distintas, incluida la distribución regional en el cuerpo, el tiempo de desarrollo de la expresión génica y la abundancia. Pero aunque todas tienen funciones distintas, estas funciones están relacionadas como proteínas portadoras de oxígeno.
Para saber un poco más sobre estos árboles y la evolución puede consultar este video:
Similitud
A primera vista, decidir que dos secuencias biológicas son similares no es diferente decidir que dos cadenas de texto son similares.
Los algoritmos de semejanza de texto se utilizan para buscar la palabra más parecida en un diccionario cuando una palabra no se encuentra en un diccionario porque está mal escrita. Es lo que conoces como corrector ortográfico básico.
Por tanto, no debe extrañar que los algoritmos de análisis de secuencias compartan mucho parecido con los algoritmos de semejanza de texto.
La distancia de Levenshtein (o distancia de edición) es un algoritmo de 1965 diseñado para comparar dos palabras, que se aplicó posteriormente a secuencias biológicas.
Este algoritmo se basa en que, para transformar una palabra en otra, sólo son necesarias tres operaciones básicas: insertar, suprimir o sustituir un carácter.
Cada vez que se aplica una de estas operaciones, la distancia se incrementa 1 punto, y cuantos más puntos mayor es la distancia de edición de las dos palabras y menos similares son.
Matriz de distancias
La distancia de edición es un algoritmo lineal que utiliza una matriz (un array de dos dimensiones) para guardar los resultados parciales que permiten ir construyendo el resultado final.
Dadas dos secuencias de ácidos nucleicos xs y ys, la matriz de distancias contiene las distancias entre todos los prefijos de la secuencia xs y todos los prefijos de la secuencia ys.
Para empezar, debes llenar la matriz con el resultado del subproblema básico a partir del cual iremos construyendo resultados parciales hasta llegar al resultado final.
Por ejemplo, si quieres comparar las secuencias CGA y AGAT, debes crear una matriz de 4 x 5 con estos valores iniciales (el símbolo - indica un string vacío "")
- A G A T
- 0 1 2 3 4
C 1 ?
G 2
A 3 Como puedes ver en este ejemplo, es evidente que la distancia de edición entre:
- "" i "A" es 1, porque para transformar "" en "A" debemos insertar una A.
- "" i "EN" es 2, porque para transformar "" en "EN" debemos insertar una A y después una G.
A continuación tienes el algoritmo inicial implementado en Python:
import numpy as np
from util import print_matrix
xs = "CGA"
ys = "AGAT"
M = np.zeros((len(xs) + 1, len(ys) + 1), dtype=np.dtype("int8"))
M[:, 0] = np.linspace(0, len(xs), len(xs)+1)
M[0, :] = np.linspace(0, len(ys), len(ys)+1)
print_matrix(xs, ys, M)Si ejecutamos el script podemos ver la tabla resultante:
- A G A T
- 0 1 2 3 4
C 1
G 2
A 3 A continuación debemos ir llenando la matriz a partir de dónde está el interrogante ?, utilizando los valores que están en las celdas M[i-1][j-1], M[i-1][j] i M[i][j-1], resolviendo subproblemas hasta que lleguemos a encontrar la solución al problema.
El algoritmo a utilizar es éste:
for i in range(1, len(xs) + 1):
for j in range(1, len(ys) + 1):
M[i][j] = min(M[i-1][j-1], M[i-1][j], M[i][j-1])
if (xs[i-1] != ys[j-1]):
M[i][j] += 1Tal y como puedes ver en el código, la distancia de edición de una celda será:
- El valor de la distancia mínima para llegar hasta esa celda a partir de la celda
M[i-1][j-1],M[i-1][j]oM[i][j-1] - Si se debe hacer una edición porque la letra de la fila
iy de la columnajno son iguales, incrementamos el valor de la celda en 1.
Por tanto, el segundo paso del algoritmo es
- A G A T
- 0 1 2 3 4
C 1 1 ?
G 2
A 3 El resultado de llenar la matriz será el que se muestra a continuación, donde tenemos el resultado de la distancia entre todos los prefijos posibles que se han ido calculando para encontrar el resultado final.
- A G A T
- 0 1 2 3 4
C 1 1 2 3 4
G 2 2 1 2 3
A 3 2 2 1 2A continuación puedes ver el camino mínino de edición marcado con *, ya que todos los demás caminos, aunque sean correctos, conllevan muchas más operaciones para convertir CGA en AGAT.
- A G A T
- 0 1 2 3 4
C 1 *1 2 3 4
G 2 2 *1 2 3
A 3 2 2 *1 *2Para convertir CGA en AGT tenemos que cambiar la C por una A y añadir una T al final.
Son 2 operaciones (tal y como indica el valor de la celda final), y por tanto la distancia de edición es 2.
Como ya hemos comentado antes, en la tabla están calculados todos los resultados de los prefijos.
Por ejemplo, la distancia entre CG y AGA es 2:
- A G A T
- 0 1 2 3 4
C 1 *1 2 3 4
G 2 2 *1 *2 3
A 3 2 2 1 2Ejecutando el script en Python puedes ver que el resultado del algoritmo es la misma tabla que hemos calculado a mano:
Actividad: Costo computacional
El algoritmo de la distancia de edición se diseñó para comparar palabras y las palabras tienen pocas letras.
Pero si comparamos secuencias de ADN, que pueden ser muy largas, el tiempo que tarda en computar el resultado es cuadrático respecto a la entrada: O(n2).
Esto significa que si tarda 0,00001 segundos en calcular la distancia de una entrada de 5x5, tardará 0,4 segundos en calcular el resultado de una entrada de 1000x1000.
Modifica el archivo para que obtenga dos grandes secuencias de Entrez y calcula la distancia de edición.
Haz un profiling del script tal y como hiciste en Profiling
Implementar este algoritmo en Python es una actividad educativa interesante aunque utilizamos Numpy.
Para poder comparar secuencias largas necesitamos un algoritmo muy eficiente, y éste debe programarse en C o Rust.
A continuación tienes el enlace a un algoritmo "Python" implementado en C: Levenshtein.
Utiliza esta librería y compara el rendimiento respecto tu solución.
Alineamiento de secuencias
Como ya hemos explicado en la introducción, la comparación de secuencias es una herramienta fundamental para identificar mutaciones que se han producido durante la evolución y han provocado divergencias de las secuencias de ambas proteínas que estamos estudiando.
Pero antes de evaluar la similitud de dos secuencias, normalmente se comienza por encontrar una alineamiento plausible entre ellas teniendo en cuenta que las secuencias en evolución acumulan inserciones y supresiones, así como sustituciones.
Sustitución
Una sustitución en una secuencia se producen cuando una mutación hace que un elemento de la secuencia se cambie en el de otro.
Por ejemplo, a continuación tienes una mutación por sustitución a nivel de ácido nucleico en el sexto elemento de la secuencia (U → A):
GUC UCU → GUC UCA = Valina-Serina → Valina-SerinaPuedes ver que esta mutación a nivel de ADN no tiene traducción a nivel de proteína porque el segundo codón está codificando el propio aminoácido.
En cambio, si la mutación por sustitución hubiera estado en el quinto aminoácido tal y como se muestra a continuación, también habría habido un cambio a nivel de aminoácido (C → A):
GUC UCU → GUC UAU = Valina-Serina → Valina-TirosinaInserción y supresión
Una inserción, o una supresión, se produce cuando se añade un nuevo elemento a la secuencia o se elimina un elemento de la secuencia.
Para mostrar este hecho introducimos un vacío (gap) en la secuencia mediante un guión para mostrar la mutación que ha habido cuando alineamos un par de secuencias.
Por ejemplo, nuestra secuencia se podría ver alterada por una supresión: GUC CU.
Para alinear las secuencias GUC UCU y GUC CU introducimos un -
GUC UCU
GUC -CUEn el caso de una inserción, por ejemplo de una A en el tercer aminoàcido del primer codón, también introducimos un -:
GU- CUC U
GUA CUC UAlineamiento
El alineamiento de secuencias es el proceso de ordenar los caracteres de un par de secuencias de forma que se maximice el número de caracteres coincidentes.
Por ejemplo, imagínate que tenemos estas dos secuencias:
GCGTAACACGTGCG
| | | |
ACAACCCGTGCGACLas barras verticales "|representan nucleótidos coincidentes, y puedes ver que hay poca coincidencia, sólo un 24%.
También que la distancia de edición es 6 si ejecutas tu código python:
C 13 10 8 9 9 7 7 7 6 6 5 4 5 6 6
G 14 11 9 9 10 8 8 8 6 7 5 5 4 5 6
Distance: 6Si introducimos mutaciones de inserción y supresión (“gaps”) de forma que las secuencias a alinear sean las que se muestran a continuación, podemos maximizar el número de aminoácidos coincidentes:
GCGTAACACGTGCG--
| ||| ||||||
AC--AACCCGTGCGACAl introducir los huecos (gaps), la coincidencia ha subido al 56%.
También puedes verificar que la distancia de edición es la misma que antes, 6, que es el número de ácidos nucleicos no coincidentes.
- 15 12 10 9 9 10 11 9 9 9 7 7 6 6 5 5 6
- 16 13 11 9 9 10 11 10 10 10 8 8 7 7 6 6 6
Distancia: 6Puedes ver que el alineamiento de secuencias no ha modificado la distancia de edición, sólo ha alineado las secuencias de forma que el número de aminoácidos coincidentes sea el mayor posible.
Alineamiento global
Uno de los primeros algoritmos para alinear dos secuencias de proteínas fue descrito por Needleman y Wunsch (1970).
Este algoritmo es importante porque produce una alineación óptima de proteínas o secuencias de ADN con la introducción de huecos (gaps).
Matriz
Lo primero que debemos hacer es crear la matriz M tal y como hemos hecho con el algoritmo de la distancia de edición.
La diferencia respecto a la distancia de Levenshtein es que los números iniciales de la primera fila y columna de la matriz son negativos, y el incremento numérico es función de la variable de penalización por brecha.
De esta forma podemos modificar la penalización por gap de forma flexible.
import numpy as np
from util import print_matrix
xs = "GCGTAACACGTGCG"
ys = "ACAACCCGTGCGAC"
gap = -1
#####
M = np.zeros((len(xs) + 1, len(ys) + 1), np.dtype("int8"))
M[:, 0] = np.linspace(0, len(xs) * gap, len(xs) + 1)
M[0, :] = np.linspace(0, len(ys) * gap, len(ys) + 1)
print_matrix(xs, ys, M)Puedes ver el resultado con un brecha de valor -2:
A C A A C C C G T G C G A C
0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14
G -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
T -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A -8 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C -9 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G -10 0 0 0 0 0 0 0 0 0 0 0 0 0 0
T -11 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G -12 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C -13 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G -14 0 0 0 0 0 0 0 0 0 0 0 0 0 0La penalización por gap es necesaria porque introduciendo gaps podemos alinear cualquier secuencia, y debemos encontrar el mejor alineamiento con el menor número de gaps.
Llenar la matriz
Desde la parte superior izquierda de la matriz procedemos a llenar la matriz de izquierda a derecha y de arriba abajo como hemos hecho con la distancia de Levenshtein.
Pero la diferencia respecto a la distancia de edición es que el algoritmo para calcular el valor de M[i][j] es diferente:
for i in range(1, len(xs)+1):
for j in range(1, len(ys)+1):
M[i, j] = max(
M[i - 1][j - 1] + (match if xs[i-1] == ys[j-1] else mismatch),
M[i-1][j] + gap,
M[i][j-1] + gap)Puedes ver que en:
- En lugar del valor mínimo utilizamos el valor máximo de las celdas adyacentes.
- En caso de que los elementos
xs[i-1]yys[j-1]sean diferentes, reducimos el valor de la celda con el valor de la variable falta de coincidencia.
A continuación puedes ver el resultado de llenar la matriz con los valores gap = -1, match = 1 y mismatch = -2:
A C A A C C C G T G C G A C
0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14
G -1 -2 -3 -4 -5 -6 -7 -8 -6 -7 -8 -9 -10 -11 -12
C -2 -3 -1 -2 -3 -4 -5 -6 -7 -8 -9 -7 -8 -9 -10
G -3 -4 -2 -3 -4 -5 -6 -7 -5 -6 -7 -8 -6 -7 -8
T -4 -5 -3 -4 -5 -6 -7 -8 -6 -4 -5 -6 -7 -8 -9
A -5 -3 -4 -2 -3 -4 -5 -6 -7 -5 -6 -7 -8 -6 -7
A -6 -4 -5 -3 -1 -2 -3 -4 -5 -6 -7 -8 -9 -7 -8
C -7 -5 -3 -4 -2 0 -1 -2 -3 -4 -5 -6 -7 -8 -6
A -8 -6 -4 -2 -3 -1 -2 -3 -4 -5 -6 -7 -8 -6 -7
C -9 -7 -5 -3 -4 -2 0 -1 -2 -3 -4 -5 -6 -7 -5
G -10 -8 -6 -4 -5 -3 -1 -2 0 -1 -2 -3 -4 -5 -6
T -11 -9 -7 -5 -6 -4 -2 -3 -1 1 0 -1 -2 -3 -4
G -12 -10 -8 -6 -7 -5 -3 -4 -2 0 2 1 0 -1 -2
C -13 -11 -9 -7 -8 -6 -4 -2 -3 -1 1 3 2 1 0
G -14 -12 -10 -8 -9 -7 -5 -3 -1 -2 0 2 4 3 2 Trace-back
El algoritmo de Needleman-Wunsch tiene como objetivo devolver las secuencias alineadas de tal forma que el número de elementos coincidentes sea el mayor posible.
El algoritmo para hacerlo consiste en empezar por la última celda que se ha computado y hacer un camino inverso de forma recursiva hasta llegar a la celda inicial:
-
De las tres celdas adyacentes de las que se derivó el valor de la celda, debes escoger la celda que tiene el valor máximo.
-
Si hay más de una celda que tiene el mismo valor máximo, la celda
M[i-1][j-1]tiene preferencia porque siempre es el camino más recto. -
Si debemos escoger entre
M[i-1][j]iM[i1][j-1]deberíamos devolver los dos caminos, pero en nuestro caso siempre escogemosM[i-1][j]para simplificar el algoritmo y que sólo devuelva un resultado.
A C A A C C C G T G C G A C
0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14
G -1 *-2 -3 -4 -5 -6 -7 -8 -6 -7 -8 -9 -10 -11 -12
C -2 -3 *-1 -2 -3 -4 -5 -6 -7 -8 -9 -7 -8 -9 -10
G -3 -4 *-2 -3 -4 -5 -6 -7 -5 -6 -7 -8 -6 -7 -8
T -4 -5 *-3 -4 -5 -6 -7 -8 -6 -4 -5 -6 -7 -8 -9
A -5 -3 -4 *-2 -3 -4 -5 -6 -7 -5 -6 -7 -8 -6 -7
A -6 -4 -5 -3 *-1 -2 -3 -4 -5 -6 -7 -8 -9 -7 -8
C -7 -5 -3 -4 -2 *0 -1 -2 -3 -4 -5 -6 -7 -8 -6
A -8 -6 -4 -2 -3 *-1 -2 -3 -4 -5 -6 -7 -8 -6 -7
C -9 -7 -5 -3 -4 -2 *0 *-1 -2 -3 -4 -5 -6 -7 -5
G -10 -8 -6 -4 -5 -3 -1 -2 *0 -1 -2 -3 -4 -5 -6
T -11 -9 -7 -5 -6 -4 -2 -3 -1 *1 0 -1 -2 -3 -4
G -12 -10 -8 -6 -7 -5 -3 -4 -2 0 *2 1 0 -1 -2
C -13 -11 -9 -7 -8 -6 -4 -2 -3 -1 1 *3 2 1 0
G -14 -12 -10 -8 -9 -7 -5 -3 -1 -2 0 2 *4 *3 *2 A continuación tienes el código:
i, j = len(xs), len(ys)
rxs, rys = [], []
while i > 0 or j > 0:
v = max(M[i-1, j-1], M[i-1][j], M[i][j-1])
if M[i-1, j-1] == v:
rxs.append(xs[i - 1])
rys.append(ys[j - 1])
i -= 1
j -= 1
elif M[i-1, j] == v:
rxs.append(xs[i - 1])
rys.append("-")
i -= 1
else:
rxs.append("-")
rys.append(ys[j - 1])
j -= 1
# Reverse the strings.
rxs = "".join(rxs)[::-1]
rys = "".join(rys)[::-1]
print("\n".join([rxs, rys]))Puedes ver que en función de qué celda tenga el valor máximo:
M[i-1][j-1] |
Aceptamos los dos aminoácidos porque son coincidentes. |
M[i-1][j] |
Introducimos un gap en la secuencia ys |
M[i][j-i] |
Introducimos un gap en la secuencia xs |
El resultado del algoritmo es óptimo, aunque no es tan bueno como el que hemos mostrado al principio porque aunque alinea el mismo número de aminoácidos lo hace introduciendo dos gaps de más:
Actividad
1.- Modifica el algoritmo para que indique el porcentaje inicial de coincidencias, el porcentaje final de coincidencias y el porcentaje de gaps introducidos.
2.- Modifica el algoritmo para que obtenga ambas secuencias a comparar de Entrez o Uniprot. Que no sean excesivamente largas o habrá problemas 😂
3.- La universidad McGill ha creado un juego online para ayudarte a encontrar los alineamientos óptimos entre 2 y más secuencias: PHYLO.
Mira com funciona el juego, y piensa como lo implementarías con TypeScript
Alineamiento local
La versión local del algoritmo de alineación de secuencias se desarrolló a principios de los años ochenta, y se conoce como algoritmo de Smith-Waterman (1981).
Hasta ahora hemos supuesto que sabemos qué secuencias queremos alinear y que busquemos la mejor concordancia entre ellas de un extremo a otro.
Una situación mucho más habitual es cuando buscamos la mejor alineación entre subsecuencias de xs y ys.
Esto ocurre, por ejemplo, cuando se sospecha que dos secuencias de proteínas pueden compartir un dominio común, o cuando se comparan secciones extendidas de la secuencia de ADN genómico.
También suele ser la forma más sensible de detectar similitudes cuando se comparan dos secuencias muy divergentes, incluso cuando pueden tener un origen evolutivo compartido a lo largo de toda su longitud.
Esto se debe a que normalmente en estos casos sólo una parte de la secuencia ha estado bajo una selección lo suficientemente fuerte como para preservar la similitud detectable: el resto de la secuencia ha acumulado tanto "ruido" a través de la mutación que ya no se puede alinear.
El alineamiento con mejor puntuación de las subsecuencias de xs e ys se llama mejor alineamiento local.
El algoritmo para encontrar alineaciones locales óptimas está muy relacionado con lo descrito en la sección anterior para alineaciones globales, con dos diferencias importantes.
En primer lugar, a cada celda de la tabla se añade una posibilidad adicional que permite a M[i][j] tomar el valor 0 si todas las demás opciones tienen un valor inferior a 0
import numpy as np
from util import print_matrix
xs = "AAUGCCAUUGACGG"
ys = "CAGCCUCGCUUAG"
gap = -4
match = 3
mismatch = -1
#####
M = np.zeros((len(xs)+1, len(ys) + 1), np.dtype("int8"))
for i in range(1, len(xs) + 1):
for j in range(1, len(ys) + 1):
M[i, j] = max(
0,
M[i-1][j-1] + (match if xs[i-1] == ys[j-1] else mismatch),
M[i-1][j] + gap,
M[i][j-1] + gap
)Tomar la opción 0 corresponde a iniciar un nuevo alineamiento. Si el mejor alineamiento en ese momento converge a una puntuación negativa, es mejor empezar un nuevo alineamiento que continuar con el antiguo.
Por tanto, una consecuencia de este 0 es que la primera fila y columna tienen el número 0, a diferencia del alineamiento global que se llena con los valores -i * gap y -j * gap.
C A G C C U C G C U U A G
0 0 0 0 0 0 0 0 0 0 0 0 0 0
A 0 0 3 0 0 0 0 0 0 0 0 0 3 0
A 0 0 3 2 0 0 0 0 0 0 0 0 3 2
U 0 0 0 2 1 0 3 0 0 0 3 3 0 2
G 0 0 0 3 1 0 0 2 3 0 0 2 2 3
C 0 3 0 0 6 4 0 3 1 6 2 0 1 1
C 0 3 2 0 3 9 5 3 2 4 5 1 0 0
A 0 0 6 2 0 5 8 4 2 1 3 4 4 0
U 0 0 2 5 1 1 8 7 3 1 4 6 3 3
U 0 0 0 1 4 0 4 7 6 2 4 7 5 2
G 0 0 0 3 0 3 0 3 10 6 2 3 6 8
A 0 0 3 0 2 0 2 0 6 9 5 1 6 5
C 0 3 0 2 3 5 1 5 2 9 8 4 2 5
G 0 0 2 3 1 2 4 1 8 5 8 7 3 5
G 0 0 0 5 2 0 1 3 4 7 4 7 6 6 El segundo cambio es que ahora una alineación puede terminar en cualquier parte de la matriz, por lo que en vez de tomar el valor de la última celda que se ha llenado como punto de partida para encontrar la solución, debemos buscar la celda que tiene el valor más alto de la matriz para empezar el trazado desde allí.
max_score = 0
max_index = (1, 1)
for i in range(1, len(xs) + 1):
for j in range(1, len(ys) + 1):
score = M[i, j]
if score > max_score:
max_score = score
max_index = (i, j)El "trace-back" termina cuando encontramos una celda con un valor inferior a match, que corresponde al inicio del alineamiento:
i, j = max_index
rxs, rys = [], []
rxs.append(xs[i - 1])
rys.append(ys[j - 1])
while True:
v = max(M[i-1, j-1], M[i-1][j], M[i][j-1])
if v < match:
break
if M[i-1, j-1] == v:
rxs.append(xs[i - 2])
rys.append(ys[j - 2])
i -= 1
j -= 1
elif M[i-1, j] == v:
rxs.append(xs[i - 2])
rys.append("-")
i -= 1
else:
rxs.append("-")
rys.append(ys[j - 2])
j -= 1
# Reverse the strings.
rxs = "".join(rxs)[::-1]
rys = "".join(rys)[::-1]Aquí tenemos el resultado del alineamiento local:
GCCAUUG
CCG-UCGEn este caso, el alineamiento local es un subconjunto del alineamiento global, pero no siempre es así.
A continuació puedes ver la diferencia entre una alineamineto global y úno local del mismo par de secuencias:
Actividad
Aquí tienes un ejemplo de aplicación web que muestra de forma interactiva cómo funciona el algoritmo de Smith-Waterman: Teaching - Smith-Waterman.
Mira los scripts que se utilizan para ver un ejemplo de cómo podrías replicar esta página: Source code
Afinidad de brecha
Gotoh formuló la versión eficiente de coste de la brecha (“gap”) afín que se utiliza normalmente.
Hasta ahora sólo hemos considerado el modelo de gap más sencillo, en el que la puntuación del gap es múltiplo de la longitud de los elementos que lo componen.
Este tipo de esquema de puntuación no es ideal para secuencias biológicas: la mayoría de los gaps no son de un solo elemento, y estamos penalizando por la extensión de los elementos que forman el gap en lugar del número de gaps que estamos introduciendo.
Por ejemplo, este alineamiento hecho a mano sólo tiene dos gaps y la penalización es de 4.
GCGTAACACGTGCG--
AC--AACCCGTGCGACEn cambio el alineamiento hecho por nuestro código python tiene 4 gaps, el doble de gaps, pero sólo penaliza 6.
GCGTAACAC-GTGCG--
AC--AAC-CCGTGCGACEn lugar de utilizar un valor fijo de penalización por gap podemos utilizar una función que calcule este valor en función de algunos parámetros como en este ejemplo:
def gap_penalty(k):
a, b = 5, 2
return -(a + b*k)La función recibe como argumento el número de elementos que tiene el gap (k) y devuelve un valor de penalización de gap en función de dos variables:
-
El valor de la variable
aes la penalización por crear el gap y es independiente del número de elementos que componen el gap. -
El valor de la variable
bes la penalización por el número de elementos que componen el gap y es función del número de elementos que componen el gap (la variable k).
Los valores a i b que hemos escogido son arbitrarios y sólo sirven para mostrar cómo funciona la función:
for i in range (1,6):
print(f"{gap_penalty(i)}", end= " ")
print()Puedes ver que en penalización entre un gap de 1 elemento y uno de 6 elementos sólo es el doble:
-7 -9 -11 -13 -15Biopython
El módulo Bio.Align tiene la clase Alineador PairwiseAligner para realizar alineamientos por parejas globales y locales utilizando los algoritmos de Needleman-Wunsch, Smith-Waterman, Gotoh (de tres estados) y Waterman-Smith-Beyer, con numerosas opciones para cambiar los parámetros de alineamiento.
Introducción
Para generar alineamientos por parejas, primero debes crear un objeto PairwiseAligner:
from Bio import Align
aligner = Align.PairwiseAligner()El objeto PairwiseAligner almacena los parámetros de alineamiento que se utilizarán para los alineamientos por parejas.
Estos atributos se pueden configurar al construir el objeto:
aligner = Align.PairwiseAligner(match_score=1.0)o después de crear el objeto:
aligner.match_score = 1.0score
El método score() permite calcular la puntuación de alineamiento entre dos secuencias:
from Bio import Align
target = "GCGTAACACGTGCG"
query = "ACAACCCGTGCGAC"
aligner = Align.PairwiseAligner()
score = aligner.score(target, query)
assert score == 10alignment
El método aligment() devuelve todos los alineamientos encontrados.
Puedes iterar sobre los objetos de alineamiento e imprimirlos por pantalla para ver los alineamientos:
aligner = Align.PairwiseAligner()
aligments = aligner.align(target, query)
for aligment in aligments:
print(aligment)Puedes ver que hay más de una posible alineación:
target 0 G-CGTAACAC-GTGCG-- 14
0 --|--|||-|-|||||-- 18
query 0 -AC--AAC-CCGTGCGAC 14
target 0 -GCGTAACAC-GTGCG-- 14
0 --|--|||-|-|||||-- 18
query 0 A-C--AAC-CCGTGCGAC 14Cada alineamiento almacena la puntuación de alineamiento así como los punteros en las secuencias que se han alineado:
aligner = Align.PairwiseAligner()
aligment = aligner.align(target, query)[0]
print(f"{aligment.score}: {aligment.target} -> {aligment.query}")10.0: GCGTAACACGTGCG -> ACAACCCGTGCGACcoordinates
El alineamiento también almacena las coordenadas del alineamiento de las secuencias:
aligner = Align.PairwiseAligner()
aligment = aligner.align(target, query)[0]
print(aligment.coordinates)
print(aligment)Las coordenadas muestran que el alineamiento consta de 10 bloques:
[[ 0 1 1 2 4 7 8 9 9 14 14]
[ 0 0 1 2 2 5 5 6 7 12 14]]
target 0 G-CGTAACAC-GTGCG-- 14
--|--|||-|-|||||-- 18
query 0 -AC--AAC-CCGTGCGAC 14
0:1 1:1 1:2 2:4 4:7 7:8 8:9 9:9 9:14 14:14
0:0 0:1 1:2 2:2 2:5 5:5 5:6 6:7 7:12 12:14Cada bloque es un "slice" de la secuencia, y un slice como [2:2] equivale a "" que es lo mismo que un gap.
length
La longitud de la alineación se define como el número de columnas de la alineación tal y como se imprime. Esto es igual a la suma del número de coincidencias, el número de no coincidencias y la longitud total de los vacíos en las secuencias target y query:
target = "GCGTAACACGTGCG"
query = "ACAACCCGTGCGAC"
aligner = Align.PairwiseAligner()
aligments = aligner.align(target, query)
print(f"{len(aligments)}:", end = " ")
for aligment in aligments:
print(aligment.length, end = " ")
print()Puedes ver que tenemos 15 alineamientos posibles, con length entre 18 y 16:
15: 18 18 17 18 18 17 18 18 17 18 18 17 17 17 16Ten en cuenta que diferentes alineamientos pueden tener las mismas subsecuencias alineadas entre sí.
En particular, esto puede ocurrir si los alineamientos difieren entre sí sólo en cuanto a la colocación de los huecos tal y como puedes ver en este ejemplo:
from Bio import Align
aligner = Align.PairwiseAligner(mode="global", mismatch_score = -10)
alignments = aligner.align("AAACAAA", "AAAGAAA")
print(alignments[0])
print(alignments[1])target 0 AAAC-AAA 7
0 |||--||| 8
query 0 AAA-GAAA 7
target 0 AAA-CAAA 7
0 |||--||| 8
query 0 AAAG-AAA 7PairwiseAligner
El objeto PairwiseAligner almacena todos los parámetros de alineamiento que se utilizarán para las alineaciones por pares.
Para ver los valores de todos los parámetros puedes realizar un print del objeto:
from Bio import Align
aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
print(aligner)Puedes ver que los valores puntuación de match_score y mode corresponden a los que se han pasado como argumentos al constructor del objeto:
Pairwise sequence aligner with parameters
wildcard: None
match_score: 1.000000
mismatch_score: 0.000000
mode: local
target_internal_open_gap_score: 0.000000
...El atributo mode puede tener el valor "global" o "local", y especifica un alineamiento por parejas global o local respectivamente.
En función de los parámetros de puntuación de "gap" y del modo, el objeto PairwiseAligner elige automáticamente el algoritmo adecuado para computar el alineamiento correspondiente.
El atributo algorithm sólo tiene acceso de lectura y no puede modificarse.
aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
assert aligner.algorithm == "Smith-Waterman"
aligner.mode = "global"
assert aligner.algorithm == "Needleman-Wunsch"El objeto PairwiseAligner también almacena el atributo de precisión `epsilonp que se utiliza durante el alineamiento.
El valor predeterminado es igual a 10−6:
aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
assert aligner.epsilon == 1e-6Dos puntuaciones se considerarán iguales a efectos del alineamiento si la diferencia absoluta entre ellas es inferior a valor del atributo epsilon.
Puntuación de sustitución
Las puntuaciones de sustitución definen el valor a añadir a la puntuación total cuando dos letras (nucleótidos o aminoácidos) están alineadas entre sí.
Las alineaciones de la secuencia de nucleótidos se basan normalmente en puntuaciones de concordancia (match) y desajuste (dismatch). Por ejemplo, de forma predeterminada, Blast utiliza una puntuación de coincidencia de +1 y una puntuación de desajuste de -2 para alineaciones de nucleótidos por megablast, con una penalización de gap de 2,5.
Las puntuaciones de coincidencia y desajuste se pueden especificar estableciendo los atributos de coincidencia y desajuste del objeto PairwiseAligner:
from Bio import Align
aligner = Align.PairwiseAligner()
assert aligner.match_score == 1.0
assert aligner.mismatch_score == 0.0
score = aligner.score("ACGT", "ACAT")
assert score == 3.0
aligner.match_score = 1.0
aligner.mismatch_score = -2.0
aligner.gap_score = -2.5
score = aligner.score("ACGT", "ACAT")
assert score == 1.0Gaps
gap scores predefinidos
De forma predeterminada, un objeto PairwiseAligner se inicia con una puntuación de match de +1,0, una puntuación de dismatch de 0,0 y todas las puntuaciones de gap son iguales a 0,0.
Aunque esto tiene la ventaja de ser un esquema de puntuación sencillo, por lo general no da el mejor rendimiento. En lugar de esto, puedes utilizar el argumento de puntuación para seleccionar un esquema de puntuación predefinido al inicializar un objeto PairwiseAligner.
Los esquemas de puntuación proporcionados son:
blastnimegablast, que son adecuados para alineaciones de nucleótidos.blastpque es adecuado para alineaciones de proteínas.
La selección de estos esquemas de puntuación inicializará el objeto PairwiseAligner con los parámetros de puntuación predeterminados utilizados por BLASTN, MegaBLAST y BLASTP, respectivamente.
from Bio import Align
aligner = Align.PairwiseAligner(scoring="blastn")Affine gap scores
Las puntuaciones del gap afín se definen por una puntuación para abrir una gap y una puntuación para ampliar una gap existente:
gap_score = open_gap_score + (n−1) × extend_gap_score,donde n es la longitud del gap.
Biopython permite un control preciso del esquema de puntuación de huecos mediante doce atributos del objeto PairwiseAligner por pares tal y como puedes ver en este ejemplo:
| target | query | score |
|---|---|---|
| A | - | query_left_open_gap_score |
TODO revisar
gap function
Para un control aún más preciso de las puntuaciones de gap, puedes especificar una función de puntuación de gap.
Por ejemplo, esta función no permite un vacío después de dos nucleótidos en la secuencia de consulta:
from Bio import Align
def gap_score(start, length):
if start == 2:
return -1000
else:
return -1 * length
aligner = Align.PairwiseAligner()
print(aligner.align("AAACTT", "AAATT")[0])
aligner.query_gap_score = gap_score
print(aligner.align("AACGGTT", "AATT")[0])Puedes ver que la función no permite introducir un gap después de AAA como se hace en el primer alineamiento:
target 0 AAACTT 6
0 |||-|| 6
query 0 AAA-TT 5
target 0 AACGGTT 7
0 -|--.|| 7
query 0 -A--ATT 4Obtención de información sobre la alineación
Google Docs
El documento continua en Google Docs: Alineament de seqüències
TODO
Important!. Falta completar aspectes com extended gap, odds-log ratio, matrius blosum,etc.