• Publicidad

Lectura y extracción de datos de multiFASTA

Perl aplicado a la bioinformática

Lectura y extracción de datos de multiFASTA

Notapor Nuiter » 2019-06-13 08:14 @385

Hola a todos, este es mi primer mensaje y aprovecho para presentarme.

Lo primero decir que en programación soy cinturón blanco. Y estoy haciendo un curso y la tarea final no consigo sacarla.

No pretendo que nadie me escriba el código porque de nada me serviría, pero agradecería mucho que me indicasen dónde me equivoco y el motivo del error.

Estoy realizando un pequeño script para leer de un multiFASTA, con una particularidad, y es que al final de la secuencia proteica tiene un "*"; y posteriormente realizar un blast.

Para ello he pensado que lo mejor sería quedarme solo con el ID de cada proteína y con la secuencia, para luego pasar esta por el blast, además de eliminar el "*".

Pero llevo mucho tiempo atascado y no consigo avanzar y además veo que pierdo el tiempo.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2.  
  3. die "USAGE: capturador.pl file.fa\n" unless ( $ARGV[0] );
  4.  
  5. open( QUERY, "$ARGV[O]" );
  6. while (<QUERY>) {
  7.  
  8.     #print; primer print control.
  9.     if (/^>(\S+)/) {
  10.         @seq = ();
  11.         $id  = $1;
  12.         print "$id\n";
  13.  
  14.     }
  15.     else { push( @seq, $_ ) }
  16.  
  17.     $juntos = join( "", @seq );
  18.  
  19.     #print "$juntos";
  20.     print "@seq\n";
  21. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

Este es el código que he escrito, pero hace cosas extrañas: me da el ID como pido, pero la secuencia me la da por partes. Primero, la primera línea, después las dos primeras, tres primeras, hasta acabar, y por otra parte no sé cómo quitarle el '*', pero eso me seguiré peleando.

Adjunto el FASTA (la extensión para mi es .fa, pero la he tenido que cambiar para subirla).

Gracias mil, por adelantado.
Adjuntos
fastadeldemon.txt
(2.31 KiB) 333 veces
Nuiter
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2019-06-11 17:52 @786

Publicidad

Re: Lectura y extracción de datos de multiFASTA

Notapor explorer » 2019-06-13 12:19 @555

Bienvenido a los foros de Perl en Español, Nuiter.

Por el foro de Bioinformática encontrarás algunos ejemplos de cómo leer archivos FASTA.

La idea es la siguiente: vas leyendo por líneas. Si la línea comienza por '>' entonces sabes que has llegado a una nueva secuencia proteica. Entonces, debes pintar (o grabar o lo que sea) la secuencia anterior (si la hubiera), y reiniciar el búfer (el array @seq en tu código). Al final de la lectura de todas las líneas, debes comprobar si queda algo en el búfer, y pintarlo (este es el caso de la última secuencia del archivo).

Debes prestar atención al formato del archivo: veo que hay líneas en blanco entre las secuencias. Esas líneas deberías descartarlas. Por ejemplo con un next if /^\s*$/;

Para quitar el '*' final, lo puedes hacer después del join() en el momento de tener toda la secuencia. Es decir, si estamos completamente seguros de que TODAS las secuencias tienen un '*' y SIEMPRE terminan con él, entonces con un simple chop() lo quitas. Si no estamos seguros de que se cumplan esas condiciones, mejor usar una expresión regular de sustitución: $juntos =~ s{[*]$}{};

Un detalle importante: dentro del bucle no estás quitando los caracteres de final de línea que hay en el archivo. A la hora de hacer el push() se están grabando junto con la secuencia proteica. No sé si eso es lo que quieres, ya que la mayoría querrá tener solo la secuencia. Si quieres quitar los finales de línea, te basta con poner un chomp; después de la lectura de la línea.

¡Ánimo!
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

Re: Lectura y extracción de datos de multiFASTA

Notapor Nuiter » 2019-06-13 13:54 @621

Muchas gracias, explorer. Gracias a tus sugerencias y observaciones he conseguido que el código haga lo que necesito.

Ahora a pelearme con el BLAST.

Un abrazo.

Comparto los cambios que he hecho, por si a otro usuario pudieran servirle en un futuro.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2.  
  3. die "USAGE: capturador.pl file.fa\n" unless ( $ARGV[0] );
  4.  
  5. open( QUERY, "$ARGV[O]" );
  6.  
  7. while (<QUERY>) {
  8.     next if /^\s*$/;
  9.     my @seq = ();
  10.     if (/^>(\S+)/) {
  11.         my $id = $1;
  12.         print "$id\n";
  13.     }
  14.     else {
  15.         push( @seq, $_ );
  16.         $juntos = join( "", @seq );
  17.         $juntos =~ s/[*]//g;
  18.         print $juntos;
  19.     }
  20. }
  21.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
Nuiter
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2019-06-11 17:52 @786

Re: Lectura y extracción de datos de multiFASTA

Notapor explorer » 2019-06-15 13:56 @622

Me temo que ese código no está resolviendo el problema...

La sola presencia de

my @seq = ();

indica que algo (bastante) está mal.
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

Re: Lectura y extracción de datos de multiFASTA

Notapor Nuiter » 2019-06-16 12:38 @568

Con my @seq=(); quería inicializar el array en cada paso del bucle. ¿No es correcto? Mi problema al quitarlo, es que, el resultado vuelve a ser el anterior, es decir, empieza una y otra vez a copiar líneas, de manera acumulativa, en cada iteración. Un desastre, vaya.

Por otra parte estoy intentando que me haga un blastp de manera remota. He incluido esta línea:

blastp -query INT -db swissprot -out result.fas -remote

pero no me lo hace:

Can't locate object method "out" via package "result" (perhaps you forgot to load "result"?) at ./capturador.pl line 24, <QUERY> line 42.

Cuando utilizo esta expresión en un terminal me la ejecuta pero me da fallos de conexión:

Error: (303.7) [URL_Connect] Failed to connect to http://www.ncbi.nlm.nih.gov:443: Not supported
Error: (303.7) [URL_Connect] Failed to connect to http://www.ncbi.nlm.nih.gov:443: Not supported
Error: (303.7) [URL_Connect] Failed to connect to http://www.ncbi.nlm.nih.gov:443: Not supported
Error: (310.5) [blast4] Service not found
Error: (315.2) CConn_Streambuf::CConn_Streambuf(): NULL connector: Unknown
.
.
.

Pero al menos me ejecuta la línea.

Mil gracias.
Nuiter
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2019-06-11 17:52 @786

Re: Lectura y extracción de datos de multiFASTA

Notapor explorer » 2019-06-17 07:37 @359

No se trata de que "quites" la línea por que sí. Se trata de que encuentres el lugar donde debe ir.

Primero: el uso de "my" es para declarar variables locales. @seq es una variable que debe ser conocida dentro del bucle y nos servirá para almacenar la secuencia que vamos leyendo. Por eso, lo mejor es declarar esa variable fuera del bucle:

my @seq;

Segundo: el bucle while va leyendo líneas. Queremos que cada línea se vaya acumulando (con un push) a @seq. Solo cuando se inicia una nueva secuencia es cuando debemos sacar la secuencia anterior (almacenada en @seq), y volver a iniciar @seq (@seq=()). Esto hay que repetirlo al salir del bucle, ya que puede quedar alguna secuencia por procesar.

Sería algo así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. #
  3. # Ejemplo de lectura de archivos MultiFASTA
  4. #
  5. # Los archivos MultiFASTA tienen esta estructura:
  6. #
  7. # -----------------------------------------8<-------------------------------------------
  8. # >sp|P0AA16|OMPR_ECOLI Transcriptional regulatory protein OmpR OS=Escherichia coli (strain K12) GN=ompR PE=1 SV=1
  9. # MQENYKILVVDDDMRLRALLERYLTEQGFQVRSVANAEQMDRLLTRESFHLMVLDLMLPG
  10. # EDGLSICRRLRSQSNPMPIIMVTAKGEEVDRIVGLEIGADDYIPKPFNPRELLARIRAVL
  11. # RRQANELPGAPSQEEAVIAFGKFKLNLGTREMFREDEPMPLTSGEFAVLKALVSHPREPL
  12. # SRDKLMNLARGREYSAMERSIDVQISRLRRMVEEDPAHPRYIQTVWGLGYVFVPDGSKA*
  13. #
  14. # >sp|P10398|ARAF_HUMAN Serine/threonine-protein kinase A-Raf OS=Homo sapiens GN=ARAF PE=1 SV=2
  15. # ...
  16. # -----------------------------------------8<-------------------------------------------
  17. #
  18. # En este ejemplo, leemos las secuencias y las sacamos a la salida estándar
  19. #
  20. # Joaquín Ferrero. 20190617
  21. #
  22.  
  23. use v5.14;
  24. use autodie;                                            # nos facilita el detectar fallos de entrada/salida
  25.  
  26. # Comprobacion de los argumentos
  27. # Necesitamos el nombre de un archivo, y con tamaño distinto de cero
  28. die "USO: $0 archivo.fa\n" unless ($ARGV[0] and -s $ARGV[0]);
  29.  
  30.  
  31. # Procesamiento
  32. my $id;                                                 # identificador de cada secuencia
  33. my @seq;                                                # almacén para cada secuencia
  34. open my $FASTA, '<', $ARGV[0];
  35. while (<$FASTA>) {
  36.     next if /^\s*$/;                                    # descartamos las líneas vacías
  37.  
  38.     chomp;                                              # quitamos los caracteres de fin de línea
  39.  
  40.     if (/^>(\S+)/) {                                    # encontrada una nueva secuencia
  41.  
  42.         if ($id) {                                      # ¿hay una secuencia anterior?
  43.             say $id;                                    # sí: la sacamos en pantalla
  44.            
  45.             my $juntos = join '', @seq;                 # componemos la secuencia
  46.             $juntos =~ s/[*]//g;                        # quitamos todos los '*'
  47.             say $juntos;
  48.  
  49.             undef @seq;                                 # vaciamos el almacén para la siguiente secuencia
  50.         }
  51.  
  52.         $id = $1;                                       # recordamos el id de la nueva secuencia
  53.     }
  54.     else {                                              # encontrada una línea de secuencia
  55.         push @seq, $_;                                  # la almacenamos
  56.     }
  57. }
  58. close $FASTA;
  59.  
  60. if ($id) {                                              # ¿hay una secuencia anterior?
  61.     say $id;                                            # sí: la sacamos en pantalla
  62.    
  63.     my $juntos = join '', @seq;                         # componemos la secuencia
  64.     $juntos =~ s/[*]//g;                                # quitamos todos los '*'
  65.     say $juntos;
  66. }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
Con este programa, la salida son las cuatro secuencias contenidas en el archivo de ejemplo, con una línea por secuencia.

En cuanto al problema que tienes para ejecutar el blastp, como no vemos el código, no sabemos qué puede ser. Pero el mensaje de error te dice que está en la línea 24.

Y si lo ejecutas desde la línea de comandos, es posible que no le estés entregando los datos en el formato que la web espera.
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

Re: Lectura y extracción de datos de multiFASTA

Notapor Nuiter » 2019-06-21 13:37 @609

Hola de nuevo, estimado explorer y resto de la comunidad.

Después de haber corregido mi código en base a tus sugerencias y a nueva información que me han dado para acabar la tarea creo que casi-casi lo tengo.

Los cambios son de sintaxis (todo aquello que entendía de tu código, explorer, no entendía por qué a un filehandler se le asignaba una variable con "$", por ejemplo), y el más radical e importante es que hay que hacer varios archivos y varias consultas con blastp, ya que no se puede pasar un multiFASTA a blastp. De esto último me he ido dando cuenta con las pruebas.

Pues bien, mi código casi hace lo que le pido, es decir: me almacena para cada secuencia su ID, y la secuencia aminoacídica, pero no lo hace bien, ya que es justo a continuación, sin retorno de línea como le indico hasta en dos ocasiones.

Aún así parece que el blastp en algunas de ellas sí lo acepta aunque no he conseguido corregirlo.

Y por otro lado no consigo que pase por el blastp los archivos que he generado, dentro del bucle while.

Me dice esto: Illegal division by zero at ./capturador2.pl line 29, <QUERY> line 2.

Lo cual es muy desconcertante. El código que me queda al final es algo así, en la ruta de la base de datos (es una ruta local, del cluster que utilizamos en el curso).

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use autodie;
  3.  
  4. die "USAGE: capturador.pl file.fa\n" unless ( $ARGV[0] and -s $ARGV[0] );
  5. my $id;
  6. my $name;
  7. my @seq;
  8. open QUERY, '<', $ARGV[O];
  9.  
  10. while (<QUERY>) {
  11.     next if /^\s*$/;
  12.     chomp;
  13.     if (/^>(\S+)/) {
  14.         $id   = $1;
  15.         $name = $id;
  16.         $name =~ s/[|]/_/g;
  17.         open INT, ">workfile$name.fa";
  18.         $id =~ s/[\n]//g;
  19.         print INT "$id\n";
  20.     }
  21.     else {
  22.         push( @seq, $_ );
  23.         $juntos = join( "", @seq );
  24.         $juntos =~ s/[*]//g;
  25.         $juntos =~ s/[\n]//g;
  26.         print INT "\n$juntos";
  27.         undef @seq;
  28.         system(blastp -query workfile$name.fa -db /mnt/beegfs/bioinfo_db/Blast/swissprot_db -out result $name.fas);
  29.     }
  30.  
  31.     # system (blastp -query workfile$name.fa -db /mnt/beegfs/bioinfo_db/Blast/swissprot_db -out result$name.fas);
  32. }
  33.  
  34. close(QUERY);
  35. close(INT);
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Si os dais cuenta he intentado ejecutar el blastp, dentro y fuera del else{}, pero el fallo es el mismo.

A modo de resumen, mis problemas son dos:
  • No consigo que los archivos generados queden de la forma "ID \n secuencia", sino que quedan "IDsecuencia".
  • El blastp, no consigo que lo ejecute.
Muchas gracias, explorer. He aprendido mucho de todo lo que me has sugerido/indicado.
Nuiter
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2019-06-11 17:52 @786

Re: Lectura y extracción de datos de multiFASTA

Notapor explorer » 2019-06-21 21:11 @924

Por favor, si no entiendes algo, pregunta.

Sigues trabajando con tu código porque no entiendes el mío. Y por lo tanto, sigue teniendo los mismos errores graves que había al principio (no procesas la última secuencia del archivo, se procesa cada secuencia de forma incremental...).

Lo más grave es que no entiendes el momento en que se debe producir la salida de $juntos. Se supone que esa variable debe contener toda la secuencia, pero si te fijas bien, imprimes $juntos por cada línea del archivo multifasta que contiene parte de la secuencia. Entonces @seq va almacenando las distintas partes de la secuencia, pero a continuación, las sacas de forma repetida.

A ver si con un ejemplo se entiende mejor:

Supongamos un archivo con el siguiente aspecto:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
>Nombre
AAAAAAAAAA
BBBBBBBBBB
CCCCCCCCCC
DDDDDDDDD*
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

Cuando el bucle while lea la línea "AAAA...", @seq la almacenará (por efecto del push()). A continuación, $juntos = "AAA...", y lo imprimes hacia INT. ¿Y llamas a blastp con ese trozo?

En la siguiente vuelta del bucle, $_ vale "BBBBBB...". @seq agrega un elemento más. Y ahora haces $juntos = "AAAAABBBBB". Lo imprimes a INT (que ahora contiene "\nAAAAAAAA\nAAAAAAABBBBBBB"). Y llamas a blastp para que procese eso. ¿Ves que algo anda mal?

Al final, el archivo asociado a INT contiene esto:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
AAAAAAAAAA
AAAAAAAAAABBBBBBBBBB
AAAAAAAAAABBBBBBBBBBCCCCCCCCCC
AAAAAAAAAABBBBBBBBBBCCCCCCCCCCDDDDDDDDD
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Y mientras has llamado a blastp 4 veces...

El system() te falla porque los argumentos deben estar entrecomillados. Como hay variables escalares en su interior, entonces ese entrecomillado deberá ser doble. Así:

system("blastp -query workfile$name.fa -db /mnt/beegfs/bioinfo_db/Blast/swissprot_db -out result $name.fas");

El "\n" extra que te aparece es porque en efecto hay dos "\n". Uno en la línea 19 y otro en la 26. Te sobra uno. O ponerlos cuando realmente sabemos cómo será el contenido del archivo.
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


Volver a Bioinformática

¿Quién está conectado?

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