Alineament de seqüències

#Introducció

La forma més fiable de determinar l’estructura o funció d’una molècula biològica és mitjançant l’experimentació directa. No obstant això, és molt més fàcil obtenir la seqüència d’ADN del gen corresponent a un ARN o una proteïna, que determinar experimentalment la seva funció o estructura.

Per tant, part del repte és simplement organitzar, classificar i analitzar la immensa riquesa de dades de seqüències, i desenvolupar mètodes computacionals que puguin inferir informació biològica només a partir de la seqüència.

#Seqüències

Les noves seqüències s’adapten a partir de seqüències preexistents en lloc d’inventar-se de zero, la qual cosa ens permet analitzar les seqüències de forma computacional per reconèixer una similitud significativa entre una seqüència nova i una seqüència sobre la qual ja se sap alguna cosa.

Quan fem això podem transferir informació sobre l’estructura i/o la funció a la nova seqüència. Diem que les dues seqüències relacionades són homòlogues i que estem transferint informació per homologia.

Dues seqüències són homòlogues si es deriven d’una seqüència comuna, és a dir, tenen un ancestre comú.

L’homologia no és quantificable, les seqüències són o no són homòlogues.

Per tant, una de les preguntes més bàsiques sobre un gen o proteïna és si està relacionat amb algun altre gen o proteïna:

• La relació de dues proteïnes a nivell de seqüència suggereix que són homòlogues.

• La relació també suggereix que poden tenir funcions comunes.

Mitjançant l’anàlisi de moltes seqüències d’ADN i proteïnes, és possible identificar dominis o motius que es comparteixen entre un grup de molècules.

Aquestes anàlisis de la relació de proteïnes i gens s’aconsegueixen alineant seqüències.

La majoria dels problemes en l’anàlisi de seqüències computacionals són essencialment estadístics. Les forces evolutives estocàstiques actuen sobre els genomes. Discernir similituds significatives entre seqüències antigament divergents enmig d’un caos de mutació aleatòria, selecció natural i deriva genètica presenta seriosos problemes de senyal al soroll. Molts dels mètodes d’anàlisi més potents disponibles utilitzen la teoria de la probabilitat.

#Alineament de proteïnes

A l’hora de comparar proteïnes, podem comparar la seqüència d’aminoàcids que la componen, o comparar els gens que codifiquen la seqüència d’aminoàcids.

Les seqüències d’aminoàcids són més informatives perquè:

• Molts canvis en una seqüència d’ADN (especialment en la tercera posició d’un codó) no canvien l’aminoàcid que s’especifica.

• Molts aminoàcids comparteixen propietats biofísiques relacionades (per exemple, la lisina i l’arginina són ambdós aminoàcids bàsics).

El resultat és que les comparacions de seqüències de proteïnes poden identificar seqüències homòlogues mentre que les comparacions de seqüències d’ADN corresponents no poden.

Per tant, quan s’analitza una seqüència de codificació de nucleòtids sovint és preferible estudiar la seva proteïna traduïda.

A Blast veurem que podem moure’ns fàcilment entre els mons de l’ADN i les proteïnes.

Per exemple, l’eina TBLASTN del lloc web NCBI BLAST permet cercar proteïnes relacionades derivades d’una base de dades d’ADN amb una seqüència de proteïnes.

Aquesta opció de consulta s’aconsegueix traduint cada seqüència d’ADN a les sis proteïnes que potencialment codifica.

No obstant això, tot el que hem explicat fins ara, en molts casos és adequat comparar seqüències d’ADN.

Aquesta comparació pot ser important per confirmar la identitat d’una seqüència d’ADN en una cerca de bases de dades, per cercar polimorfismes, per analitzar la identitat d’un fragment d’ADNc clonat o per comparar regions reguladores, per posar alguns exemples.

#Homologia

Dues seqüències són homòlogues si comparteixen un avantpassat evolutiu comú.

No hi ha graus d’homologia; les seqüències són homòlogues o no ho són.

