Página 1 de 1

Alineamiento con restricciones

NotaPublicado: 2017-06-13 07:36 @358
por sergio4393
Hola. El script de abajo es el alineamiento global para dos secuencias. Ahora lo que tengo que realizar es el alineamiento con restricciones.

Por ejemplo, quiero alinear la A de la primera secuencia (posición 0) con la G de la segunda secuencia (posición 1), es decir, si tengo como restricción (0,1) quiere decir que la posición 0 de la secuencia 1 tiene que alinear con la posición 1 de la secuencia 2, independientemente de si hay coincidencia o no.

Sé que para que se lleve a cabo, desde el punto de vista de la matriz de programación dinámica, la reconstrucción tiene que pasar por (0,1) y en esa celda retroceder en diagonal.

He estado probando con if() marcando la posición cuando se recorren todos los índices pero no consigo llegar a nada. Si alguien me da alguna pista o algo para poder realizarlo se lo agradecería.

Un saludo y gracias de antemano.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. use strict;
  2.        
  3. my $MATCH=2;
  4. my $MISMATCH=-1;
  5. my $GAP=-1;
  6. my $seq1='ATGCTTA';
  7. my $seq2='TGCATTAA';
  8.        
  9. sub alineamiento {
  10.         my($seq1,$seq2)=@_;
  11.        
  12.         my @matrix;
  13.         $matrix[0][0]{score}   = 0;
  14.         $matrix[0][0]{pointer} = "none";
  15.        
  16.         for(my $j = 1; $j <= length($seq1); $j++) {
  17.             $matrix[0][$j]{score}   = $GAP * $j;
  18.             $matrix[0][$j]{pointer} = "left";
  19.         }
  20.         for (my $i = 1; $i <= length($seq2); $i++) {
  21.             $matrix[$i][0]{score}   = $GAP * $i;
  22.             $matrix[$i][0]{pointer} = "up";
  23.         }
  24.        
  25.         # fill
  26.         for(my $i = 1; $i <= length($seq2); $i++) {
  27.             for(my $j = 1; $j <= length($seq1); $j++) {
  28.                 my ($diagonal_score, $left_score, $up_score);
  29.        
  30.                 # calculate match score
  31.                 my $letter1 = substr($seq1, $j-1, 1);
  32.                 my $letter2 = substr($seq2, $i-1, 1);                            
  33.                 if ($letter1 eq $letter2) {
  34.                     $diagonal_score = $matrix[$i-1][$j-1]{score} + $MATCH;
  35.                 }
  36.                 else {
  37.                     $diagonal_score = $matrix[$i-1][$j-1]{score} + $MISMATCH;
  38.                 }
  39.        
  40.                 # calculate gap scores
  41.                 $up_score   = $matrix[$i-1][$j]{score} + $GAP;
  42.                 $left_score = $matrix[$i][$j-1]{score} + $GAP;
  43.        
  44.                 # choose best score
  45.                 if ($diagonal_score >= $up_score) {
  46.                     if ($diagonal_score >= $left_score) {
  47.                         $matrix[$i][$j]{score}   = $diagonal_score;
  48.                         $matrix[$i][$j]{pointer} = "diagonal";
  49.                     }
  50.                 else {
  51.                         $matrix[$i][$j]{score}   = $left_score;
  52.                         $matrix[$i][$j]{pointer} = "left";
  53.                     }
  54.                 } else {
  55.                     if ($up_score >= $left_score) {
  56.                         $matrix[$i][$j]{score}   = $up_score;
  57.                         $matrix[$i][$j]{pointer} = "up";
  58.                     }
  59.                     else {
  60.                         $matrix[$i][$j]{score}   = $left_score;
  61.                         $matrix[$i][$j]{pointer} = "left";
  62.                     }
  63.                 }
  64.             }
  65.         }
  66.        
  67.        
  68.         my $align1 = "";
  69.         my $align2 = "";
  70.        
  71.         # start at last cell of matrix
  72.         my $j = length($seq1);
  73.         my $i = length($seq2);
  74.        
  75.         while (1) {
  76.             last if $matrix[$i][$j]{pointer} eq "none"; # ends at first cell of matrix
  77.        
  78.             if ($matrix[$i][$j]{pointer} eq "diagonal") {
  79.                 $align1 .= substr($seq1, $j-1, 1);
  80.                 $align2 .= substr($seq2, $i-1, 1);
  81.                 $i--;
  82.                 $j--;
  83.             }
  84.             elsif ($matrix[$i][$j]{pointer} eq "left") {
  85.                 $align1 .= substr($seq1, $j-1, 1);
  86.                 $align2 .= "-";
  87.                 $j--;
  88.             }
  89.             elsif ($matrix[$i][$j]{pointer} eq "up") {
  90.                 $align1 .= "-";
  91.                 $align2 .= substr($seq2, $i-1, 1);
  92.                 $i--;
  93.             }    
  94.         }
  95.         $align1 = reverse $align1;
  96.         $align2 = reverse $align2;
  97.         return $matrix[$i-1][$j-1]{score},$align1,$align2
  98. }
  99.  
  100. my @resultado=alineamiento($seq1,$seq2);
  101.  
  102. print $resultado[1],"\n";
  103. print $resultado[2],"\n";
  104. print $resultado[0];
  105.  
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-13 15:42 @696
por explorer
Bienvenido a los foros de Perl en Español, sergio4393.

