Foro - Perl en Español

por **Alfumao** » 2012-06-11 03:11 @174

Buenas a todos,

El problema que se me presenta hoy me está costando solucionarlo y quería saber si hay una forma sencilla de hacerlo usando Perl (parece que últimamente hay una epidemia de usuarios de Bash en el campo de la bioinformática...)

He de filtrar en un archivo FASTQ las secuencias que no pasan el corte de calidad, eliminándolas del archivo. Estas secuencias se diferencian por tener un carácter "#" delante de la última línea, pero el problema es que son 4 líneas las que forman la secuencia y no sé muy bien cómo hacer para decirle a Perl que quite la línea que empieza por "#" (esta es fácil: /^#/) y las 3 anteriores a ella (esto es lo que no sé cómo hacer, porque el comando grep en Bash tiene las opciones -B y -A para quitar líneas anteriores y posteriores, pero en Perl no encuentro esa posibilidad).

Aquí os pongo el formato FASTQ normal :

Sintáxis: [ Descargar ] [ Ocultar ]

Using text Syntax Highlighting

@ILLUMINA-GA:48:FC64L4UAAXX:1:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G
Coloreado en 0.000 segundos,  usando GeSHi 1.0.8.4

Y aquí el formato a filtrar:

Sintáxis: [ Descargar ] [ Ocultar ]

Using text Syntax Highlighting

@ILLUMINA-GA:48:FC64L4UAAXX:1:1:10257:1000 1:N:0:CGATGT

NAGTTATAGGCAGAAAGTATTCCCAGCACCTTAACCCTTACACAAGGAAAAGTTCTTTCTTTAAAAACACATTTATTTGTGTGTGTGTGT

+

#(.((''-+,:1:::888858889379898@@@@@@@@@@@@@5@5775788585850080001770500@@@@@###############
Coloreado en 0.000 segundos,  usando GeSHi 1.0.8.4

Y aquí el script de Bash que supuestamente lo haría (según la casa comercial que vende los secuenciadores que producen estos archivos FASTQ) pero con sus expresiones regulares me hago un lío considerable porque parecen como las de Perl, pero noooooooo, ¡ja,ja,ja!

Sintáxis: [ Descargar ] [ Ocultar ]

Using bash Syntax Highlighting

grep -A 4 '^@.* [^:]*:N:[^:]*:'
Coloreado en 0.002 segundos,  usando GeSHi 1.0.8.4

Gracias por adelantado.

Por que desconocen BioPerl. O porque saben qué es pero tienen miedo a meterse con él. O porque son incompetentes que dependen sólo de las herramientas prefabricadas que compran: prefieren comprarse un caro programa de gestión de secuencias, que no saber que su problema se resuelve en unos minutos con BioPerl... pero claro... tendrían que *pensar*. Y cuando su maravilloso programa no hace todo lo que quieren, echan mano de lo que único que aprendieron, como por ejemplo Bash.

Esa es la diferencia entre un biólogo, un informático y un bioinformático: se supone que los biólogos no quieren saber nada de informática, y que sus problemas los resuelvan "mágicamente" las máquinas, aunque sea comprando programas comerciales. El informático no sabe nada más que tiene que manejar archivos con muchas letras. Y el bioinformático debe saber de todo eso a la vez.

Y es justo en el *debe* lo que distingue a los buenos de los malos: o resuelve un problema comprando una cara herramienta, o hace un script de mil líneas en Bash, o un programa Perl de un centenar de líneas en puro Perl, o un programa con BioPerl de apenas una docena de líneas.

Alfumao escribiste:He de filtrar en un archivo FASTQ las secuencias que no pasan el corte de calidad, eliminándolas del archivo. Estas secuencias se diferencian por tener un carácter "#" delante de la última línea, pero el problema es que son 4 líneas las que forman la secuencia y no sé muy bien cómo hacer para decirle a Perl que quite la línea que empieza por "#" (esta es fácil: /^#/) y las 3 anteriores a ella.

Esta es una manera...

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#!/usr/bin/env perl
use Tie::File;
 
tie my @lineas, 'Tie::File', 'kk.txt';  # consideramos el archivo como un array
 
for (my $i = 0; $i < @lineas; $i++) {   # recorremos todas las líneas
   if ($lineas[$i] =~ /^@/) {           # detectado comienzo de secuencia
      if ($lineas[$i+3] =~ /^#/) {      # detectada mala calidad en la cuarta línea
          splice @lineas, $i, 4;        # rebanamos las 4 líneas
          redo;                         # reinicimos el bucle, en la posición $i
      }
   }
}
 
untie @lineas;
Coloreado en 0.002 segundos,  usando GeSHi 1.0.8.4

por **Alfumao** » 2012-06-11 07:41 @362

Muchas gracias una vez más, explorer.

Respecto a lo del Bioperl, he de decirte que está bastante bien, para lo que está bastante bien... pero hay veces que si te sales de los supuestos para los que se creo algún módulo concreto, es más rápido hacerte un parser en Perl puro (sobre todo después de la cantidad de cosas que se aprenden en este foro) que intentar pegarte con el código del módulo de Bioperl o preguntar en su mailing list... te lo digo por experiencia, parece que si encuentras alguna "debilidad" en "su modulo" algunos de los autores se agarran un pique importante y se ponen bastante bordes

He intentado introducir tu código en un programa más grande y no parece que funciona... no acierto a ver dónde metí la pata, porque tampoco me da ningún error, solo me da un archivo de 0 kb en donde habría de salir el archivo procesado...

Este es el programa donde lo inserté:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#!usr/bin/perl -w
use Tie::File;
 
#usage program.pl [full dir path]
 
my $path = shift @ARGV;
print "$path\n";
chdir $path or die "ERROR: Unable to enter $path: $!\n";
opendir( TEMP, "." );
my @files = readdir(TEMP);
closedir TEMP;
 
#print"@files\n";
 
for my $file (@files) {
 
    if ( $file =~ /(\w+)\.fastq/ ) {
 
        print "This is $file processing...\n";
 
        tie my @lineas, 'Tie::File', '$file';    # consideramos el archivo como un array
 
        for ( my $i = 0; $i < @lineas; $i++ ) {  # recorremos todas las líneas
            if ( $lineas[$i] =~ /^@/ ) {         # detectado comienzo de secuencia
                if ( $lineas[ $i + 3 ] =~ /^#/ ) {    # detectada mala calidad en la cuarta línea
                    splice @lineas, $i, 4;            # rebanamos las 4 líneas
                    redo;              # reinicimos el bucle, en la posición $i
                }
            }
        }
 
        untie @lineas;
 
        print "$file untied\n";
    }
}
 
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

por **explorer** » 2012-06-11 16:31 @730

Te sobran las comillas simples que rodean a $file, en la línea 21. Quítalas.

Con ellas, estás diciendo que se abra el archivo '$file' (es decir: el signo del dolar, más "file").

por **Alfumao** » 2012-06-12 02:24 @141

Ya lo había probado con las comillas simples... pero no funcionaba tampoco.

He de comentar que son archivos de extensión .fastq (como se ve en la expresión regular de la línea 17) y que su tamaño es de más de 1Gb.

El resultado que debería obtener, ¿no serían los archivos originales modificados?
Porque me aparece un archivo llamado $file y que tiene 0 kb... pero no modifica los originales.

por **explorer** » 2012-06-12 02:38 @151

Claro, te ha generado el archivo $file porque usaste comillas simples...

Tie::File intentará leer todo el archivo a memoria. Quizás el problema sea el tamaño del archivo, pero 1Gb es un tamaño asequible para la memoria de un ordenador moderno.

Una comprobación sería poner una línea 22 que informe del número de líneas que se han leído.

Te sobran las paréntesis de la expresión regular.

Otra forma de hacerte con todos los archivos sería esto:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

my $path = shift @ARGV;
print "Ruta: [$path]\n";
chdir $path or die "ERROR: Unable to enter $path: $!\n";
 
for my $file (glob("*.fastq")) {
 
        print "This is $file processing...\n";
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

por **Alfumao** » 2012-06-12 03:28 @186

Hola de nuevo, explorer.
Ya ni me acordaba del comando glob...

Lo de los paréntesis en la expresión regular, supongo que te refieres a los de (\w+), los necesito porque en algún caso pretenderé obtener el nombre del archivo fastq y así paso el contenido de los paréntesis a una variable mediante $1, por eso los pongo.

No tengo muy claro cómo hacer lo de la cuenta de las líneas procesadas que me comentas...

por **explorer** » 2012-06-12 04:04 @211

La idea es sacar el número de líneas para que estemos seguros de que se ha leído correctamente el archivo:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

        print "Líneas: ", scalar(@lineas), "\n";
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

por **Alfumao** » 2012-06-12 04:50 @243

Ok, explorer, hecho y confirmado: no lee los ficheros. No sé si por la extensión (.fastq) o por el tamaño tan grande...

por **explorer** » 2012-06-12 06:12 @300

No, la extensión no es problema.

A mí sí que me funciona, pero no he probado con archivos tan grandes.

Si el mensaje que sale de print "This is $file processing...\n"; es correcto, y tu ves que salen bien los nombres de los archivos, entonces hay que pensar que Tie::File no es capaz de leer archivos tan grandes.

Haz una prueba con un archivo de unas pocas líneas, para ver si filtra las malas entradas

Yo probé con este:

Sintáxis: [ Descargar ] [ Ocultar ]

Using text Syntax Highlighting

@ILLUMINA-GA:48:FC64L4UAAXX:1:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G

@ILLUMINA-GA:48:FC64L4UAAXX:2:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G

@ILLUMINA-GA:48:FC64L4UAAXX:3:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G

@ILLUMINA-GA:48:FC64L4UAAXX:4:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G

@ILLUMINA-GA:48:FC64L4UAAXX:1:1:10257:1000 1:N:0:CGATGT

NAGTTATAGGCAGAAAGTATTCCCAGCACCTTAACCCTTACACAAGGAAAAGTTCTTTCTTTAAAAACACATTTATTTGTGTGTGTGTGT

+

#(.((''-+,:1:::888858889379898@@@@@@@@@@@@@5@5775788585850080001770500@@@@@###############

@ILLUMINA-GA:48:FC64L4UAAXX:5:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@G

@ILLUMINA-GA:48:FC64L4UAAXX:2:1:10257:1000 1:N:0:CGATGT

NAGTTATAGGCAGAAAGTATTCCCAGCACCTTAACCCTTACACAAGGAAAAGTTCTTTCTTTAAAAACACATTTATTTGTGTGTGTGTGT

+

#(.((''-+,:1:::888858889379898@@@@@@@@@@@@@5@5775788585850080001770500@@@@@###############

@ILLUMINA-GA:48:FC64L4UAAXX:6:1:4045:7422 1:N:0:CGATGT

CGAGGAAAACTGACAAAGGTGGAACATTTAGAAATATCCACTGTAGGTAGTGGAATATGGCAAGAAAACTGAAAATCATGGAAAATGAGA

+

IIIIIHIIIIIIHIIIGIIIIIIIIHIIIHBEGGGFIHIIHIIIFEIGIIGGGGEEBGDGIGBEIGGIIIHGIIDIIIFII@GIDGG@@GColoreado en 0.000 segundos,  usando GeSHi 1.0.8.4

Si te funciona, pues entonces el problema es el tamaño... así que habrá que pensar en otra solución.

Foro - Perl en Español

Eliminar secuencias de un archivo FASTQ

Eliminar secuencias de un archivo FASTQ

Publicidad

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

Re: Eliminar secuencias de un archivo FASTQ

¿Quién está conectado?