Foro - Perl en Español

por **Kryban** » 2017-02-25 00:56 @080

Saludos, comunidad.

Mi problema es el siguiente, y si alguien me puede ayudar le estaré eternamente agradecido.

Tengo un archivo que contiene secuencias de 40 nucleótidos (alrededor de 2000 secuencias) y lo que debo hacer es saber si estos fragmentos de secuencias se encuentran en 6 archivos del formato fastq que traen cerca de 90 millones de secuencias.

En primera instancia, el script solicita el archivo de entrada, el nombre de salida y la dirección de la carpeta donde se encuentran los fastq para crear una lista y, de forma automática, ir cargando los demás archivos fastq cuando se termina de trabajar con uno.

Luego, a cada secuencia se busca si tiene los dinucleótidos GT (en la parte izquierda) o AG (parte derecha) y, si tiene la presencia de este patrón, se procede a buscar la secuencia que forman los 20 nucleótidos finales e iniciales de cada secuencia (izquierda y derecha) y esto lo va buscando y almacenando la cantidad de aciertos en un arreglo para luego ser mostrado en el archivo de salida.

He probado también con el módulo Bio::SeqIO::fastq pero el tiempo de ejecución es lento, ya que va secuencia por secuencia y la función grep de perl me analiza el archivo completo pero igual es lento.

Quedo atento a vuestros comentarios, sugerencias y/o críticas.

De antemano, gracias.

PD: Soy un usuario promedio de Perl.

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

!/usr/bin/perl -w
 
my $input = $ARGV[0]; #archivo de entrada
my $output = $ARGV[1];#archivo de salida
my $directory = $ARGV[2]; #dar el directorio de donde se encuentran los .fastq
 
open (IN, "< $input") or die ("no such file!\n");
open (DATAS,"> $output") or die ("don't create $output!\n");
 
print DATAS "ID_LEFT\tSTART_LEFT\tEND_LEFT\tSTRAND_LEFT\tID_RIGHT\tSTART_RIGHT\tEND_RIGHT\tSTRAND_RIGHT\tDIRECTION\tLEFT\tRIGHT\tJUNCTION\t1\t2\t3\t4\t5\t6\n";
my $linea = 0;
opendir(d, "$directory") or die "Can't open $directory: $!\n";
my @flist= grep {/.fastq/} readdir(d);
closedir(d);
 
my $cantidad=@flist;
 
while (defined($eachline=<IN>)){
        print "dentro while\n";
        if($linea!=0){
                my @array=(0,0,0,0,0,0);
                
                my @line=split(/\t/,$eachline);#columna 10 y 14
                
                my $left = $line[10];
                my $right = $line[15];
                my $GT = rindex $left, "GT";
                my $AG = index $right, "AG";
 
                if($GT!=(-1) and $AG!=(-1)){
                        my @array_main;
                        push(@array_main, $line[1], $line[2], $line[3], $line[4], $line[5],$line[6],$line[7], $line[8], $line[9],$left, $right);
 
                        my $junction_point= substr ($left, -20).substr ($right, 0, 20);#concateno las secuencias        
                        
                        foreach my $f(@flist){
                                open INFILE,$f;
                                my @result= grep {/$junction_point/} <INFILE>;
 
                                my $size=@result;                       
                                my @dor=split(/_/,$f);
                                if($size>0){
                                        @array[($dor[1]-1)]=$size;
                                }
                                close INFILE;
                        }
                        
                        print DATAS "@array_main $junction_point @array\n";
                }
        }
        $linea=($linea+1);
}
close DATAS;
close IN;
Coloreado en 0.004 segundos,  usando GeSHi 1.0.8.4

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

Hay unas limitaciones importantes: los datos de las columnas etiquetadas 1 a 6 y a sacar la información por filas (una por nucleótido) obliga a leer los archivos en un determinado orden: en el primer nivel, por cada secuencia de nucleótidos, y en el segundo nivel, por cada uno de los seis archivos grandes.

