• Publicidad

Puntuación y consenso

Perl aplicado a la bioinformática

Puntuación y consenso

Notapor perl junior » 2017-01-08 19:07 @838

Hola, estoy realizando el siguiente script:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #Enumeración de motivos
  2.  
  3.  
  4.  
  5. sub motifs {
  6.         my $t =$_[0];
  7.         my $k = $_[1];
  8.         my @s =@{$_[2]};#Lista pasada por referencia
  9.         my @secuencias = @{$_[3]};
  10.        
  11.         my @motifs=();
  12.         for (my $i=0; $i< scalar @s; $i++) {
  13.                 my $indice =$s[$i];
  14.                 my $secuencia =$secuencias[$i];
  15.                 my $k_mero = substr ($secuencia,$indice,$k);
  16.                
  17.                 push @motifs,$k_mero;
  18.         }
  19.         return @motifs;
  20. }
  21.  
  22.  
  23.  
  24.  
  25. sub profile{
  26.         my ($t,$k,@motifs) =@_;
  27.         my %tabla =();
  28.         my @bases = ("A","C","G","T");
  29.         foreach my $base (@bases){
  30.                 my @lista_ceros =();
  31.                 for (my $i =0; $i<$k;$i++){
  32.                         push @lista_ceros,0;
  33.                 }
  34.                 $tabla{$base} = \@lista_ceros;
  35.         }
  36.        
  37.     return  %tabla;
  38.          #LO primero que tengo que hacer es inicializar la tabla, un buen punto de partida es poner ocho ceros
  39.          #declarar una tabla con ocho ceros
  40.        
  41. }
  42.  
  43. my $t1 = 5;
  44. my $k1 = 8;
  45. my @s1 =(6,20,23,10,11);
  46. my @secuencias1 = (
  47.     "CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA",
  48.     "GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG",
  49.     "TAGTACCGAGACCGAAAGAAGTATACAGGCGT",
  50.     "TAGATCAAGTTTCAGGTGCACGTCGGTGAACC",
  51.     "AATCCACCAGCTCCACGTGCAATGTTGGCCTA",
  52. );
  53.  
  54.  
  55.  
  56. my  @resultado1 = motifs($t1,$k1,\@s1,\@secuencias1);
  57. print join (" ", @resultado1, "\n");
  58. my %tabla1 = profile($t1,$k1, @resultado1);
  59. foreach my $base (keys%tabla1){
  60.         my @lista = @{$tabla1{$base}};#Me quedo con una copia de la lista a la que referencia
  61.         print "Base" ,$base," ",join(" ",@lista), "\n";
  62.  
  63. }
Coloreado en 0.020 segundos, usando GeSHi 1.0.8.4


Mi duda es que una vez llegado a este punto me piden que calcule la máxima puntuación que devuelve la matriz y que devuelva la secuencia consenso, pero no consigo enfocarlo.

Si alguien me pudiera dar algunas indicaciones para hacer lo primero...

Muchas gracias.
perl junior
Perlero Nuevo
Perlero Nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Publicidad

Re: Puntuación y consenso

Notapor explorer » 2017-01-10 17:23 @766

Supongo que tendrás que recorrer la matriz por filas y columnas, y recordar la posición que contenga el valor más pequeño.
JF^D Perl Programming
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 13907
Registrado: 2005-07-24 18:12 @800
Ubicación: Madrid, España

Re: Puntuación y consenso

Notapor perl junior » 2017-01-10 18:27 @811

Quizás no me haya explicado bien.

El resultado que devuelve la segunda función es este:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
BaseC 1 3 3 1 1 0 1 0
BaseT 4 1 1 0 0 0 3 0
BaseG 0 0 0 0 4 5 1 5
BaseA 0 1 1 4 0 0 0 0
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Siendo las bases las claves, y el valor asociado a la clave la lista de números, entonces tengo que acceder a la primera posición de la lista de cada valor para cada clave y calcular el mayor y así para todas las columnas. Por ejemplo, de la primera columna el mayor es el 4.

El problema, y es en lo que me estoy atascando, es cómo acceder a la primera posición de la lista que es el valor de la clave, para cada una de las claves.

¡¡Muchas gracias!!
perl junior
Perlero Nuevo
Perlero Nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509

Re: Puntuación y consenso

Notapor explorer » 2017-01-11 03:54 @204

