Las matrices de sustitución más utilizadas son las PAM y las BLOSUM. Éstas matrices son las que utiliza el famoso programa de alineamiento BLAST.
En el siguiente ejemplo se mostrará una sencilla subrutina en Perl que evalúa un alineamiento de proteínas utilizando la matriz BLOSUM62. Las dos secuencias alineadas se pasan como sendas cadenas de caracteres:
1: # prog3.pl Evalúa un alineamiento
2: sub score_alignment {
3: my ($type,$score_matrix,$seq1,$seq2)=@_;
4: my (@indices,@matrix,@score_results);
5: if ($score_matrix eq 'blosum62'){
6: @indices = ('C', 'S', 'T', 'P', 'A', 'G', 'N', 'D', 'E', 'Q', 'H', 'R', 'K', 'M', 'I', 'L', 'V', 'F', 'Y', 'W');
7: @matrix = (
8: [9, -1, -1, -3, 0, -3, -3, -3, -4, -3, -3, -3, -3, -1, -1, -1, -1, -2, -2, -2],
9: [-1, 4, 1, -1, 1, 0, 1, 0, 0, 0, -1, -1, 0, -1, -2, -2, -2, -2, -2, -3],
10: [-1, 1, 4, 1, -1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -2, -2, -2, -2, -2, -3],
11: [-3, -1, 1, 7, -1, -2, -1, -1, -1, -1, -2, -2, -1, -2, -3, -3, -2, -4, -3, -4],
12: [0, 1, -1, -1, 4, 0, -1, -2, -1, -1, -2, -1, -1, -1, -1, -1, -2, -2, -2, -3],
13: [-3, 0, 1, -2, 0, 6, -2, -1, -2, -2, -2, -2, -2, -3, -4, -4, 0, -3, -3, -2],
14: [-3, 1, 0, -2, -2, 0, 6, 1, 0, 0, -1, 0, 0, -2, -3, -3, -3, -3, -2, -4],
15: [-3, 0, 1, -1, -2, -1, 1, 6, 2, 0, -1, -2, -1, -3, -3, -4, -3, -3, -3, -4],
16: [-4, 0, 0, -1, -1, -2, 0, 2, 5, 2, 0, 0, 1, -2, -3, -3, -3, -3, -2, -3],
17: [-3, 0, 0, -1, -1, -2, 0, 0, 2, 5, 0, 1, 1, 0, -3, -2, -2, -3, -1, -2],
18: [-3, -1, 0, -2, -2, -2, 1, 1, 0, 0, 8, 0, -1, -2, -3, -3, -2, -1, 2, -2],
19: [-3, -1, -1, -2, -1, -2, 0, -2, 0, 1, 0, 5, 2, -1, -3, -2, -3, -3, -2, -3],
20: [-3, 0, 0, -1, -1, -2, 0, -1, 1, 1, -1, 2, 5, -1, -3, -2, -3, -3, -2, -3],
21: [-1, -1, -1, -2, -1, -3, -2, -3, -2, 0, -2, -1, -1, 5, 1, 2, -2, 0, -1, -1],
22: [-1, -2, -2, -3, -1, -4, -3, -3, -3, -3, -3, -3, -3, 1, 4, 2, 1, 0, -1, -3],
23: [-1, -2, -2, -3, -1, -4, -3, -4, -3, -2, -3, -2, -2, 2, 2, 4, 3, 0, -1, -2],
24: [-1, -2, -2, -2, 0, -3, -3, -3, -2, -2, -3, -3, -2, 1, 3, 1, 4, -1, -1, -3],
25: [-2, -2, -2, -4, -2, -3, -3, -3, -3, -3, -1, -3, -3, 0, 0, 0, -1, 6, 3, 1],
26: [-2, -2, -2, -3, -2, -3, -2, -3, -2, -1, 2, -2, -2, -1, -1, -1, -1, 3, 7, 2],
27: [-2, -3, -3, -4, -3, -2, -4, -4, -3, -2, -2, -3, -3, -1, -3, -2, -3, 1, 2, 11]
28: );
29: }
30: if (length($seq1)==length($seq2)){
31: my @seq1 = split(//,$seq1);
32: my @seq2 = split(//,$seq2);
33: for (my $i=0; $i<=$#seq1; $i++) {
34: # Find letters position in matrix
35: my ($pos1) = grep { $indices[$_] eq $seq1[$i] } 0 .. $#indices;
36: my ($pos2) = grep { $indices[$_] eq $seq2[$i] } 0 .. $#indices;
37: if (defined($pos1) && defined($pos2)){
38: push (@score_results, $matrix[$pos1][$pos2]);
39: } else {
40: push (@score_results, '-');
41: }
42: }
43: } else {
44: die "# Both sequences must have the same length to score the alignment\n";
45: }
46: return @score_results;
47: }
Esto se podría mejorar un poco cambiando la forma de trabajar de los índices, pasando de un array a un hash:
ResponderEliminarmy %indices = (
'C' => 0, 'S' => 1, 'T' => 2, 'P' => 3, 'A' => 4,
'G' => 5, 'N' => 6, 'D' => 7, 'E' => 8, 'Q' => 9,
'H' => 10, 'R' => 11, 'K' => 12, 'M' => 13, 'I' => 14,
'L' => 15, 'V' => 16, 'F' => 17, 'Y' => 18, 'W' => 19,
);
y luego, podemos sustituir la búsqueda con grep(), a una asignación directa:
my $pos1 = $indices{ $seq1[$i] };
Interesante blog! Tengo una pregunta sobre Blosum en Hash. Mi idea era hacer un hash de hashes con la matriz de blosum, de manera que tuviera una pinta así:
ResponderEliminar%HoH = (A => {A => '1', C => '3', N => '2' etc hasta 20 },
C => {A => '1', C => '5', N => '4' hasta 20} );
y hasta 20.
Es decir 20x20 = 400 valores que dificilmente se habrán de introducir a mano. Mi pregunta es, partiendo de un file así:
Es decir, con la matriz de blosum,
- A C T R N E ...hasta 20
A 2 5 8 1 2 4 ...
C 5 2 3 1 5 2 ...
T 8 3 2 2 1 0 ..
R 1 1 ...
N 2
E ...
...
hasta 20
Como puedo hacer que el script de perl lea ésto y construya la Hash de hashes?
Muchisimas gracias!
Se me ocurre primero dividir en filas la matriz:
ResponderEliminar@rows = split("\n",$matrix)
Después leer la primera fila y extraer las posiciones de los amino ácidos:
@aas = split("\s",shift(@rows))
y luego recorrer el resto de filas desde la columna 2 anotándolos en el Hash:
forearch $row (@rows) {
@cols = split("\s",$row)
$aa1 = $cols[0]
for ($i=1;$i<=$#cols;$i++) {
$aa2 = $aas[$i]
$hash_blosum->{$aa1}{$aa2} = $cols[$i]
}
}
He omitido alguna cosilla como los fines de línea y los 'my' imagino que sabes algo de perl...
Gracias por tu sugerencia, Alvaro.
ResponderEliminarTengo un problema un poco de perogrullo, seguramente. El archivo con la matriz tiene tantas lineas como filas tiene la matriz (comprobado con el vi: "blosum62.txt" 25L, 1698C). Pero esto:
open (BLOSUM, $blosum) || die "Warning: Can't open $blosum: $!\n";
while ()
{
chomp;
@rows=split("\n",$_);
}
print "$rows[0]\n";
print "@rows\n";
Los dos prints me dan el mismo resultado: la ultima fila de la matriz (* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 1
). He probado "$line=" sustituyendo $_ por $line, también "\n" por /\n/, he probado con otro archivo de 3 lineas, y también, he probado sin chomp. No sé. Es tan sencillo que debe de ser una tontería, pero no la sé cual.
uy, dentro del while había los "<"BLOSUM">" pero el comentario los ha omitido. Idem en $line= "<"BLOSUM">".
ResponderEliminar¡Hola a todos!
ResponderEliminarEstoy iniciándome en esto del Perl, y estoy haciendo una pequeña aplicación para evaluar alinamientos de secuencias de aminoácidos.
Estoy retocando el código de ejemplo que habéis puesto para que se adapte a lo que mi aplicación ha de hacer.
¿Alguien podría explicarme esta sentencia?
my ($pos1) = grep { $indices[$_] eq $seq1[$i] } 0 .. $#indices;
Comprendo el objetivo de la instrucción pero no como está construida. ¿Que significa el 0?¿Por que se marca la variable indices como $#indices?
Gracias de antemano por la aclaración, y un saludo .
La sentencia por la que preguntas simplemente sirve para encontrar la posición de la letra del aminoácido que se va a comparar ($seq1[$i]) en el array @indices. Y una vez conocida la posición de los 2 AA a comparar ($pos1 y $pos2) se busca el valor en la matriz BLOSUM: $matrix[$pos1][$pos2]
ResponderEliminarAunque creo que lo que preguntas es la construcción de la sentencia:
'grep' lo que hace es recorrer todos los elementos de un array con una sola línea de código. En este caso el array es 0 .. $#indices (con .. se generan sucesiones de números: 0,1,2,3,4... hasta el tamaño del array indices obtenido con $#). Los valores del array se pasan a la variable $_ que es la que se evalúa hasta obtener TRUE, entonces se para 'grep' y devuelve el valor del array donde se ha detenido.
Si estás empezando con perl puedes sustituir esta lína por un bucle 'for' o 'foreach', mucho más claro:
my $pos1;
for (my $k=0; $k<=$#indices; $k++) {
if ($indices[$k] eq $seq1[i]) {
$pos1 = $k;
last;
}
}
Más info: http://perldoc.perl.org/functions/grep.html
¡Muchas gracias por la aclaración Álvaro!
ResponderEliminarHola a todos,
ResponderEliminarhe modificado el código para utilizarlo en una aplicación que compara alineamientos múltiples de secuencia. Funciona bien, el problema es que si tiene que evaluar un MSA muy grande se ralentiza mucho el proceso. He pensado que podría paralelizarse la subrutina para poder realizar de una forma más eficiente los cálculos, pero no consigo dar con una solución satisfactoria...
¿A alguien se le ocurre como podría paralelizarla para mejorar la velocidad de ejecución? Gracias y un saludo.
Tienes varios posts de paralelización en el blog, pero el problema del alineamiento múltiple es complejo y con muchas secuencias no tiene solución exacta en un tiempo de cálculo realista. Además el código aquí publicado es muy sencillo, deberías usar algoritmos de programación dinámica para optimizar el tiempo de cálculo, te recomiendo que leas los artículos de la Wikipedia sobre alineamientos y el libro de J. Pevsner "Bioinformatics and Functional Genomics". Mi consejo es que uses programas que han hecho otros como Clustal Omega, MUSCLE o T-Coffee, funcionan muy bien en línea de comandos, no te rompas la cabeza en algo que ya está muy optimizado.
ResponderEliminarHola. En primer lugar, gracias a los autores del blog por sus interesantes apuntes sobre bioinformática. En segundo lugar, quisiera comentarle que soy principiante en este campo y mi consulta es: hice un alineamiento de 22 proteinas en T-Coffee, hay varios bloques separados por gaps y dos bloques grandes con muchos aminoácidos conservados. Mis preguntas son: ¿qué textos me recomendaría para aprender a interpretar un MSA y saber si es correcto o no?. Tengo entendido que para elegir una matriz PAM o BLOSUM (que son las que muestran en T-Coffee) es necesario conocer el % de divergencia entre las proteínas, ¿cómo haría yo para saber este valor?. Gracias
ResponderEliminarGracias a ti por leernos ;) Tenemos otra entrada del blog que habla sobre las matrices de sustitución (http://bioinfoperl.blogspot.com.es/2012/07/matrices-de-sustitucion-y-alineamiento.html), también te recomendaría el libro 'Bioinformatics For Dummies' que trata 95% sobre alineamiento de secuencias y está escrito por el creador de T-Coffee y para un conocimiento más avanzado "Bioinformatics and Functional Genomics" de J. Pevsner. El tema de los alineamientos múltiples es complejo y normalmente la diferencia entre usar una u otra matriz es poca (salvo casos extremos).
ResponderEliminarMuchas gracias por sus recomendaciones, estaré leyendo esos libros. Sólo otra pregunta más: ¿por qué dice que el tema de los alineamientos múltiples es complejo?. Gracias :)
Eliminar