Eso tiene la penalización de que estamos leyendo los seis archivos unas 2000 veces. Es decir, aprox. 90E6 * 2E3 = 180E9.

Eso son muchas lecturas...

Hay otra opción: cambiando la forma de leer los datos. Si invertimos el orden indicado antes, solo leeremos los archivos una vez.

Primero leemos el archivo con las secuencias de nucleótidos, y vamos guardando en un array toda la información que luego usaremos en la salida.

Luego, hacemos el bucle por los seis archivos fastq. Por cada uno, hacemos el segundo bucle por todas las secuencias de nucleótidos, anotando los resultados en una matriz.

La matriz tiene, como filas, el índice de la secuencia; y como columna, el número de archivo fastq. Es decir, vamos rellenando la matriz por filas, y luego por columnas.

Al final, estamos haciendo 2E3 + 90E6 lecturas. Quedaría por sacar el resultado, recorriendo la matriz por filas.

El código quedaría algo así (no probado):

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#!/usr/bin/perl
use v5.14;
use strict;
use warnings;
use autodie;
 
my $input     = $ARGV[0]; # archivo de entrada
my $output    = $ARGV[1]; # archivo de salida
my $directory = $ARGV[2]; # dar el directorio de donde se encuentran los .fastq
 
# Lectura secuencias nucleótidos
open my $IN, '<', $input;
my $cabecera = <$IN>;
 
my @junctions;
 
while (my $eachline = <$IN>) {
 
    my @line = split /\t/, $eachline;
 
    my $left  = $line[10];
    my $right = $line[15];
 
    if (rindex($left, 'GT') != -1  and  index($right, 'AG') != -1) {
 
        # hemos encontrado un "GT"~"AG". Guardamos la información de la junction
        push @junctions, [ @line[1..9], $left, $right, (substr($left, -20) . substr($right, 0, 20)) ];
    }
}
 
close $IN;
 
# Listado de archivos fastq
opendir(my $DIR, $directory);
my @flist = grep { /[.]fastq/ } readdir $DIR;
closedir $DIR;
 
# Lectura de los archivos fasq
my @matriz;
 
foreach my $file (@flist){
    # columna a rellenar
    my @dor = split /_/, $file;
    my $col = $dor[1] - 1;
 
    # leemos el archivo fastq
    open my $INFILE, '<', $file;
    my @archivo_fastq = <$INFILE>;
    close $INFILE;
 
    for my $j (0 .. $#junctions) {
 
        # extraemos la junction_point guardada antes
        # (está en la última posición)
        my $junction_point = $junctions[$j][-1];
 
        # guardamos el número de coincidencias en la fila $j, columna $col
        my @result = grep {/$junction_point/} @archivo_fastq;
 
        $matriz[$j][$col] = @result;
    }
}
 
# Salida del resultado
open my $DATAS, '>', $output;
say     $DATAS  join "\t",
        qw(ID_LEFT  START_LEFT  END_LEFT  STRAND_LEFT
           ID_RIGHT START_RIGHT END_RIGHT STRAND_RIGHT
           DIRECTION
           LEFT RIGHT
           JUNCTION
        ),
        1..6
        ;
 
for my $j (0 .. $#junctions) {
 
    # sacamos la información, por líneas
    say $DATAS  join "\t", @{$junctions[$j]}, @{$matriz[$j]};
}
 
close $DATAS;
Coloreado en 0.003 segundos,  usando GeSHi 1.0.8.4

por **Kryban** » 2017-02-25 19:02 @834

¡wooow! Muchísimas gracias por la ayuda, realmente me sorprendió ver otra forma de plantearlo y de escribirlo.

Lo probaré. Repito los agradecimientos.

Foro - Perl en Español

Buscar secuencias pequeñas en un fastq

Buscar secuencias pequeñas en un fastq

Publicidad

Re: Buscar secuencias pequeñas en un fastq

Re: Buscar secuencias pequeñas en un fastq

¿Quién está conectado?