Les proteïnes homòlogues gairebé sempre comparteixen una estructura tridimensional significativament relacionada. Un exemple d’homologia és el de la mioglobina humana (NP_005359.1) i la beta globina (NP_000509.1), dues proteïnes que estan relacionades llunyanament, però de manera significativa.

Es creu que la mioglobina i les cadenes d’hemoglobina (alfa, beta i altres) van divergir fa uns 450 milions d’anys, prop del moment en què els llinatges dels humans i dels peixos cartilaginosos van divergir.

La mioglobina i la beta globina tenen estructures molt semblants, tal com determina la cristal·lografia de raigs X i com pots veure en aquesta imatge:

Mioglobina humana (3RGK)	Beta globina humana (subunitat de 2H35)

Aquestes proteïnes són homòlogues (descendeixen d’un avantpassat comú) i comparteixen estructures tridimensionals molt similars.

Tanmateix, l’alineament per parelles de les seqüències d’aminoàcids d’aquestes proteïnes revela que comparteixen una identitat d’aminoàcids molt limitada.

#Identitat

Quan dues seqüències són homòlogues, les seves seqüències d’aminoàcids o de nucleòtids solen compartir una identitat significativa.

Si bé l’homologia és una inferència qualitativa (les seqüències són o no són homòlogues), la identitat i la semblança són quantitats que descriuen la relació de les seqüències.

En particular, dues molècules poden ser homòlogues sense compartir una identitat estadísticament significativa d’aminoàcids (o de nucleòtids).

A la família de la globina, tots els membres són homòlegs, però alguns tenen seqüències que han divergit tant que no comparteixen cap identitat de seqüència reconeixible.

Per exemple:

• La beta globina humana i la neuroglobina humana només comparteixen un 22% d’identitat d’aminoàcids.
• Les cadenes de globina individuals comparteixen la mateixa forma general que la mioglobina, tot i que les proteïnes de mioglobina i d’alfa globina només comparteixen un 26% d’identitat d’aminoàcids.

En general, les estructures tridimensionals divergeixen molt més lentament que la identitat de la seqüència d’aminoàcids entre dues proteïnes.

I és l’estructura la que defineix la funció de la proteïna!

Reconèixer aquest tipus d’homologia és un problema bioinformàtic especialment complex.

#Ortologia

Les proteïnes que són homòlogues poden ser ortòlogues o paràlogues.

Els ortòlegs són seqüències homòlogues d’espècies diferents que tenen un gen avantpassat comú i que han divergit per especiació.

La història del gen reflecteix la història de l’espècie, per això els gens s’anomenen ortòlegs (orto = exacte). Se suposa que els ortòlegs tenen funcions biològiques similars, i per tant pots deduir la funció d’una seqüència a partir d’una altra seqüència homòloga.

Per exemple, els humans i els rosegadors van divergir fa uns 90 milions d’anys (MYA), moment en què un únic gen de la mioglobina ancestral va divergir per especiació, i en tots dos casos la funció del gen és transportar oxigen a les cèl·lules musculars.

A continuació tens un dibuix amb l’arbre dels ortòlegs de la mioglobina:

#Paralogia

Els paràlegs són seqüències homòlogues sorgides per un mecanisme com la duplicació de gens.

En aquest cas l’homologia és el resultat de la duplicació gènica, de manera que totes dues còpies han baixat l’una al costat de l’altra durant la història d’un organisme.

Els gens s’anomenen paràlegs (para = en paral·lel) perquè evolucionen dins de la mateixa espècie.

Per exemple, la globina alfa 1 humana (NP_000549.1) és paràloga de la globina alfa 2 (NP_000508.1). De fet, aquestes dues proteïnes comparteixen una identitat d’aminoàcids del 100%.

La globina alfa 1 i la globina beta humanes també són paràlogues, com totes les proteïnes humanes que es mostren a continuació i que són membres de la família globina:

