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.
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.003 segundos, usando GeSHi 1.0.8.4