• Publicidad

Extraer secuencias genómicas de forma ordenada

Perl aplicado a la bioinformática

Re: Extraer secuencias genómicas de forma ordenada

Notapor Alfumao » 2011-04-07 09:44 @447

Hola de nuevo.

Si transformo el programa de este post en una subrutina (que quiero usar para extraer sólo una secuencia de nombre pasado en la llamada a la subrutina):

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. &EXP_RET($linea);# ¿Cómo imprimo el resultado de esta expresión?
  2.                  #¿ $sec_result=&EXP_RET($linea);
  3.                  # print "$sec_result";?
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4


¿Cómo lo puedo sacar para incorporarlo al programa donde ésta está incluida?, es decir, ¿qué le asigno al comando "return" de la subrutina para que luego me devuelva la secuencia y la pueda imprimir? Llevo horas con esto y no lo puedo solucionar...

Aquí dejo el código de la subrutina derivado del tuyo. Éste funciona perfectamente si lo uso como programa individual, pero una vez pasado a subrutina e incluido en un CGI, no puedo sacar la secuencia resultado en pantalla. ¿Ves algún problema con este código?

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. sub EXP_RET {
  2.     my ($linea) = @_;
  3.  
  4.     # Lectura del fichero de secuencias
  5.     my %es_interesante;
  6.  
  7.     $es_interesante{$linea} = 1;       # la almacenamos en la memoria asociativa
  8.  
  9.     # Lectura del fichero FASTA
  10.     my $nombre_secuencia;              # Guarda el nombre de cada secuencia
  11.     my $secuencia = '';                # Contenido de la secuencia interesante
  12.  
  13.     open FASTA,  '<PC91.fasta'     or die "ERROR: $!\n";
  14.     open SALIDA, '>secuencias.txt' or die "ERROR: $!\n";
  15.  
  16.     while ( my $linea = <FASTA> ) {
  17.  
  18.         if ( $linea =~ /^>(\w+)/ ) {   # Comienzo de una secuencia nueva
  19.  
  20.             if ($secuencia) {          # Si tenemos una ya leída...
  21.                 procesar_secuencia( $nombre_secuencia, $secuencia )
  22.                     ;                  # la procesamos
  23.             }
  24.  
  25.             $nombre_secuencia
  26.                 = $1;                  # Nombre de la secuencia extraída desde la exp. reg.
  27.             $secuencia = '';           # Reiniciamos nuestra secuencia a capturar
  28.  
  29.             next;                      # pasamos a la siguiente línea
  30.         }
  31.  
  32.         if ( $es_interesante{$nombre_secuencia} )
  33.         {                              # Si estamos en una secuencia interesante
  34.             $secuencia .= $linea;      # la guardamos
  35.         }
  36.     }
  37.  
  38.     close FASTA;
  39.     close SALIDA;
  40.  
  41.     #if ($secuencia) { # Si aún queda alguna $secuencia por sacar
  42.     #    procesar_secuencia($nombre_secuencia, $secuencia); # la procesamos
  43.     #}
  44.  
  45.     #return ¿QUE VARIBLE ASIGNO PARA PODER IMPRIMIR LA SECUENCIA RESULTANTE?
  46.  
  47.     sub procesar_secuencia {
  48.         my ( $nombre, $secuencia ) = @_;
  49.         print ">$nombre\n$secuencia";
  50.         print SALIDA ">$nombre\n$secuencia";
  51.     }
  52. }
  53.  
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
Alfumao
Perlero nuevo
Perlero nuevo
 
Mensajes: 178
Registrado: 2009-12-10 11:20 @514

Publicidad

Re: Extraer secuencias genómicas de forma ordenada

Notapor explorer » 2011-04-07 10:17 @470

El resultado es un escalar: el nombre de la secuencia junto con un carácter de fin de línea, y seguido por la propia secuencia. Así que con guardar el resultado en una variable escalar, vale. La diferencia es que en mi versión, ese escalar se envía a la salida estándar y al controlador SALIDA. Tú solo tienes que ponerlo en el return().
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Extraer secuencias genómicas de forma ordenada

Notapor Alfumao » 2011-04-07 10:39 @486

Hola explorer,

Eso ya lo intenté de la siguiente forma,

en la subrutina puse como retorno:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my$ret = "$nombre_secuencia\n $secuencia";
  2. return ($ret);
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Y en el programa principal pongo:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my $sol= &EXP_RET($linea);
  2.  
  3. print"$sol";
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Me sale el nombre, pero en vez del valor de la secuencia, me sale el valor en blanco...
Alfumao
Perlero nuevo
Perlero nuevo
 
Mensajes: 178
Registrado: 2009-12-10 11:20 @514

Re: Extraer secuencias genómicas de forma ordenada

Notapor explorer » 2011-04-07 12:17 @553

Pues si sale en blanco, eso es que $secuencia no contiene nada.

Por lo que veo en el código, hay cosas muy raras... ¿qué hace $linea en la segunda línea de la subrutina?

Tienes además comentadas las líneas 41 a 43. Y no haces el filtrado de la secuencia que dices vas a pasar como argumento...

Te recomiendo que uses el truco de ir poniendo print() a lo largo de la subrutina para que vayas viendo cuándo y cómo va variando $secuencia.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Extraer secuencias genómicas de forma ordenada