Totes les globines tenen propietats diferents, incloent-hi la distribució regional al cos, el moment de desenvolupament de l’expressió gènica i l’abundància. Però encara que totes tenen funcions diferents, aquestes funcions estan relacionades com a proteïnes transportadores d’oxigen.

Per saber una mica més sobre aquests arbres i l’evolució pots consultar aquest vídeo:

#Similitud

Mira “/python/dynamic-programming”

A primer cop d’ull, decidir que dues seqüències biològiques són similars no és diferent de decidir que dues cadenes de text són similars.

Els algorismes de semblança de text s’utilitzen per buscar la paraula més semblant d’un diccionari quan una paraula no s’hi troba perquè està mal escrita. És el que coneixes com a corrector ortogràfic bàsic.

Per tant, no t’ha d’estranyar que els algorismes d’anàlisi de seqüències s’assemblin molt als algorismes de semblança de text.

La distància de Levenshtein (o distància d’edició) és un algorisme de 1965 dissenyat per comparar dues paraules, que es va aplicar posteriorment a seqüències biològiques.

Aquest algorisme es basa en el fet que, per transformar una paraula en una altra, només calen tres operacions bàsiques: inserir, suprimir o substituir un caràcter.

Cada vegada que s’aplica una d’aquestes operacions, la distància s’incrementa 1 punt: com més punts, més gran és la distància d’edició de les dues paraules i menys similars són.

#Matriu de distàncies

La distància d’edició és un algorisme lineal que utilitza una matriu (un array de dues dimensions) per guardar els resultats parcials que permeten anar construint el resultat final.

Donades dues seqüències d’àcids nucleics xs i ys, la matriu de distàncies conté les distàncies entre tots els prefixos de la seqüència xs i tots els prefixos de la seqüència ys.

Per començar, has d’omplir la matriu amb el resultat del subproblema bàsic a partir del qual anirem construint resultats parcials fins a arribar al resultat final.

Per exemple, si vols comparar les seqüències CGA i AGAT, has de crear una matriu de 4 x 5 amb aquests valors inicials (el símbol - indica un string buit "")

    -  A  G  A  T
 -  0  1  2  3  4
 C  1  ? 
 G  2
 A  3 

Com pots veure en aquest exemple, és evident que la distància d’edició entre:

• “” i “A” és 1, perquè per transformar “” en “A” hem d’inserir una A.
• “” i “AG” és 2, perquè per transformar “” en “AG” hem d’inserir una A i després una G.

A continuació tens l’algorisme inicial implementat 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 executem l’script podem veure la taula resultant:

    -  A  G  A  T
 -  0  1  2  3  4
 C  1   
 G  2
 A  3 

A continuació hem d’anar omplint la matriu a partir d’on hi ha l’interrogant ?, utilitzant els valors que hi ha a les cel·les M[i-1][j-1], M[i-1][j] i M[i][j-1], i resolent subproblemes fins que arribem a trobar la solució del problema.

L’algorisme que hem d’utilitzar és aquest:

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] += 1

Tal com pots veure al codi, la distància d’edició d’una cel·la serà:

• El valor de la distància mínima per arribar fins a aquesta cel·la a partir de la cel·la M[i-1][j-1], M[i-1][j] o M[i][j-1]
• Si s’ha de fer una edició perquè la lletra de la fila i i la de la columna j no són iguals, incrementem el valor de la cel·la en 1.

Per tant, el segon pas de l’algorisme és

    -  A  G  A  T
 -  0  1  2  3  4
 C  1  1 ? 
 G  2
 A  3 

El resultat d’omplir la matriu serà el que es mostra a continuació, on tenim el resultat de la distància entre tots els prefixos possibles que s’han anat calculant per trobar el resultat 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  2

A continuació pots veure el camí mínim d’edició marcat amb *, ja que tots els altres camins, encara que siguin correctes, comporten moltes més operacions per 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  *2

Per convertir CGA en AGT hem de canviar la C per una A i afegir una T al final.