No me queda claro cuál es el resultado que quieres conseguir. ¿Algo como esto?
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
-ATGCTTA-
TG-CATTAA
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Pero, en ese caso, hay pocos alineamientos, ¿no?

Hay algunos hilos en este foro que hablan de los alineamientos. En el primero de ellos pongo un procedimiento para pintar la tabla de alineamientos.

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-14 04:08 @214
por sergio4393
Sì, exacto, eso es lo que debería devolver.

Forzar el alineamiento de una posición de la secuencia 1 con otra posición de la secuencia 2, independientemente de si hay coincidencia o no.

En este ejemplo no hay muchas posibilidades, la verdad, pero lo suyo es que poniendo secuencias más largas y con más restricciones funcione bien.

Voy a echarle un vistazo a esos enlaces.

Un saludo.

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-14 06:55 @329
por explorer
Entonces, el poner una restricción inicial, ¿consiste en recorrer la matriz de alineamiento desde una determinada celda?

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-14 07:28 @352
por sergio4393
No desde una determinada celda, si no que tu recorres la matriz de manera normal, y cuando llegue a esas posiciones tiene que alinear en diagonal. Tiene que alinearse ya sea igual la letra (match) o distinta (mismatch), pero no con gap.

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-14 08:59 @416
por sergio4393
He estado probando y por ejemplo una manera de forzar el alineamiento es dándole un valor alto para esa posición pero, claro, luego la puntuación no da bien, como es lógico.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. use strict;
  2.  
  3. my $MATCH    = 2;
  4. my $MISMATCH = -1;
  5. my $GAP      = -1;
  6. my $seq1     = 'ATGCTTA';
  7. my $seq2     = 'TGCATTAA';
  8.  
  9. sub alineamiento {
  10.     my ( $seq1, $seq2 ) = @_;
  11.  
  12.     my @matrix;
  13.     $matrix[0][0]{score}   = 0;
  14.     $matrix[0][0]{pointer} = "none";
  15.  
  16.     for ( my $j = 1; $j <= length($seq1); $j++ ) {
  17.         $matrix[0][$j]{score}   = $GAP * $j;
  18.         $matrix[0][$j]{pointer} = "left";
  19.     }
  20.     for ( my $i = 1; $i <= length($seq2); $i++ ) {
  21.         $matrix[$i][0]{score}   = $GAP * $i;
  22.         $matrix[$i][0]{pointer} = "up";
  23.     }
  24.  
  25.     # fill
  26.     for ( my $i = 1; $i <= length($seq2); $i++ ) {
  27.         for ( my $j = 1; $j <= length($seq1); $j++ ) {
  28.             my ( $diagonal_score, $left_score, $up_score );
  29.             if ( $i == 1 and $j == 1 ) {
  30.                 $diagonal_score = $matrix[ $i - 1 ][ $j - 1 ]{score} + 99;
  31.  
  32.             }
  33.             else {
  34.                 # calculate match score
  35.                 my $letter1 = substr( $seq1, $j - 1, 1 );
  36.                 my $letter2 = substr( $seq2, $i - 1, 1 );
  37.                 if ( $letter1 eq $letter2 ) {
  38.                     $diagonal_score = $matrix[ $i - 1 ][ $j - 1 ]{score} + $MATCH;
  39.                 }
  40.                 else {
  41.                     $diagonal_score = $matrix[ $i - 1 ][ $j - 1 ]{score} + $MISMATCH;
  42.                 }
  43.             }
  44.  
  45.             # calculate gap scores
  46.             $up_score   = $matrix[ $i - 1 ][$j]{score} + $GAP;
  47.             $left_score = $matrix[$i][ $j - 1 ]{score} + $GAP;
  48.  
  49.             # choose best score
  50.             if ( $diagonal_score >= $up_score ) {
  51.                 if ( $diagonal_score >= $left_score ) {
  52.                     $matrix[$i][$j]{score}   = $diagonal_score;
  53.                     $matrix[$i][$j]{pointer} = "diagonal";
  54.                 }
  55.                 else {
  56.                     $matrix[$i][$j]{score}   = $left_score;
  57.                     $matrix[$i][$j]{pointer} = "left";
  58.                 }
  59.             }
  60.             else {
  61.                 if ( $up_score >= $left_score ) {
  62.                     $matrix[$i][$j]{score}   = $up_score;
  63.                     $matrix[$i][$j]{pointer} = "up";
  64.                 }
  65.                 else {
  66.                     $matrix[$i][$j]{score}   = $left_score;
  67.                     $matrix[$i][$j]{pointer} = "left";
  68.                 }
  69.             }
  70.         }
  71.     }
  72.  
  73.     my $align1 = "";
  74.     my $align2 = "";
  75.  
  76.     # start at last cell of matrix
  77.     my $j = length($seq1);
  78.     my $i = length($seq2);
  79.  
  80.     while (1) {
  81.         last if $matrix[$i][$j]{pointer} eq "none";    # ends at first cell of matrix
  82.  
  83.         if ( $matrix[$i][$j]{pointer} eq "diagonal" ) {
  84.             $align1 .= substr( $seq1, $j - 1, 1 );
  85.             $align2 .= substr( $seq2, $i - 1, 1 );
  86.             $i--;
  87.             $j--;
  88.         }
  89.         elsif ( $matrix[$i][$j]{pointer} eq "left" ) {
  90.             $align1 .= substr( $seq1, $j - 1, 1 );
  91.             $align2 .= "-";
  92.             $j--;
  93.         }
  94.         elsif ( $matrix[$i][$j]{pointer} eq "up" ) {
  95.             $align1 .= "-";
  96.             $align2 .= substr( $seq2, $i - 1, 1 );
  97.             $i--;
  98.         }
  99.     }
  100.     $align1 = reverse $align1;
  101.     $align2 = reverse $align2;
  102.     return $matrix[ $i - 1 ][ $j - 1 ]{score}, $align1, $align2;
  103. }
  104.  
  105. my @resultado = alineamiento( $seq1, $seq2 );
  106.  
  107. print $resultado[1], "\n";
  108. print $resultado[2], "\n";
  109. print $resultado[0];
  110.  
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

