Foro - Perl en Español

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.


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.

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!

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.

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

La sola presencia de

my @seq = ();

indica que algo (bastante) está mal.

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.

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í:
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.

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).


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.

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:

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:
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.