• Publicidad

Puntuación errónea en alineamiento

Perl aplicado a la bioinformática

Puntuación errónea en alineamiento

Notapor perl junior » 2017-01-26 05:00 @250

Buenas, he hecho un programa para calcular la puntuación en un alineamiento, y no entiendo por qué me devuelve el valor 6 y no 4 como debería ser.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. use List::Util qw[min max];
  2. sub programacion_dinamica{
  3.         my($secuencia1,$secuencia2) =@_;
  4.        
  5.         my $len1 =length ($secuencia1);
  6.         my $len2 =length($secuencia2);
  7.         my @secuencia1 = split("",$secuencia1);
  8.         my @secuencia2= split("",$secuencia2);
  9.        
  10.        
  11.         #Construimos la tabla s y la inicializamos a 0
  12.         my @s=();
  13.         my @decisiones=();
  14.         for (my $i =0; $i <=$len1 ;$i++){
  15.        
  16.                 my @ fila_s=();
  17.                 my @fila_decision=();
  18.                 for (my $j=0;$j <=$len2;$j++){
  19.                         push @fila_s,0;
  20.                         push@fila_decision,'';
  21.                         print @fila_decision;
  22.                 }
  23.                 push @s,\@fila_s; #por referencia
  24.                 push @decisiones,\@fila_decision;
  25.                
  26.         }
  27.  
  28.        
  29.        
  30.         for (my $i =0;$i<$len1;$i++){
  31.                 $s[0][0]= 0;
  32.                 $s[$i][0] =$s[$i][0];
  33.                
  34.                 $decisiones[$i][0]="abajo";
  35.         }
  36.  
  37.         for (my $j=0;$j<=$len2;$j++){
  38.                
  39.                 $s[0][$j] =$s[0][$j];
  40.                
  41.                 $decisiones[0][$j]="derecha";
  42.         }
  43.  
  44.        
  45.         for (my $j=1;$j<$len2;$j++){
  46.                 for (my $i=1;$i<$len1;$i++){
  47.                        
  48.                        
  49.                        
  50.                         if ($secuencia1[$i-1]== $secuencia2[$j-1]){
  51.                        
  52.                                 $suma=1;
  53.                         }
  54.                         else{
  55.                                 $suma=0;}
  56.                                
  57.                         $s[$i][$j]=max(
  58.                                 $s[$i-1][$j]+0,
  59.                                 $s[$i][$j-1]+0,
  60.                                 $s[$i-1][$j-1]+$suma);
  61.                        
  62.                        
  63.                     if($s[$i][$j]== $s[$i][$j-1]+0){
  64.                                 $decisiones[$i][$j] ="abajo";
  65.                         }
  66.                        
  67.                         elsif ($s[$i][$j]==$s[$i-1][$j-1]+$suma){
  68.                                 $decisiones[$i][$j] ="diagonal";
  69.                         }
  70.                         else{
  71.                                 $decisiones[$i][$j]="derecha";
  72.                     }
  73.                
  74.                 }
  75.         }
  76.        
  77.        
  78.         return$s[$len1-1][$len2-1];
  79.        
  80. }
  81.  
  82.  
  83. #Construimos las tablas y decisiones
  84.  
  85.  
  86.  
  87. $secuencia1="ATGTTATA";
  88. $secuencia2="ATCGTCC";
  89.  
  90.  
  91. my $resultado=programacion_dinamica($secuencia1,$secuencia2);
  92. print $resultado;
  93.  
Coloreado en 0.022 segundos, usando GeSHi 1.0.8.4

Muchas gracias a cualquier pista que me ayude a ver el fallo.
perl junior
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Publicidad

Re: Puntuación errónea en alineamiento

Notapor explorer » 2017-01-26 07:10 @340