Resultado:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
ATGC-TT-A
T-GCATTAA
106
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-15 13:08 @588
por explorer
Pero, si le sumas 99 en esa casilla, ¿no vale con restar 99 (o 98 o 100, según el contenido original de la celda) al resultado final?

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-16 07:43 @363
por sergio4393
Sí, tienes razón, funciona a la perfección, pero, no sé, la forma es un poco cutre :P

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-16 11:20 @514
por explorer
Humm... vamos a ver...

La resolución del problema del alineamiento se hace en dos fases: en la primera rellenamos la matriz, viendo las coincidencias, y en la segunda, recorremos la matriz en orden inverso.

El problema que yo veo es que si hacemos solo esto, no está garantizado que el recorrido pase por la celda que nos interesa, ¿verdad?

Entonces... hago un par de suposiciones. Desde la celda hasta el final del recorrido, sí que tenemos claro qué tenemos que hacer (pues partimos de una celda en concreto). El problema está en cómo conectar el final de la matriz de alineamiento, y comienzo del recorrido, con la celda que nos interesa.

Yo -no lo he comprobado- estoy suponiendo que la solución sería partir de la celda en sentido contrario al recorrido normal, y con ello llegar a la celda de partida (la que sería la primera celda de partida). No tengo claro, ahora, si eso se puede hacer, ya que el recorrido se va formando en la primera fase.

Re: Alineamiento con restricciones

NotaPublicado: 2017-06-19 04:57 @248
por sergio4393
Ya no sé...

Yo creo que lo dejaré así de la manera mencionada arriba y listo.

¡Gracias por todo!

Un saludo.