FASTQ sort + parallel

Buenas,
recientemente en el laboratorio hemos estado manipulando archivos de pares de secuencias (pair-end reads) en formato FASTQ. Una de las tareas habituales ha sido limpiar los archivos de secuencias de baja calidad, ya sea recortando o eliminando directamente, para luego volver a definir parejas entre las secuencias que superaron el corte. Una estrategia posible para esta tarea es simplemente linearizar las secuencias, de manera que cada una ocupe ahora una sola línea, separando con tabuladores la cabecera, la secuencia, el separador y las calidades. Por ejemplo, la siguiente secuencia:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTA...
+
!''*((((***+))%%%++)(...

quedaría así:

@SEQ_ID  GATTTGGGGTTCAAAGCAGTA...   +   !''*((((***+))%%%++)(...

Tras esta transformación ya sí es posible ordenar un archivo FASTQ con GNU sort, que viene instalado en cualquier sistema linux (y se puede instalar en Windows). GNU sort es ideal para ordenar conjuntos de datos que no caben en memoria (M), como a menudo ocurre con los archivos FASTQ, porque de manera implícita divide el problema inicial, de tamaño N, en N/M trozos que luego mezcla (merge) de manera externa.

En nuestra experiencia GNU sort es significativemente más eficiente que nuestros scripts para este tipo de problemas, puesto que ya trae de fábrica toda la lógica para partir el problema en trozos y luego mezclar las soluciones parciales. Sólo hay que tener cuidado de asignar la variable de ambiente LC_ALL, por ejemplo con:

$ export LC_ALL=POSIX

y echar a andar. Muy bien. Pero te quedas con la duda de si estás sacando el máximo partido a tu CPU multicore, podremos optimizar sort en paralelo? Y si invocamos a GNU parallel (mira el vídeo)?

Nos ponemos manos a la obra y hacemos pruebas con un archivo FASTQ linearizado real, de 576Mb:

$ ls -lh /tmp/unsortedXXpJ6CaB
-rw-------. 1 576M Jun  5 10:15 /tmp/unsortedXXpJ6CaB

Ahora lo ordenamos con GNU sort, dándole un máximo de 500Mb de área en RAM para trabajar (la M de antes):

$ time sort -k 1,1 -u -S 500M /tmp/unsortedXXpJ6CaB > /tmp/unsortedXXpJ6CaB.S
real 0m7.628s
user 0m5.143s
sys 0m2.373s

Finalmente probamos ahora con parallel, agrupando las secuencias en grupos de 100.000 elementos (ojo con esto, puedes llegar a obtener resultados parcialmente desordenados porque la segunda llamada a parallel puede recibir más argumentos de los que el shell soporta). Cambiando este valor a 10000 o a 10E6 los resultados son similares:

$ time cat /tmp/unsortedXXpJ6CaB | parallel -N 100000 --pipe --files sort -k 1,1 | 

   parallel -Xj1 sort -k 1,1 -u -m {} ';' rm {} > /tmp/unsortedXXpJ6CaB.P
real 0m15.451s
user 0m9.919s
sys 0m8.371s

Comprobamos que los resultados son idénticos:

$ diff /tmp/unsortedXXpJ6CaB.S /tmp/unsortedXXpJ6CaB.P

Conclusión de estas pruebas: no vale la pena complicarse con parallel para ordenar grandes archivos FASTQ, ya que probablemente el cuello de botella sea el merge final, y eso parece resolverlo mejor directamente GNU sort. De todos modos es posible que haya otras maneras de invocar a parallel más ventajosas
Si en vuestras pruebas obtenéis resultados distintos por favor escribid,

un saludo,

Bruno

perl + GNU Scientific Library

Hola,
leyendo el Linux Journal me he encontrado con este artículo que explica como usar la GNU Scientific Library en C, y me he acordado de que tenemos por el laboratorio, sin estrenar, la versión 1.12 del manual. Para qué vale esta librería? La respuesta corta es: para hacer cálculos de todo tipo con algoritmos eficientes que minimizan los errores de redondeo de los procesadores digitales.
A veces una librería es la mejor forma de emplear algoritmos extensivamente optimizados y que tal vez no sepamos cómo programar. A día de hoy la SGL incluye funciones en estas áreas:

Complex Numbers	Roots of Polynomials
Special Functions	Vectors and Matrices
Permutations	Sorting
BLAS Support	Linear Algebra
Eigensystems	Fast Fourier Transforms
Quadrature	Random Numbers
Quasi-Random Sequences	Random Distributions
Statistics	Histograms
N-Tuples	Monte Carlo Integration
Simulated Annealing	Differential Equations
Interpolation	Numerical Differentiation
Chebyshev Approximation	Series Acceleration
Discrete Hankel Transforms	Root-Finding
Minimization	Least-Squares Fitting
Physical Constants	IEEE Floating-Point
Discrete Wavelet Transforms	Basis splines

En algunos aspectos puede ser más limitada que las Numerical Recipes, pero al liberarse bajo GPL tenemos mayor libertad de uso, y sin pagar un duro. Además, gracias al trabajo de Jonathan Leto, disponemos en CPAN de una interfaz para Perl, llamada Math::GSL. Algunas de las aplicaciones de la librería se pueden conseguir con otros módulos de CPAN, pero a priori la GSL tiene a su favor: i) menor tiempo de ejecución y ii) mayor precisión numérica, una limitación de Perl ya discutida.

Pasemos a la práctica: cómo se instala esto? En mi Ubuntu 10.4 fue tan sencillo como decir $ sudo cpan -i Math::SGL y luego ir aceptando la instalación de sus dependencias.

Podemos ahora probar un ejemplo donde calcularemos una derivada numérica, mostrando algunas de las capacidades de la librería:

 #!/usr/bin/perl -w  
 # Ejemplo de derivada de ln(x) con GSL  
 # Adaptado de ejemplos de Jonathan Leto en:  
 # https://github.com/leto/math--gsl/blob/master/examples/deriv/basic  
   
 use strict;  
 use Math::GSL::Deriv qw/:all/;  
 use Math::GSL::Errno qw/:all/;  
   
 # incremento de derivada  
 my $h = 0.01;   
   
 # queremos X = pi, ejemplo constante matematica en GSL   
 my $x = $Math::GSL::Deriv::M_PI;   
    
 # derivada numerica, recuerda en perl log = ln  
 my ($status,$val,$err) = gsl_deriv_central ( sub { log($_[0]) }, $x, $h);  
   
 # derivada analitica, demostracion en  
 # http://www.math.com/tables/derivatives/more/es-ln.htm  
 sub dlndx{ return 1/$_[0] }  
   
 if($status == $GSL_SUCCESS)   
 {  
   printf("deriv(ln((%g)) = %.18g, error maximo esperado = %.18g\n", $x, $val, $err);  
   printf("dlndx(%g)   = %.18g\n" , $x, dlndx($x));  
   printf("error observado   = %.18g\n",abs($val-dlndx($x)));  
 }   
 else   
 {  
   my $gsl_error = gsl_strerror($status);  
   print "ERROR: $gsl_error (es derivable log(x) en ese punto?)\n";  
 }  

Se obtiene algo como:

deriv(ln((3.14159)) = 0.31830988618650752, error maximo esperado = 5.7425388097458009e-11
dlndx(3.14159)      = 0.31830988618379069
error observado     = 2.7168267635602206e-12

El autor del módulo tiene en la web una colección de ejemplos que sin duda amplian la descripción de las capacidad de GSL. Por ejemplo, en este ejemplo se explica como calcular productos escalares de dos maneras, encontrando que la implementación del módulo Math::GSL::BLAS es al menos el doble de rápida.

Ya que este ejemplo no representa las tareas típicas de la programación en Bioinformática os muestro un segundo mucho más terrenal, donde ordenamos vectores de números reales, comparando la implementación de mergesort de Perl con la heapsort inestable de GGL. Hay un capítulo del manual dedicado sólo a esto, pero probablemente la documentación del módulo Math::GSL::Sort sea suficiente:

 use Math::GSL::Sort qw/:all/;  
   
 my $numbers = [ map { rand(100) } (1..100000) ];  
 my $sorted = gsl_sort( $numbers,1,scalar(@$numbers) );     
 my @sortedp = sort {$a<=>$b}( @$numbers );  

Si hacéis pruebas veréis que la implementación GSL resulta ser más rápida, en mi caso aproximadamente un 20%.

Otras aplicaciones posibles incluyen el cálculo de permutaciones y combinaciones, algo ya discutido en este blog,
hasta otra,
Bruno

PD Si alguien ha probado a instalar este módulo en sistemas Windows le agradeceremos sus comentarios

Ordenamiento de resultados de BLAST

Un problema con el que me he encontrado recientemente es el de ordenar resultados de BLAST contenidos en diferentes ficheros. Para definir mejor el problema supongamos que tenemos 4 genomas A,B,C y D y queremos comparar, todas contra todas, sus secuencias de proteínas por medio del programa blastp. Una manera de hacerlo sería hacer 4x4 comparaciones por parejas (AA,AB,AC,AD,BA,BB,BC,BD,CA,CB,CC,CD,DA,DB,DC,DD) teniendo en cuenta que la dirección en BLAST normalmente importa.

Una vez completados estos 16 trabajos de BLAST encontraremos que cada uno de los archivos de salida están internamente ordenados en términos de E-value. Sin embargo, dada la secuencia 1 del genoma A, si queremos averiguar en qué genomas se encuentra su secuencia más similar (best hit), deberemos mirar en todos los archivos (7 de 16) donde aparece el genoma A. Otra alternativa, la que da pie a este artículo, es fusionar, intercalar y ordenar (merge-sort) esos 16 archivos internamente ordenados para crear uno sólo, que facilite posteriores consultas.
Sirvan de muestra los archivos AA.blast:

16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8

y AB.blast:

16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103

que contienen, en formato tabular, los alineamientos significativos de la secuencia 16 dentro del genoma A (contra las secuencias 16,18,78) y dentro del genoma B (582,637), con valores esperados de 0.0, 2e-24,1e-19,5e-23,3e-21 respectivamente.
El problema que vamos a resolver es como combinar estos dos archivos (o en general N archivos) en uno sólo hasta obtener una lista ordenada por E-value:

16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8

Sin más preámbulos, el siguiente código Perl implementa dos posibles soluciones:

 package BlastSort;   
   
 use strict;  
 use Symbol qw(gensym);  
 use sort 'stable';   
   
 require Exporter;  
 use vars qw(@ISA @EXPORT);  
 @ISA = 'Exporter';  
 @EXPORT = qw(merge_BLAST_files merge_BLAST_files_GNUsort);  
   
 sub merge_BLAST_files   
 {  
    # Adapted from File-Sort-1.01(http://search.cpan.org/perldoc?File::Sort)  
    # Assumes infiles are BLAST output files in tabular format with sequences  
    # identified by natural numbers, such as 11,12,1439,1440 in the sample:  
    #11   1439   78.24   625   136   0   4   628   6   630   0.0    993  
    #12   1440   80.88   272   52   0   1   272   1   272   7e-125    446  
    # Order of infiles is RELEVANT as merging is stable, so sequences from  
    # the first files will be given sorting priority  
      
    my ($outfile,@infiles) = @_;  
      
    my (%fh,%order,@fhorder,$filein,$id,$first,$n_of_fhs,$curr,$line);  
   
    # for Schwartzian transform (ST) see  
    # http://www.hidemail.de/blog/perl_tutor.shtml#sort_orcish_schwartzian  
    sub blastsort { $a->[2] <=> $b->[2] || $a->[3] <=> $b->[3] }  
    sub blastmap   
    {   
       my @tmp = split(/\s+/,$fh{$_});   
       [$_,$fh{$_},$tmp[0],$tmp[10]]   
       # returns anonymous array[4]: filehandle, line,query_id,E-value  
    }  
   
    ## open all input BLAST files and keep filehandles, ordered 0..N  
    $n_of_fhs = 0;  
   foreach $filein (@infiles)   
    {  
      $id = gensym(); # get valid id for filehandle  
         
       open($id,$filein) ||   
          die "# merge_BLAST_files : cannot read $filein: $!";  
   
       $order{$n_of_fhs} = $id;  
       $n_of_fhs++;  
    }  
    @fhorder = (0 .. $n_of_fhs-1);  
      
     
    ## open outfile and ensure IO buffer is used  
    $| = 0;   
   unlink($outfile) if(-s $outfile);  
    open(OUT,">$outfile") ||   
       die "# merge_BLAST_files : cannot create $outfile: $!";  
   
   ## get first BLAST line from all filehandles    
   %fh = map {  
     my $fh = $order{$_};  
     ($_ => scalar <$fh>);  
   } @fhorder;  
   
   ## start merging BLAST lines   
    while(scalar(@fhorder)>1)   
    {        
       ($first) = (map {$_->[0]} sort blastsort map &blastmap, @fhorder); #ST  
         
       print OUT $fh{$first};  
   
     $curr = $order{$first};  
     $line = scalar <$curr>;  
     if(defined($line)) # update current filehandle  
       {   
          $fh{$first} = $line;   
       }  
       else # exhausted filehandle  
       {   
          @fhorder = grep { $_ ne $first } @fhorder;  
       }  
   }  
      
    ## take care of last filehandle left and close file  
    print OUT $fh{$fhorder[0]};  
    $curr = $order{$fhorder[0]};  
    while(<$curr>){ print OUT }  
   close(OUT);    
      
    ## close all input files  
    foreach $id (0 .. $n_of_fhs-1){ close($order{$id}); }  
      
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
   
 sub merge_BLAST_files_GNUsort  
 {  
    my ($outfile,$tmpdirpath,$maxbufferMb,@infiles) = @_;  
   
    # local sort -k11g fails with: 0.00,0.006 and 1e-23 (coreutils v6.10)  
    # probably as LC_ALL is not set to 'POSIX'  
    # http://www.gnu.org/software/coreutils/faq/coreutils-faq.html   
   # #Sort-does-not-sort-in-normal-order_0021  
    $ENV{'LC_ALL'} = 'POSIX';  
     
   unlink($outfile) if(-s $outfile);  
   
   my $sort_command = "sort --temporary-directory=$tmpdirpath " .  
     "--buffer-size=$maxbufferMb -s -k1g -k11g -m " .  
     join(' ',@infiles)." > $outfile ";  
   
   system("$sort_command");    
   
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
 __END__  

El módulo contiene dos subrutinas,  merge_BLAST_files y merge_BLAST_files_GNUsort; mientras la primera muestra explícitamente cómo se desarrolla el ordenamiento externo en disco, manteniendo los N ficheros de entrada abiertos y sin guardar nada en memoria, la segunda subrutina es realmente una llamada a la utilidad sort del shell (GNU coreutils 6.10 en mi sistema), donde sí usamos la memoria para acelerar el ordenamiento, 500Mb en el siguiente ejemplo:

 #!/usr/bin/perl  
 use strict;  
 use Benchmark;  
 use BlastSort;  
   
 my @infiles = ( 'A-A.blast','A-B.blast' );    
   
 my ($outfile,$gnu_outfile,$start_time,$end_time,$sortOK) =   
    ('out.merge-sort.txt','out.gnu-sort.txt');  
   
 print "# number of BLAST files to merge-sort = ".scalar(@infiles)."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files($outfile,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime (BlastSort): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files_GNUsort($gnu_outfile,'./',500,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime    (GNU): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  

Como se observa en la gráfica, aunque la versión Perl es autoexplicativa, es exponencialmente más lenta que GNUsort. Por tanto, en general convendrá usar la segunda opción y la primera sólo tendrá interés si queremos reservar la RAM, a costa de tiempo de ejecución. La única posible complicación del GNUsort es que depende de la variable de ambiente LC_ALL, a la que deberemos dar el valor 'POSIX' para tener resultados correctos. De no ser así no ordena bien los E-valores.