Són 2 operacions (tal com indica el valor de la cel·la final), i per tant la distància d’edició és 2.

Com ja hem comentat abans, a la taula hi ha calculats tots els resultats dels prefixos.

Per exemple, la distància entre CG i AGA és 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   2

Executant l’script en Python pots veure que el resultat de l’algorisme és la mateixa taula que hem calculat a mà:

#Activitat: Cost computacional

L’algorisme de la distància d’edició es va dissenyar per comparar paraules, i les paraules tenen poques lletres.

Però si comparem seqüències d’ADN, que poden ser molt llargues, el temps que triga a computar el resultat és quadràtic respecte a l’entrada: O(n²).

Això vol dir que si triga 0,00001 segons a calcular la distància d’una entrada de 5x5, trigarà 0,4 segons a calcular el resultat d’una entrada de 1000x1000.

Modifica el fitxer perquè obtingui dues seqüències grans d’Entrez i calcula la distància d’edició.

Fes un profiling de l’script tal com vas fer a “/python/profiling/”

Implementar aquest algorisme en Python és una activitat educativa interessant, encara que fem servir “/python/numpy”.

Per poder comparar seqüències llargues necessitem un algorisme molt eficient, i s’ha de programar en C o Rust.

A continuació tens l’enllaç a un algorisme “Python” implementat en C: Levenshtein.

Utilitza aquesta llibreria i compara’n el rendiment amb la teva solució.

#Alineament de seqüències

Com ja hem explicat a la introducció, la comparació de seqüències és una eina fonamental per identificar mutacions que s’han produït durant l’evolució i han provocat divergències de les seqüències de les dues proteïnes que estem estudiant.

Però abans d’avaluar la similitud de dues seqüències, normalment es comença per trobar un alineament plausible entre elles tenint en compte que les seqüències en evolució acumulen insercions i supressions, així com substitucions.

#Substitució

Una substitució en una seqüència es produeix quan una mutació fa que un element de la seqüència es canviï pel d’un altre.

Per exemple, a continuació tens una mutació per substitució a nivell d’àcid nucleic en el sisè element de la seqüència (U → A):

GUC UCU → GUC UCA   =   Valina-Serina → Valina-Serina

Pots veure que aquesta mutació a nivell d’ADN no té traducció a nivell de proteïna perquè el segon codó codifica el mateix aminoàcid.

En canvi, si la mutació per substitució hagués estat al cinquè aminoàcid tal com es mostra a continuació, també hi hauria hagut un canvi a nivell d’aminoàcid (C → A):

GUC UCU → GUC UAU   =   Valina-Serina → Valina-Tirosina

#Inserció i supressió

Una inserció, o una supressió, es produeix quan s’afegeix un element nou a la seqüència o se n’elimina un.

Per mostrar aquest fet introduïm un forat (gap) a la seqüència mitjançant un guionet, per indicar la mutació que hi ha hagut quan alineem un parell de seqüències.

Per exemple, la nostra seqüència es podria veure alterada per una supressió: GUC CU.

Per alinear les seqüències GUC UCU i GUC CU introduïm un -

GUC UCU
GUC -CU

En el cas d’una inserció, per exemple d’una A en el tercer aminoàcid del primer codó, també introduïm un -:

GU- CUC U
GUA CUC U

#Alineament

L’alineament de seqüències és el procés d’ordenar els caràcters d’un parell de seqüències de manera que es maximitzi el nombre de caràcters coincidents.

Per exemple, imagina’t que tenim aquestes dues seqüències:

GCGTAACACGTGCG
 |    |  | |
ACAACCCGTGCGAC

Les barres verticals “|” representen nucleòtids coincidents, i pots veure que hi ha poca coincidència, només un 24%.

També que la distància d’edició és 6 si executes el teu codi 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: 6

Si introduïm mutacions d’inserció i supressió («gaps») de manera que les seqüències a alinear siguin les que es mostren a continuació, podem maximitzar el nombre d’aminoàcids coincidents:

