• Publicidad

Buscar secuencias pequeñas en un fastq

Perl aplicado a la bioinformática

Buscar secuencias pequeñas en un fastq

Notapor 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
  1. !/usr/bin/perl -w
  2.  
  3. my $input = $ARGV[0]; #archivo de entrada
  4. my $output = $ARGV[1];#archivo de salida
  5. my $directory = $ARGV[2]; #dar el directorio de donde se encuentran los .fastq
  6.  
  7. open (IN, "< $input") or die ("no such file!\n");
  8. open (DATAS,"> $output") or die ("don't create $output!\n");
  9.  
  10. 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";
  11. my $linea = 0;
  12. opendir(d, "$directory") or die "Can't open $directory: $!\n";
  13. my @flist= grep {/.fastq/} readdir(d);
  14. closedir(d);
  15.  
  16. my $cantidad=@flist;
  17.  
  18. while (defined($eachline=<IN>)){
  19.         print "dentro while\n";
  20.         if($linea!=0){
  21.                 my @array=(0,0,0,0,0,0);
  22.                
  23.                 my @line=split(/\t/,$eachline);#columna 10 y 14
  24.                
  25.                 my $left = $line[10];
  26.                 my $right = $line[15];
  27.                 my $GT = rindex $left, "GT";
  28.                 my $AG = index $right, "AG";
  29.  
  30.                 if($GT!=(-1) and $AG!=(-1)){
  31.                         my @array_main;
  32.                         push(@array_main, $line[1], $line[2], $line[3], $line[4], $line[5],$line[6],$line[7], $line[8], $line[9],$left, $right);
  33.  
  34.                         my $junction_point= substr ($left, -20).substr ($right, 0, 20);#concateno las secuencias       
  35.                        
  36.                         foreach my $f(@flist){
  37.                                 open INFILE,$f;
  38.                                 my @result= grep {/$junction_point/} <INFILE>;
  39.  
  40.                                 my $size=@result;                      
  41.                                 my @dor=split(/_/,$f);
  42.                                 if($size>0){
  43.                                         @array[($dor[1]-1)]=$size;
  44.                                 }
  45.                                 close INFILE;
  46.                         }
  47.                        
  48.                         print DATAS "@array_main $junction_point @array\n";
  49.                 }
  50.         }
  51.         $linea=($linea+1);
  52. }
  53. close DATAS;
  54. close IN;
Coloreado en 0.022 segundos, usando GeSHi 1.0.8.4
Kryban
Perlero nuevo
Perlero nuevo
 
Mensajes: 3
Registrado: 2017-02-24 22:45 @989

Publicidad

Re: Buscar secuencias pequeñas en un fastq

Notapor explorer » 2017-02-25 14:37 @651

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
  1. #!/usr/bin/perl
  2. use v5.14;
  3. use strict;
  4. use warnings;
  5. use autodie;
  6.  
  7. my $input     = $ARGV[0]; # archivo de entrada
  8. my $output    = $ARGV[1]; # archivo de salida
  9. my $directory = $ARGV[2]; # dar el directorio de donde se encuentran los .fastq
  10.  
  11. # Lectura secuencias nucleótidos
  12. open my $IN, '<', $input;
  13. my $cabecera = <$IN>;
  14.  
  15. my @junctions;
  16.  
  17. while (my $eachline = <$IN>) {
  18.  
  19.     my @line = split /\t/, $eachline;
  20.  
  21.     my $left  = $line[10];
  22.     my $right = $line[15];
  23.  
  24.     if (rindex($left, 'GT') != -1  and  index($right, 'AG') != -1) {
  25.  
  26.         # hemos encontrado un "GT"~"AG". Guardamos la información de la junction
  27.         push @junctions, [ @line[1..9], $left, $right, (substr($left, -20) . substr($right, 0, 20)) ];
  28.     }
  29. }
  30.  
  31. close $IN;
  32.  
  33. # Listado de archivos fastq
  34. opendir(my $DIR, $directory);
  35. my @flist = grep { /[.]fastq/ } readdir $DIR;
  36. closedir $DIR;
  37.  
  38. # Lectura de los archivos fasq
  39. my @matriz;
  40.  
  41. foreach my $file (@flist){
  42.     # columna a rellenar
  43.     my @dor = split /_/, $file;
  44.     my $col = $dor[1] - 1;
  45.  
  46.     # leemos el archivo fastq
  47.     open my $INFILE, '<', $file;
  48.     my @archivo_fastq = <$INFILE>;
  49.     close $INFILE;
  50.  
  51.     for my $j (0 .. $#junctions) {
  52.  
  53.         # extraemos la junction_point guardada antes
  54.         # (está en la última posición)
  55.         my $junction_point = $junctions[$j][-1];
  56.  
  57.         # guardamos el número de coincidencias en la fila $j, columna $col
  58.         my @result = grep {/$junction_point/} @archivo_fastq;
  59.  
  60.         $matriz[$j][$col] = @result;
  61.     }
  62. }
  63.  
  64. # Salida del resultado
  65. open my $DATAS, '>', $output;
  66. say     $DATAS  join "\t",
  67.         qw(ID_LEFT  START_LEFT  END_LEFT  STRAND_LEFT
  68.            ID_RIGHT START_RIGHT END_RIGHT STRAND_RIGHT
  69.            DIRECTION
  70.            LEFT RIGHT
  71.            JUNCTION
  72.         ),
  73.         1..6
  74.         ;
  75.  
  76. for my $j (0 .. $#junctions) {
  77.  
  78.     # sacamos la información, por líneas
  79.     say $DATAS  join "\t", @{$junctions[$j]}, @{$matriz[$j]};
  80. }
  81.  
  82. close $DATAS;
Coloreado en 0.021 segundos, usando GeSHi 1.0.8.4
JF^D Perl Programming Language
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14102
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Buscar secuencias pequeñas en un fastq

Notapor 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.
Kryban
Perlero nuevo
Perlero nuevo
 
Mensajes: 3
Registrado: 2017-02-24 22:45 @989


Volver a Bioinformática

¿Quién está conectado?

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