Pues eso... recorrer, para todas las columnas, y por cada fila de la tabla, todos los valores, y llevar la cuenta del máximo encontrado en cada columna:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/env perl
  2. #
  3. # Enumeración de motivos
  4.  
  5. my @bases = qw(A C G T);
  6.  
  7. sub motifs {
  8.     my($longitud, $inicios_ref, $secuencias_ref) = @_;
  9.  
  10.     my @motifs;
  11.  
  12.     for (my $i = 0; $i < scalar @$inicios_ref; $i++) {
  13.         push @motifs, substr $secuencias_ref->[$i], $inicios_ref->[$i], $longitud;
  14.     }
  15.  
  16.     return @motifs;
  17. }
  18.  
  19. sub profile{
  20.     my ($longitud, $motifs_ref) = @_;
  21.     my %tabla;
  22.  
  23.     # Inicialización
  24.     foreach my $base (@bases) {
  25.         $tabla{ $base } = [ (0) x $longitud ];
  26.     }
  27.  
  28.     # Estadísticas
  29.     foreach my $motifs (@$motifs_ref) {
  30.  
  31.         for (my $i = 0; $i < $longitud; $i++) {
  32.  
  33.             $tabla{ substr($motifs, $i, 1) }->[$i]++;           # en la fila de la base, columna $i-ésima, sumamos uno
  34.         }
  35.     }
  36.  
  37.     return \%tabla;
  38. }
  39.  
  40. my $longitud = 8;
  41. my @inicios  = (6, 20, 23, 10, 11);
  42.  
  43. my @secuencias = qw(
  44.     CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA,
  45.     GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG,
  46.     TAGTACCGAGACCGAAAGAAGTATACAGGCGT,
  47.     TAGATCAAGTTTCAGGTGCACGTCGGTGAACC,
  48.     AATCCACCAGCTCCACGTGCAATGTTGGCCTA,
  49. );
  50.  
  51.  
  52. my @motifs = motifs($longitud, \@inicios, \@secuencias);
  53. print "@motifs\n";
  54.  
  55. my $tabla_ref = profile($longitud, \@motifs);
  56.  
  57. foreach my $base (sort @bases) {
  58.     print "Base $base: @{ $tabla_ref->{$base} }\n";
  59. }
  60.  
  61. # Calcular máximos
  62. my @maximos_columnas;
  63.  
  64. for (my $i = 0; $i < $longitud; $i++) {         # para todas las columnas
  65.  
  66.     my $maximo = 0;                             # vamos a encontrar un máximo
  67.  
  68.     foreach my $base (@bases) {                 # para todas las filas
  69.  
  70.         my $veces = $tabla_ref->{$base}->[$i];  # extraemos el número correspondiente a esa fila y columna
  71.  
  72.         if ($maximo < $veces) {                 # si el número es mayor que el máximo,
  73.             $maximo = $veces;                   # nuevo máximo encontrado, lo recordaremos
  74.         }
  75.     }
  76.  
  77.     push @maximos_columnas, $maximo;            # guardamos el máximo encontrado en la lista
  78. }
  79.  
  80. print "Máx.  : @maximos_columnas\n";
  81.  
Coloreado en 0.028 segundos, usando GeSHi 1.0.8.4
La salida es:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
TCTCGGGG CCAAGGTG TACAGGCG TTCAGGTG TCCACGTG
Base A: 0 1 1 4 0 0 0 0
Base C: 1 3 3 1 1 0 1 0
Base G: 0 0 0 0 4 5 1 5
Base T: 4 1 1 0 0 0 3 0
Máx.  : 4 3 3 4 4 5 3 5
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
JF^D Perl Programming
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 13907
Registrado: 2005-07-24 18:12 @800
Ubicación: Madrid, España

Re: Puntuación y consenso

Notapor perl junior » 2017-01-11 08:18 @387