En la línea 50, debes cambiar la comparación numérica '==' por la alfanumérica 'eq', ya que lo que estás comparando son (las letras de) las bases.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/env perl
  2. use v5.14;
  3.  
  4. use List::Util qw[max];
  5.  
  6. sub programacion_dinamica {
  7.     my($secuencia1, $secuencia2) = @_;
  8.  
  9.     my @secuencia1 = split //, $secuencia1;
  10.     my @secuencia2 = split //, $secuencia2;
  11.  
  12.     my @s;
  13.     my @decisiones;
  14.  
  15.     print "  ", join("  " => @secuencia1), "\n";
  16.  
  17.     for (my $j = 0; $j < @secuencia2; $j++) {
  18.         print "$secuencia2[$j] ";
  19.  
  20.         for (my $i = 0; $i < @secuencia1; $i++) {
  21.  
  22.             my $iguales = $secuencia1[$i] eq $secuencia2[$j];               # 1 si hay coincidencia
  23.  
  24.             $s[$i+1][$j+1] = max(
  25.                 0        + $s[$i]  [$j+1],
  26.                 0        + $s[$i+1][$j],
  27.                 $iguales + $s[$i]  [$j]
  28.             );
  29.  
  30.             if ($s[$i+1][$j+1] == $s[$i+1][$j]) {
  31.                 $decisiones[$i+1][$j+1] = "v";
  32.             }
  33.  
  34.             elsif ($s[$i+1][$j+1] == $s[$i][$j]+$iguales) {
  35.                 $decisiones[$i+1][$j+1] = "\\";
  36.             }
  37.             else {
  38.                 $decisiones[$i+1][$j+1] = ">";
  39.             }
  40.  
  41.             print "$s[$i+1][$j+1]$decisiones[$i+1][$j+1] ";
  42.         }
  43.         print "\n";
  44.     }
  45.  
  46.     return $s[@secuencia1][@secuencia2];
  47. }
  48.  
  49.  
  50. #my $secuencia1 = "ATGTTATA";
  51. #my $secuencia2 = "ATCGTCC";
  52. my $secuencia1 = "ATGCTTA";
  53. my $secuencia2 = "TGCATTAA";
  54.  
  55. say programacion_dinamica($secuencia1, $secuencia2);
Coloreado en 0.017 segundos, usando GeSHi 1.0.8.4

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  A  T  G  T  T  A  T  A
A 1\ 1> 1> 1> 1> 1\ 1> 1\
T 1v 2\ 2> 2\ 2\ 2> 2\ 2>
C 1v 2v 2v 2v 2v 2v 2v 2v
G 1v 2v 3\ 3> 3> 3> 3> 3>
T 1v 2v 3v 4\ 4\ 4> 4\ 4>
C 1v 2v 3v 4v 4v 4v 4v 4v
C 1v 2v 3v 4v 4v 4v 4v 4v
4
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Puntuación errónea en alineamiento

Notapor perl junior » 2017-01-26 09:09 @423

Muchas gracias, no había caído. Efectivamente, si lo pruebo da el resultado deseado.

Sin embargo si pongo estas 2 secuencias:

secuencia1: ATGCTTA
secuencia2: TGCATTAA

¿El resultado que debe salir no es 6? ¿Por qué devuelve 5?
perl junior
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Re: Puntuación errónea en alineamiento

Notapor explorer » 2017-01-26 09:38 @443

Es que me parece que no haces todas las comparaciones de todas las bases. Es decir, los bucles for() van de 1 a longitud(secuencia)-1. La comparación se hace con ($i-1, $j-1), por lo que entonces estamos comparando las bases que van de 0 a longitud(secuencia)-2, lo cual nos deja fuera la última base de cada secuencia.

Entonces... en los bucles for, en lugar de poner '<', pondría '<=', para que hagan una vuelta más.

Y en el return,
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         return$s[$len1-1][$len2-1];
Coloreado en 0.006 segundos, usando GeSHi 1.0.8.4

cambiarlo para que sea
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         return$s[$len1][$len2];
Coloreado en 0.006 segundos, usando GeSHi 1.0.8.4


Mi código también fallaba en el return, así que lo he cambiado y ya sale bien.
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  A  T  G  C  T  T  A
T 0v 1\ 1> 1> 1\ 1\ 1>
G 0v 1v 2\ 2> 2> 2> 2>
C 0v 1v 2v 3\ 3> 3> 3>
A 1\ 1v 2v 3v 3v 3v 4\
T 1v 2\ 2v 3v 4\ 4\ 4v
T 1v 2v 2v 3v 4v 5\ 5>
A 1v 2v 2v 3v 4v 5v 6\
A 1v 2v 2v 3v 4v 5v 6v
6
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Puntuación errónea en alineamiento

Notapor perl junior » 2017-01-26 10:12 @466

Siento ser pesado e insistente. Lo de la puntuación lo tengo arreglado, pero no cuando quiero obtener el camino y que me devuelva las bases en las que coinciden.