GCGTAACACGTGCG--
 |  ||| ||||||   
AC--AACCCGTGCGAC

En introduir els forats (gaps), la coincidència ha pujat al 56%.

També pots verificar que la distància d’edició és la mateixa que abans, 6, que és el nombre d’àcids nucleics no coincidents.

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

Pots veure que l’alineament de seqüències no ha modificat la distància d’edició, només ha alineat les seqüències de manera que el nombre d’aminoàcids coincidents sigui el més gran possible.

#Alineament global

Un dels primers algorismes per alinear dues seqüències de proteïnes el van descriure Needleman i Wunsch (1970).

Aquest algorisme és important perquè produeix un alineament òptim de proteïnes o seqüències d’ADN amb la introducció de forats (gaps).

#Matriu

El primer que hem de fer és crear la matriu M tal com hem fet amb l’algorisme de la distància d’edició.

La diferència respecte a la distància de Levenshtein és que els nombres inicials de la primera fila i columna de la matriu són negatius, i l’increment numèric és funció de la variable de penalització per forat.

D’aquesta manera podem modificar la penalització per gap de manera 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)

Pots veure el resultat amb un forat 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   0

La penalització per gap és necessària perquè introduint gaps podem alinear qualsevol seqüència, i hem de trobar el millor alineament amb el mínim nombre de gaps.

#Omplir la matriu

Des de la part superior esquerra de la matriu, procedim a omplir-la d’esquerra a dreta i de dalt a baix com hem fet amb la distància de Levenshtein.

Però la diferència respecte a la distància d’edició és que l’algorisme per calcular el valor de M[i][j] és diferent:

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)

Pots veure que:

• En lloc del valor mínim, fem servir el valor màxim de les cel·les adjacents.
• En cas que els elements xs[i-1] i ys[j-1] siguin diferents, reduïm el valor de la cel·la amb el valor de la variable de manca de coincidència.

A continuació pots veure el resultat d’omplir la matriu amb els valors gap = -1, match = 1 i 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

L’algorisme de Needleman-Wunsch té com a objectiu retornar les seqüències alineades de manera que el nombre d’elements coincidents sigui el més gran possible.

L’algorisme per fer-ho consisteix a començar per l’última cel·la que s’ha computat i fer un camí invers de manera recursiva fins a arribar a la cel·la inicial:

1. De les tres cel·les adjacents de les quals es va derivar el valor de la cel·la, has d’escollir la que té el valor màxim.

2. Si hi ha més d’una cel·la amb el mateix valor màxim, la cel·la M[i-1][j-1] té preferència perquè sempre és el camí més recte.

3. Si hem d’escollir entre M[i-1][j] i M[i1][j-1] hauríem de retornar els dos camins, però en el nostre cas sempre escollim M[i-1][j] per simplificar l’algorisme i que només retorni un resultat.

        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ó tens el codi:

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

# Invertim les cadenes.
rxs = "".join(rxs)[::-1]
rys = "".join(rys)[::-1]

print("\n".join([rxs, rys]))

Pots veure que, en funció de quina cel·la tingui el valor màxim:

El resultat de l’algorisme és òptim, encara que no és tan bo com el que hem mostrat al principi perquè, tot i que alinea el mateix nombre d’aminoàcids, ho fa introduint dos gaps de més:

#Activitat

1.- Modifica l’algorisme perquè indiqui el percentatge inicial de coincidències, el percentatge final de coincidències i el percentatge de gaps introduïts.

2.- Modifica l’algorisme perquè obtingui les dues seqüències a comparar d’Entrez o Uniprot. Que no siguin excessivament llargues o hi haurà problemes 😂

3.- La universitat McGill ha creat un joc en línia per ajudar-te a trobar els alineaments òptims entre 2 i més seqüències: PHYLO.

Mira com funciona el joc, i pensa com l’implementaries amb “/ts/”

#Alineament local

La versió local de l’algorisme d’alineament de seqüències es va desenvolupar a principis dels anys vuitanta, i es coneix com a algorisme de Smith-Waterman (1981).