Muchas gracias, explorer. Como el objetivo era obtener la puntuación consenso y la puntuación máxima dejo aquí el programa por si hubiera alguien que quisiera saberlo

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #Enumeración de motivos
  2.  
  3. sub motifs {
  4.         my $t =$_[0];
  5.         my $k = $_[1];
  6.         my @s =@{$_[2]};#LIsta pasada por referencia
  7.         my @secuencias = @{$_[3]};
  8.        
  9.         my @motifs=(); # lista vacía  donde iré metiendo las motivos de tamaño k que extraiga de cada secuencia
  10.         for (my $i=0; $i< scalar @s; $i++) {
  11.                 my $indice =$s[$i];
  12.                 my $secuencia =$secuencias[$i];
  13.                 my $k_mero = substr ($secuencia,$indice,$k);
  14.                
  15.                 push @motifs,$k_mero;
  16.         }
  17.         return @motifs;
  18. }
  19.  
  20. my $t1 = 5;
  21. my $k1 = 8;
  22. my @s1 =(6,20,23,10,11);
  23. my @secuencias1 = (
  24.     "CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA",
  25.     "GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG",
  26.     "TAGTACCGAGACCGAAAGAAGTATACAGGCGT",
  27.     "TAGATCAAGTTTCAGGTGCACGTCGGTGAACC",
  28.     "AATCCACCAGCTCCACGTGCAATGTTGGCCTA",
  29. );
  30.  
  31.  
  32. sub profile{
  33.         my ($t,$k,@motifs) =@_;
  34.         my %tabla_profile =();
  35.         my @bases = ("A","C","G","T");
  36.         foreach my $base (@bases){
  37.                 my @lista_ceros =();
  38.                 for (my $i =0; $i<$k;$i++){
  39.                         push @lista_ceros,0;
  40.                 }
  41.                 $tabla_profile{$base} = \@lista_ceros;
  42.                 #print @{$tabla_profile{"A"}};
  43.                
  44.                
  45.         }
  46.         foreach my $motif(@motifs){ #Queremos acceder a los motivos,vamos a recorrerlos/Al ser irrevelante el índice de la secuencia usar foreach
  47.                 for (my $i=0;$i<$k;$i++){
  48.                         my $base = substr($motif,$i,1);#Vamos a extraer cada base de la secuencia, en este caso la posicion es relevante
  49.                         #$base, es la clave de acceso a la tabla/ $tabla_profile{$base}= referencia a la lista
  50.                         #@{$tabla{$base}}, hago una copia
  51.                         #${$tabla{$base}}, acceso al valor de esa lista
  52.                         my $ref_lista = $tabla_profile{$base}; #COn esto tengo el identificador de la lista
  53.                        
  54.                         ${$ref_lista}[$i] = ${$ref_lista}[$i]+1;
  55.                        
  56.                      
  57.                 }
  58.        
  59.                
  60.         }
  61.     return  %tabla_profile;
  62.          #LO primero que tengo que hacer es inicializar la tabla, un buen punto de partida es poner ocho ceros
  63.          #declarar una tabla con ocho ceros
  64.        
  65. }
  66.  
  67.  
  68.  
  69.  
  70.  
  71. my  @resultado1 = motifs($t1,$k1,\@s1,\@secuencias1);
  72. print join (" ", @resultado1, "\n");
  73. my %tabla= profile($t1,$k1, @resultado1);
  74. foreach my $base (keys%tabla){
  75.         my @lista = @{$tabla{$base}};#Me quedo con una copia de la lista a la que referencia
  76.         print "Base" ,$base," ",join(" ",@lista), "\n";
  77.  
  78. }
  79. #if($contador < $num_max_apariciones){$num_max_apariciones =$contador;}
  80.                 #$num_min_apariciones=$num_max_apariciones
  81.  
  82.  
  83. sub puntuacion {
  84.         my($k,%profile)=@_;
  85.         my @maximos_columnas= ();
  86.         my $consenso="";
  87.         for( my $i=0;$i<$k;$i++){                        # para todas las columnas (en este caso si que interesa la posicion)
  88.                 my $maximo= 0 ;                              #vamos a encontrar un máximo de esas columnas
  89.                 foreach my $base (keys%tabla){               # para todas las filas
  90.                         @lista = @{$profile{$base}};
  91.                        
  92.                         my $veces= $lista[$i];
  93.                                                                  #si el número es mayor que el máximo
  94.                         if ($maximo < $veces){
  95.                                 $maximo=$veces;
  96.                                
  97.                        
  98.                        
  99.                         }
  100.                
  101.                
  102.                
  103.                 }
  104.                  
  105.          push @maximos_columnas,$maximo;                  #guardamos el máximo encontrado en una lista
  106.          foreach my $base (keys%tabla){                    # para todas las filas
  107.                         @lista = @{$profile{$base}};
  108.                         my$repeticiones= $lista[$i];
  109.                 if ($repeticiones eq $maximo){$consenso=$consenso.$base;}
  110.              
  111.              
  112.              }
  113.         }
  114.         $puntuacion_total=0;
  115.         foreach $puntuacion (@maximos_columnas){
  116.         $puntuacion_total =$puntuacion+$puntuacion_total;
  117.    
  118. }
  119.         return  $consenso,$puntuacion_total
  120. }
  121.  
  122. @resultado=puntuacion($k1,%tabla);
  123.  
  124. print "La mejor puntuación obtenida es: ",$resultado[1],"\n";
  125. print "La secuencia consenso es: ",$resultado[0];
Coloreado en 0.028 segundos, usando GeSHi 1.0.8.4


Salida
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
TCTCGGGG CCAAGGTG TACAGGCG TTCAGGTG TCCACGTG
BaseC 1 3 3 1 1 0 1 0
BaseT 4 1 1 0 0 0 3 0
BaseA 0 1 1 4 0 0 0 0
BaseG 0 0 0 0 4 5 1 5
La mejor puntuación obtenida es: 31
La secuencia consenso es: TCCAGGTG
 
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
perl junior
Perlero Nuevo
Perlero Nuevo
 
Mensajes: 10
Registrado: 2017-01-03 11:13 @509


Volver a Bioinformática

¿Quién está conectado?

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