Si pongo estas secuencias, me lo da perfectamente:
S1= ATGTTATA
S2= ATCGTCC

Pero si pongo estas
S1= ATGCTTA
S2= ATCGTCC

me devuelve este error: Modification of non-creatable array value attempted, subscript -9 at puntuacion.pl line 96.

Y no lo entiendo por qué. Siento las molestias y muchísimas gracias por la ayuda.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. use List::Util qw[min max];
  2. sub programacion_dinamica{
  3.         my($secuencia1,$secuencia2) =@_;
  4.        
  5.         my $len1 =length ($secuencia1);
  6.         my $len2 =length($secuencia2);
  7.         my @secuencia1 = split("",$secuencia1);
  8.         my @secuencia2= split("",$secuencia2);
  9.        
  10.        
  11.         #Construimos la tabla s y la inicializamos a 0
  12.         my @s=();
  13.         my @decisiones=();
  14.         for (my $i =0; $i <=$len1 ;$i++){
  15.        
  16.                 my @ fila_s=();
  17.                 my @fila_decision=();
  18.                 for (my $j=0;$j <=$len2;$j++){
  19.                         push @fila_s,0;
  20.                         push@fila_decision,'';
  21.                         print @fila_decision;
  22.                 }
  23.                 push @s,\@fila_s; #por referencia
  24.                 push @decisiones,\@fila_decision;
  25.                
  26.         }
  27.  
  28.        
  29.        
  30.         for (my $i =0;$i<$len1;$i++){
  31.                 $s[0][0]= 0;
  32.                 $s[$i][0] =$s[$i][0];
  33.                
  34.                 $decisiones[$i][0]="abajo";
  35.         }
  36.  
  37.         for (my $j=0;$j<=$len2;$j++){
  38.                
  39.                 $s[0][$j] =$s[0][$j];
  40.                
  41.                 $decisiones[0][$j]="derecha";
  42.         }
  43.  
  44.        
  45.         for (my $j=1;$j<= $len2;$j++){
  46.                 for (my $i=1;$i<= $len1;$i++){
  47.                        
  48.                        
  49.                        
  50.                         if ($secuencia1[$i-1]eq $secuencia2[$j-1]){
  51.                        
  52.                                 $suma=1;
  53.                         }
  54.                         else{
  55.                                 $suma=0;}
  56.                                
  57.                         $s[$i][$j]=max(
  58.                                 $s[$i-1][$j]+0,
  59.                                 $s[$i][$j-1]+0,
  60.                                 $s[$i-1][$j-1]+$suma);
  61.                        
  62.                        
  63.                     if($s[$i][$j]== $s[$i][$j-1]+0){
  64.                                 $decisiones[$i][$j] ="abajo";
  65.                         }
  66.                        
  67.                         elsif ($s[$i][$j]==$s[$i-1][$j-1]+$suma){
  68.                                 $decisiones[$i][$j] ="diagonal";
  69.                         }
  70.                         else{
  71.                                 $decisiones[$i][$j]="derecha";
  72.                     }
  73.                
  74.                 }
  75.         }
  76.        
  77.        
  78.        
  79.         @camino=();
  80.         @secuencia=();
  81.         @camino =construir_camino($len1-1,$len2-1,\@decisiones,\@secuencia1,@camino,@secuencia);
  82.        
  83.         return $s[$len1][$len2],@camino;
  84.        
  85. }
  86.  
  87.  
  88.  
  89.  
  90. sub construir_camino{
  91.         my ($i,$j,$decisiones_ref,$secuencia1_ref,@camino,@secuencia)=@_;
  92.         if($i ==0 and $j ==0){
  93.                 return (@camino);
  94.                
  95.         }
  96.         elsif (${$decisiones_ref}[$i][$j] eq "derecha"){
  97.                
  98.                 unshift @camino,"derecha";
  99.                 return construir_camino($i-1,$j,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  100.         }
  101.         elsif (${$decisiones_ref}[$i][$j] eq "abajo"){
  102.                
  103.                 unshift @camino,"abajo";
  104.                 return construir_camino($i,$j-1,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  105.         }
  106.         else{
  107.                 unshift @camino,"diagonal";
  108.                 unshift @secuencia,${$secuencia1_ref}[$i];
  109.                 return construir_camino($i-1,$j-1,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  110.         }
  111. }
  112.  
  113.  
  114.  
  115.  
  116.  
  117.  
  118.  
  119.  
  120. $secuencia1="ATGCTTA";
  121. $secuencia2="TGCATTAA";
  122.  
  123.  
  124.  
  125. my ($resultado,@secuencia)=programacion_dinamica($secuencia1,$secuencia2);
  126. print $resultado,"\n";
  127. print join(" ",@secuencia);
  128.  
  129.  
Coloreado en 0.026 segundos, usando GeSHi 1.0.8.4
perl junior
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Re: Puntuación errónea en alineamiento

Notapor explorer » 2017-01-26 11:45 @531

El problema está en la condición de parada de la creación del camino.

En la línea 92 compruebas si AMBOS índices, el $i y el $j son IGUALES a 0, y en ese caso, regresar con la solución.

Pero... hay situaciones en las cuales los valores se vuelven negativos (aún no sé por qué), por lo que es complicado que se dé la condición de igualdad a 0 en los dos índices a la vez. Cuando se produce el error que muestras, $i vale -9, y $j vale -12. En Perl, los valores negativos se refieren a los índices del array empezando por el final (-1 se refiere al índice del último elemento del array, el -2, el penúltimo, etc.), así que cuando ponemos un valor muy negativo, Perl no puede ir "más allá" del principio del array, por lo que salta el error.

Cambiando la expresión por

if($i <= 0 and $j <= 0) {

Con ese cambio, ya sale bien (aunque con la secuencia de las bases al revés. Para darle la vuelta, usa el comando reverse).
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Puntuación errónea en alineamiento

Notapor perl junior » 2017-01-26 12:05 @545

Pero la función reverse solo me valdría para ese caso, si se diera el de las anteriores secuencias no valdría ¿no?
perl junior
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Re: Puntuación errónea en alineamiento

Notapor explorer » 2017-01-26 12:08 @547

P.D. Estoy casi seguro que es un problema de índices. Quiero decir que hay que tener claro siempre a qué corresponde la primera magnitud (las filas), y la segunda magnitud (las columnas) y que sus respectivas variables de índice no cambian ($i, $j). Y conforme a eso, asignar lo de 'derecha', 'abajo'...

Cuando llegue a casa lo miro con más calma, pero es justo eso: no está garantizado que el algoritmo comienza justo en la esquina inferior derecha, y por eso nunca llega a parar en el (0,0), que sería lo lógico.
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Puntuación errónea en alineamiento

Notapor explorer » 2017-01-26 12:12 @550

Si antes lo digo...

Con "arreglar" lo de 'abajo' y 'derecha', ya funciona todo.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/env perl
  2. #
  3. use List::Util qw[min max];
  4.  
  5. sub programacion_dinamica{
  6.         my($secuencia1,$secuencia2) =@_;
  7.  
  8.         my $len1 =length ($secuencia1);
  9.         my $len2 =length($secuencia2);
  10.         my @secuencia1 = split("",$secuencia1);
  11.         my @secuencia2= split("",$secuencia2);
  12.  
  13.         #Construimos la tabla s y la inicializamos a 0
  14.         my @s=();
  15.         my @decisiones=();
  16.         for (my $i =0; $i <=$len1 ;$i++){
  17.  
  18.                 my @ fila_s=();
  19.                 my @fila_decision=();
  20.                 for (my $j=0;$j <=$len2;$j++){
  21.                         push @fila_s,0;
  22.                         push@fila_decision,'';
  23.                         print @fila_decision;
  24.                 }
  25.                 push @s,\@fila_s; #por referencia
  26.                 push @decisiones,\@fila_decision;
  27.         }
  28.  
  29.  
  30.         for (my $i =0;$i<$len1;$i++){
  31.                 $s[0][0]= 0;
  32.                 $s[$i][0] =$s[$i][0];
  33.  
  34.                 $decisiones[$i][0]="abajo";
  35.         }
  36.  
  37.         for (my $j=0;$j<=$len2;$j++){
  38.  
  39.                 $s[0][$j] =$s[0][$j];
  40.  
  41.                 $decisiones[0][$j]="derecha";
  42.         }
  43.  
  44.         for (my $j=1;$j<= $len2;$j++){
  45.                 for (my $i=1;$i<= $len1;$i++){
  46.  
  47.                         if ($secuencia1[$i-1] eq $secuencia2[$j-1]){
  48.                                 $suma=1;
  49.                         }
  50.                         else{
  51.                                 $suma=0;
  52.                         }
  53.  
  54.                         $s[$i][$j]=max(
  55.                                 $s[$i-1][$j]+0,
  56.                                 $s[$i][$j-1]+0,
  57.                                 $s[$i-1][$j-1]+$suma);
  58.  
  59.                         if($s[$i][$j]== $s[$i][$j-1]+0){
  60.                                 $decisiones[$i][$j] ="derecha";
  61.                         }
  62.                         elsif ($s[$i][$j]==$s[$i-1][$j-1]+$suma){
  63.                                 $decisiones[$i][$j] ="diagonal";
  64.                         }
  65.                         else{
  66.                                 $decisiones[$i][$j]="abajo";
  67.                     }
  68.                 }
  69.         }
  70.  
  71.         @camino=();
  72.         @secuencia=();
  73.  
  74.         use Data::Dump;
  75.         dd @s;
  76.         #exit;
  77.  
  78.  
  79.         @camino =construir_camino($len1,$len2,\@decisiones,\@secuencia1,@camino,@secuencia);
  80.  
  81.         return $s[$len1][$len2],@camino;
  82. }
  83.  
  84.  
  85.  
  86.  
  87. sub construir_camino{
  88.     my ($i,$j,$decisiones_ref,$secuencia1_ref,@camino,@secuencia)=@_;
  89.  
  90.     use Data::Dumper;
  91.     print "$i:$j: \n";# , Dumper $decisiones_ref;
  92.  
  93.     if($i <=0 and $j <=0){
  94.         return (@camino);
  95.  
  96.     }
  97.     elsif ($decisiones_ref->[$i][$j] eq "abajo"){
  98.  
  99.         unshift @camino,"derecha";
  100.         return construir_camino($i-1,$j,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  101.     }
  102.     elsif ($decisiones_ref->[$i][$j] eq "derecha"){
  103.  
  104.         unshift @camino,"abajo";
  105.         return construir_camino($i,$j-1,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  106.     }
  107.     else{
  108.         unshift @camino,"diagonal";
  109.         unshift @secuencia,${$secuencia1_ref}[$i];
  110.         return construir_camino($i-1,$j-1,$decisiones_ref,$secuencia1_ref,@camino,@secuencia);
  111.     }
  112. }
  113.  
  114.  
  115. $secuencia1="ATGCTTA";
  116. $secuencia2="TGCATTAA";
  117.  
  118.  
  119.  
  120. my ($resultado,@secuencia)=programacion_dinamica($secuencia1,$secuencia2);
  121. print $resultado,"\n";
  122. print join(" ",@secuencia);
Coloreado en 0.026 segundos, usando GeSHi 1.0.8.4


Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
(
  [0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 1, 1, 1, 1, 1],
  [0, 1, 1, 1, 1, 2, 2, 2, 2],
  [0, 1, 2, 2, 2, 2, 2, 2, 2],
  [0 .. 3, 3, 3, 3, 3, 3],
  [0 .. 3, 3, 4, 4, 4, 4],
  [0 .. 3, 3, 4, 5, 5, 5],
  [0 .. 4, 4, 5, 6, 6],
)
7:8:
7:7:
6:6:
5:5:
4:4:
4:3:
3:2:
2:1:
1:0:
0:0:
6
derecha diagonal diagonal diagonal abajo diagonal diagonal diagonal abajo  A T T C G
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Puntuación errónea en alineamiento

Notapor perl junior » 2017-01-26 12:49 @576

De acuerdo, tengo un par de preguntas.

¿Por qué si viene de la derecha añades a @camino abajo?

Además, si vemos el resultado, estaría mal, ¿no debería salir "TGCTTA"?

Y por otra parte, viendo ahora el código tanto en la línea 30 como en la 37 creo que lo tengo mal y que los índices empiecen por 1.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         for (my $i =1;$i<$len1;$i++){
  2.                
  3.                 $s[$i][0] =$s[$i-1][0];
  4.                
  5.                 $decisiones[$i][0]="abajo";
  6.         }
  7.  
  8.         for (my $j=1;$j<$len2;$j++){
  9.                
  10.                 $s[0][$j] =$s[0][$j-1];
  11.                
  12.                 $decisiones[0][$j]="derecha";
  13.         }
Coloreado en 0.007 segundos, usando GeSHi 1.0.8.4


Gracias de nuevamente, estaré atento.
perl junior
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Siguiente

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 1 invitado