Fins ara hem suposat que sabem quines seqüències volem alinear i que en busquem la millor concordança d’un extrem a l’altre.

Una situació molt més habitual és quan busquem el millor alineament entre subseqüències de xs i ys.

Això passa, per exemple, quan se sospita que dues seqüències de proteïnes poden compartir un domini comú, o quan es comparen seccions extenses de la seqüència d’ADN genòmic.

També acostuma a ser la manera més sensible de detectar similituds quan es comparen dues seqüències molt divergents, fins i tot quan poden tenir un origen evolutiu compartit al llarg de tota la seva longitud.

Això és així perquè normalment, en aquests casos, només una part de la seqüència ha estat sota una selecció prou forta com per preservar la similitud detectable: la resta de la seqüència ha acumulat tant “soroll” a través de la mutació que ja no es pot alinear.

L’alineament amb millor puntuació de les subseqüències de xs i ys s’anomena millor alineament local.

L’algorisme per trobar alineaments locals òptims està molt relacionat amb el que hem descrit a la secció anterior per als alineaments globals, amb dues diferències importants.

En primer lloc, a cada cel·la de la taula s’hi afegeix una possibilitat addicional que permet que M[i][j] prengui el valor 0 si totes les altres opcions tenen 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
        )

Prendre l’opció 0 correspon a iniciar un alineament nou. Si el millor alineament en aquell moment convergeix a una puntuació negativa, és millor començar un alineament nou que continuar amb l’antic.

Per tant, una conseqüència d’aquest 0 és que la primera fila i la primera columna tenen el nombre 0, a diferència de l’alineament global, que s’omple amb els valors -i * gap i -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 segon canvi és que ara un alineament pot acabar en qualsevol part de la matriu, de manera que en comptes de prendre el valor de l’última cel·la que s’ha omplert com a punt de partida per trobar la solució, hem de buscar la cel·la que té el valor més alt de la matriu per començar el traçat des d’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” acaba quan trobem una cel·la amb un valor inferior a match, que correspon a l’inici de l’alineament:

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

# Invertim les cadenes.
rxs = "".join(rxs)[::-1]
rys = "".join(rys)[::-1]

Aquí tenim el resultat de l’alineament local:

GCCAUUG
CCG-UCG

En aquest cas, l’alineament local és un subconjunt de l’alineament global, però no sempre és així.

A continuació pots veure la diferència entre un alineament global i un de local del mateix parell de seqüències:

#Activitat

Aquí tens un exemple d’aplicació web que mostra de manera interactiva com funciona l’algorisme de Smith-Waterman: Teaching - Smith-Waterman.

Mira els scripts que s’utilitzen per veure un exemple de com podries replicar aquesta pàgina: Source code

#Gap afí

Gotoh va formular la versió eficient del cost del forat («gap») afí que s’utilitza normalment.

Fins ara només hem considerat el model de gap més senzill, en què la puntuació del gap és múltiple de la longitud dels elements que el componen.

Aquest tipus d’esquema de puntuació no és ideal per a seqüències biològiques: la majoria dels gaps no són d’un sol element, i estem penalitzant per l’extensió dels elements que formen el gap en lloc del nombre de gaps que estem introduint.

Per exemple, aquest alineament fet a mà només té dos gaps i la penalització és de 4.

GCGTAACACGTGCG-- 
AC--AACCCGTGCGAC

En canvi, l’alineament fet pel nostre codi python té 4 gaps, el doble de gaps, però només penalitza 6.

GCGTAACAC-GTGCG--   
AC--AAC-CCGTGCGAC

En lloc d’utilitzar un valor fix de penalització per gap podem utilitzar una funció que calculi aquest valor en funció d’alguns paràmetres, com en aquest exemple:

def gap_penalty(k):
    a, b = 5, 2
    return -(a + b*k)

La funció rep com a argument el nombre d’elements que té el gap (k) i retorna un valor de penalització de gap en funció de dues variables:

• El valor de la variable a és la penalització per crear el gap i és independent del nombre d’elements que el componen.

• El valor de la variable b és la penalització pel nombre d’elements que componen el gap i és funció d’aquest nombre d’elements (la variable k).

Els valors a i b que hem escollit són arbitraris i només serveixen per mostrar com funciona la funció:

for i in range (1,6):
    print(f"{gap_penalty(i)}", end= " ")
print()

Pots veure que la penalització entre un gap d’1 element i un de 6 elements només és el doble:

-7 -9 -11 -13 -15

#Biopython

El mòdul Bio.Align té la classe alineadora PairwiseAligner per fer alineaments per parelles globals i locals utilitzant els algorismes de Needleman-Wunsch, Smith-Waterman, Gotoh (de tres estats) i Waterman-Smith-Beyer, amb nombroses opcions per canviar els paràmetres d’alineament.

#Introducció

Per generar alineaments per parelles, primer has de crear un objecte PairwiseAligner:

from Bio import Align
aligner = Align.PairwiseAligner()

L’objecte PairwiseAligner emmagatzema els paràmetres d’alineament que s’utilitzaran per als alineaments per parelles.

Aquests atributs es poden configurar en construir l’objecte:

aligner = Align.PairwiseAligner(match_score=1.0)

o després de crear l’objecte:

aligner.match_score = 1.0

#score

El mètode score() permet calcular la puntuació d’alineament entre dues seqüències:

from Bio import Align

target = "GCGTAACACGTGCG"
query = "ACAACCCGTGCGAC"

aligner = Align.PairwiseAligner()
score = aligner.score(target, query)
assert score == 10

#alignment

El mètode aligment() retorna tots els alineaments trobats.

Pots iterar sobre els objectes d’alineament i imprimir-los per pantalla per veure’ls:

aligner = Align.PairwiseAligner()
aligments = aligner.align(target, query)
for aligment in aligments:
    print(aligment)

Pots veure que hi ha més d’un alineament possible:

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 14

Cada alineament emmagatzema la puntuació d’alineament així com els punters a les seqüències que s’han alineat:

aligner = Align.PairwiseAligner()
aligment = aligner.align(target, query)[0]

print(f"{aligment.score}: {aligment.target} -> {aligment.query}")

10.0: GCGTAACACGTGCG -> ACAACCCGTGCGAC

#coordinates

L’alineament també emmagatzema les coordenades de l’alineament de les seqüències:

aligner = Align.PairwiseAligner()
aligment = aligner.align(target, query)[0]


print(aligment.coordinates)
print(aligment)

Les coordenades mostren que l’alineament consta de 10 blocs:

