Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Alineamiento con restricciones

Perl aplicado a la bioinformática
Publicar una respuesta

Alineamiento con restricciones

Notapor sergio4393 » 2017-06-13 07:36 @358

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.006 segundos, usando GeSHi 1.0.8.4
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba

Publicidad

Re: Alineamiento con restricciones

Notapor explorer » 2017-06-13 15:42 @696

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.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Alineamiento con restricciones

Notapor sergio4393 » 2017-06-14 04:08 @214

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.
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba

Re: Alineamiento con restricciones

Notapor explorer » 2017-06-14 06:55 @329

Entonces, el poner una restricción inicial, ¿consiste en recorrer la matriz de alineamiento desde una determinada celda?
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Alineamiento con restricciones

Notapor sergio4393 » 2017-06-14 07:28 @352

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.
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba

Re: Alineamiento con restricciones

Notapor sergio4393 » 2017-06-14 08:59 @416

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.004 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
Última edición por explorer el 2017-06-14 18:03 @794, editado 1 vez en total
Razón: Formateado de código con perltidy
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba

Re: Alineamiento con restricciones

Notapor explorer » 2017-06-15 13:08 @588

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?
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Alineamiento con restricciones

Notapor sergio4393 » 2017-06-16 07:43 @363

Sí, tienes razón, funciona a la perfección, pero, no sé, la forma es un poco cutre :P
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba

Re: Alineamiento con restricciones

Notapor explorer » 2017-06-16 11:20 @514

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.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Alineamiento con restricciones

Notapor sergio4393 » 2017-06-19 04:57 @248

Ya no sé...

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

¡Gracias por todo!

Un saludo.
sergio4393
Perlero nuevo
Perlero nuevo
 
Mensajes: 6
Registrado: 2017-06-09 07:46 @365
Arriba
Publicar una respuesta

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 0 invitados

Powered by phpBB® Forum Software © phpBB Group
Traducción al Español Mexicano por nextgen en colaboración con phpbb-es.com
phpBB SEO
Advertisements by Advertisement Management