4 de junio de 2010

Evaluación de un alineamiento de proteínas

Un alineamiento de secuencias proteicas se puede evaluar numéricamente comparando las identidades entre los aminoácidos alineados. Se puede asignar una puntuación a cada par de aminoácidos usando una matriz de sustitución, que describe el ritmo al que un aminoácido en una secuencia cambia a otro a lo largo de la evolución natural.

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

13 comentarios:

  1. Esto se podría mejorar un poco cambiando la forma de trabajar de los índices, pasando de un array a un hash:

    my %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] };

    ResponderEliminar
  2. 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í:

    %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!

    ResponderEliminar
  3. Se me ocurre primero dividir en filas la matriz:
    @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...

    ResponderEliminar
  4. Gracias por tu sugerencia, Alvaro.
    Tengo 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.

    ResponderEliminar
  5. uy, dentro del while había los "<"BLOSUM">" pero el comentario los ha omitido. Idem en $line= "<"BLOSUM">".

    ResponderEliminar
  6. ¡Hola a todos!
    Estoy 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 .

    ResponderEliminar
  7. 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]

    Aunque 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

    ResponderEliminar
  8. ¡Muchas gracias por la aclaración Álvaro!

    ResponderEliminar
  9. Hola a todos,
    he 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.

    ResponderEliminar
  10. 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.

    ResponderEliminar
  11. Hola. 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

    ResponderEliminar
  12. Gracias 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).

    ResponderEliminar
    Respuestas
    1. Muchas 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