[[ 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:14

Cada bloc és un “slice” de la seqüència, i un slice com [2:2] equival a “”, que és el mateix que un gap.

#length

La longitud de l’alineament es defineix com el nombre de columnes de l’alineament tal com s’imprimeix. Això és igual a la suma del nombre de coincidències, el nombre de no coincidències i la longitud total dels forats a les seqüències target i 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()

Pots veure que tenim 15 alineaments possibles, amb length entre 18 i 16:

15: 18 18 17 18 18 17 18 18 17 18 18 17 17 17 16

Tingues en compte que alineaments diferents poden tenir les mateixes subseqüències alineades entre si.

En particular, això pot passar si els alineaments només difereixen entre si quant a la col·locació dels forats, tal com pots veure en aquest exemple:

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 7

#PairwiseAligner

L’objecte PairwiseAligner emmagatzema tots els paràmetres d’alineament que s’utilitzaran per als alineaments per parelles.

Per veure els valors de tots els paràmetres pots fer un print de l’objecte:

from Bio import Align
aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
print(aligner)

Pots veure que els valors de match_score i mode corresponen als que s’han passat com a arguments al constructor de l’objecte:

Pairwise sequence aligner with parameters
  wildcard: None
  match_score: 1.000000
  mismatch_score: 0.000000
  mode: local
  target_internal_open_gap_score: 0.000000
  ...

L’atribut mode pot tenir el valor “global” o “local”, i especifica un alineament per parelles global o local, respectivament.

En funció dels paràmetres de puntuació de “gap” i del mode, l’objecte PairwiseAligner tria automàticament l’algorisme adequat per computar l’alineament corresponent.

L’atribut algorithm només té accés de lectura i no es pot modificar.

aligner = Align.PairwiseAligner(match_score=1.0, mode="local")
assert aligner.algorithm == "Smith-Waterman"

aligner.mode = "global"
assert aligner.algorithm == "Needleman-Wunsch"

L’objecte PairwiseAligner també emmagatzema l’atribut de precisió epsilon que s’utilitza durant l’alineament.

El valor predeterminat és igual a 10⁻⁶:

aligner = Align.PairwiseAligner(match_score=1.0, mode="local")

assert aligner.epsilon == 1e-6

Dues puntuacions es consideraran iguals a efectes de l’alineament si la diferència absoluta entre elles és inferior al valor de l’atribut epsilon.

#Puntuació de substitució

Les puntuacions de substitució defineixen el valor que s’ha d’afegir a la puntuació total quan dues lletres (nucleòtids o aminoàcids) estan alineades entre si.

Els alineaments de seqüències de nucleòtids es basen normalment en puntuacions de coincidència (match) i de desajust (mismatch). Per exemple, de manera predeterminada, Blast utilitza una puntuació de coincidència de +1 i una puntuació de desajust de -2 per als alineaments de nucleòtids amb megablast, amb una penalització de gap de 2,5.

Les puntuacions de coincidència i de desajust es poden especificar establint els atributs corresponents de l’objecte 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.0

#Gaps

#gap scores predefinits

De manera predeterminada, un objecte PairwiseAligner s’inicia amb una puntuació de match de +1,0, una puntuació de mismatch de 0,0 i totes les puntuacions de gap iguals a 0,0.

Encara que això té l’avantatge de ser un esquema de puntuació senzill, en general no dona el millor rendiment. En lloc d’això, pots utilitzar l’argument de puntuació per seleccionar un esquema de puntuació predefinit en inicialitzar un objecte PairwiseAligner.

Els esquemes de puntuació que es proporcionen són:

• blastn i megablast, que són adequats per a alineaments de nucleòtids.
• blastp, que és adequat per a alineaments de proteïnes.

La selecció d’aquests esquemes de puntuació inicialitzarà l’objecte PairwiseAligner amb els paràmetres de puntuació predeterminats que fan servir BLASTN, MegaBLAST i BLASTP, respectivament.

from Bio import Align

aligner = Align.PairwiseAligner(scoring="blastn")

#Affine gap scores

Les puntuacions del gap afí es defineixen amb una puntuació per obrir un gap i una puntuació per ampliar un gap existent:

gap_score = open_gap_score + (n−1) × extend_gap_score,

on n és la longitud del gap.

Biopython permet un control precís de l’esquema de puntuació de forats mitjançant dotze atributs de l’objecte PairwiseAligner per parelles, tal com pots veure en aquest exemple:

TODO revisar

#gap function

Per a un control encara més precís de les puntuacions de gap, pots especificar una funció de puntuació de gap.

Per exemple, aquesta funció no permet un forat després de dos nucleòtids a la seqüència 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])

Pots veure que la funció no permet introduir un gap després d’AAA com es fa en el primer alineament:

target        	0 AAACTT 6
              	0 |||-|| 6
query         	0 AAA-TT 5

target        	0 AACGGTT 7
              	0 -|--.|| 7
query         	0 -A--ATT 4

#Obtenció d’informació sobre l’alineament

#Google Docs

El document continua a Google Docs: Alineament de seqüències

#TODO

• Important!. Falta completar aspectes com extended gap, odds-log ratio, matrius blosum,etc.
• Scoring protein sequence alignments using deep learning