Notapor Alfumao » 2011-04-07 15:19 @680

Hola explorer,

$linea = @_; es para recibir el argumento que le paso a la subrutina, ¿es que no se hace así?

Las líneas 41-43 están comentadas porque, si uso la subrutina a modo de programa (metiendo el valor de $linea desde <STDIN>, me devuelve el resultado correcto con esas líneas comentadas...)

En fin, seguiré probando el truco del print() a ver si al final lo consigo, pero es que cuando lo pregunté fue porque mi frustración y ofuscación llegaron al grado máximo.

Siento ser un pesado.

Un saludo
Alfumao
Perlero nuevo
Perlero nuevo
 
Mensajes: 178
Registrado: 2009-12-10 11:20 @514

Re: Extraer secuencias genómicas de forma ordenada

Notapor explorer » 2011-04-07 19:25 @850

¿El fichero de entrada siempre es el mismo? Porque eso es lo que veo en la línea 13.

Si en la línea 2 recibes el nombre de la secuencia, luego no lo usas para saber si la has encontrado o no, en la línea 25.

Habría que tener una bandera que se activase cuando hayamos encontrado el nombre que buscamos, y que permitiese almacenar en $secuencia las siguientes líneas, hasta el comienzo de la siguiente secuencia o el final de fichero.

¿El nombre a buscar es exacto al que encontraremos en el fichero o puede tener alguna variación?

Todas las líneas de SALIDA y de %es_interesante, sobran. Solo vamos a almacenar una $secuencia, así que con una sola variable, nos vale. La cuestión es saber cuándo empezar a capturar.

Una opción muy buena sería: leer todo el fichero en una sola variable escalar, y luego usar una expresión regular con el nombre de la secuencia a buscar, y que obtuviéramos como resultado, la secuencia que le sigue. En una sola línea de programa lo tendrías resuelto.

Bueno, sería más rápido usar index(), naturalmente:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use common::sense;              # sentido común
  3. use autodie;                    # sacrificio
  4. use File::Slurp 'slurp';        # chupa, chupa
  5. use open ':locale';
  6.  
  7. # Buscamos por una determinada secuencia, en un determinado fichero
  8. my $nombre_secuencia_a_buscar = 'jgi|Schco1|102220'; # 'jgi|Schco1|230323';
  9. my $fichero_de_secuencias     = 'PC91.fasta';
  10.  
  11. say "Buscando secuencia $nombre_secuencia_a_buscar en fichero $fichero_de_secuencias";
  12.  
  13. if (my $secuencia = buscar_secuencia($nombre_secuencia_a_buscar, $fichero_de_secuencias)) {
  14.     say '¡Encontrado!';
  15.     say $secuencia;
  16. }
  17. else {
  18.     say 'Lo siento, no lo he encontrado';
  19. }
  20.  
  21.  
  22. sub buscar_secuencia {
  23.     my ($nombre, $fichero) = @_;
  24.  
  25.     # Leemos el fichero
  26.     my $secuencias = slurp($fichero);
  27.  
  28.     # Buscamos el nombre de la secuencia
  29.     if (my $posición_inicial = index $secuencias, $nombre) {
  30.  
  31.         # Buscamos el final de la secuencia
  32.         my $posición_final = index $secuencias, "\n\n", $posición_inicial;
  33.                  
  34.         # Caso de que encontremos un final de fichero
  35.         $posición_final = length($secuencias)-1 if $posición_final == -1;
  36.  
  37.         # Extraemos toda la información de la secuencia
  38.         my $secuencia = substr $secuencias, $posición_inicial, $posición_final - $posición_inicial + 1;
  39.                  
  40.         # Extraemos solo la parte de la secuencia
  41.         $secuencia = substr $secuencia, index($secuencia, "\n") + 1;
  42.                  
  43.         return $secuencia;                              # Éxito...
  44.     }
  45.     else {
  46.         return;                                         # Fallo...
  47.     }
  48. }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4

Esta solución busca por un nombre de secuencia en un determinado fichero. Usa index() para localizar las distintas partes: "\n\n" como separador entre distintas secuencias, y "\n" como separador entre el nombre y la propia secuencia. Es una solución muchísimo más rápida que hacer bucles o usando expresiones regulares.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Extraer secuencias genómicas de forma ordenada

Notapor Alfumao » 2011-04-26 13:21 @598

Hola de nuevo explorer.

Cuando ya no sabía que más hacer, "usando el truco del print" me dí cuenta del fallo. El valor que pasaba a la subrutina era "1" en vez del nombre de la secuencia y era todo culpa de la siguiente línea:

chomp(my $linea=@_);

No se porqué en mi desesperación escribí semejante aberración, porque no necesitaba el chomp. Haciéndolo así es lógico que el resultado del valor pasado sea "1" (que es el valor escalar de un array de un solo elemento, si no me equivoco). Así que era tan fácil como escribir:

my($linea)=@_;

Que es como "reciben" los datos las subrutinas...

Siento mucho las molestias causadas.

Un saludo
Alfumao
Perlero nuevo
Perlero nuevo
 
Mensajes: 178
Registrado: 2009-12-10 11:20 @514

Anterior

Volver a Bioinformática

¿Quién está conectado?

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