clustal one-liner con parallel

Hola,
hoy comparto un comando que a veces utilizo cuando necesito calcular muchos alineamientos múltiples a partir de una colección de ficheros de secuencias en formato FASTA. Como mi máquina, igual que la de casi todos, tiene amplia RAM y muchos cores, es un trabajo ideal para parallel. Supongamos que los archivos de salida están en la carpeta 'entrada' y queremos guardar los ficheros de salida en la carpeta 'path/to/salida', y que tenemos 20 cores disponibles:

$ mkdir /path/to/salida/

$ cd entrada

$ ls -1 *fasta | parallel --gnu -j 20 ~/soft/clustal-omega-1.2.1/src/clustalo \

--threads=1 -i {} -o /path/to/salida/{} :::

Este comando pondrá a trabajar 20 cores del sistema hasta que todos los archivos FASTA de la carpeta entrada estén alineados, con ganancias de tiempo de ejecución importantes en un experimento con 100 ficheros:

| cores (-j) | time(real) | time(user) | time(sys) |
|   1        | 4m34.440s  | 4m5.180s   | 0m2.168s  |
|  10        | 0m29.358s  | 3m57.768s  | 0m2.400s  |
|  20        | 0m23.248s  | 5m6.204s   | 0m3.364s  |

Un saludo,
Bruno

C++ STL en paralelo

Hola,
como contábamos en una entrada anterior, últimamente hemos estado liados trabajando con archivos FASTQ de varios GBs, con decenas de millones de lecturas o reads. Para algunas de las tareas que hemos tenido que realizar, como ordenar o renombrar, nos hemos dado cuenta que lenguajes interpretados como Perl o Python son varias veces más lentos que caballos de carreras compilados como C/C++, como también comenta el autor de readfq, por ejemplo.
Por eso he estado refrescando la Standard Template Library (STL) de C++, que para alguien acostumbrado a programar en lenguajes de alto nivel es esencial, como recurso para crear estructuras de datos flexibles y eficientes en programas escritos en C/C++. Como dice Scott Meyers en su libro Effective STL, probablemente los contenedores más populares de la STL son vectores y strings. Sólo por estas dos estructuras dinámicas, a las que yo añado también las tablas asociativas (maps), vale la pena aprender a usar la STL. Sin embargo, hay mucho más que esto, y por eso escribo esta entrada, porque la STL implementada en el compilador GNU gcc (libstdc++) incluye una biblioteca de algoritmos en paralelo de propósito general, que podemos usar fácilmente en nuestros programas para sacarle el jugo a los procesadores multicore y resolver de manera más eficiente múltiples problemas. Entre otros, la versión actual de libstdc++ paraleliza los siguientes algoritmos de la STL, todos ellos útiles para tareas habituales en programación:

• std::count
• std::find
• std::search
• std::replace
• std::max_element
• std::merge
• std::min_element
• std::nth_element
• std::set_union
• std::set_intersection
• std::sort

Esquema de sort parelelo y merge posterior,
tomado de http://javaero.org/tag/parallel-merge-sort

El siguiente ejemplo de sort en paralelo, que depende de la librería OpenMP, se compila con $ g++ -O3 -fopenmp -o testP testP.cc para generar un ejecutable que empleará 4 cores para ordenar un millón de palabras de 25 caracteres de longitud:

 
#include <stdlib.h>  
#include <omp.h>  
#include <vector>  
#include <parallel/algorithm>  
using namespace std;  

// g++ -O3 -fopenmp -o testP testP.cc  
// http://gcc.gnu.org/onlinedocs/libstdc++/manual/parallel_mode.html  
#define MAXTHREADS 4;  
 
string randomStrGen(int length) 
{
   //http://stackoverflow.com/questions/2146792/how-do-you-generate-random-strings-in-c    
   static string charset = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:|.@";  
   string result;  
   result.resize(length);  
   for (int i = 0; i < length; i++) result[i] = charset[rand() % charset.length()];  
   return result;  
}  

int main()  
{  
   unsigned threads = MAXTHREADS;  
   omp_set_dynamic(false);  
   omp_set_num_threads(threads);  
   srand(time(NULL));  
   std::vector<string> data(1000000);  
   for(int i=0;i<data.size();i++) data[i] = randomStrGen(25);  
   __gnu_parallel::sort(data.begin(), data.end()); //std::less<string>() , std::smaller<string>() );  
   printf("%s\n%s\n%s\n%s\n",data[0].c_str(),data[1].c_str(),data[2].c_str(),data[data.size()-1].c_str());  
   return 0;  
}  

Referencias más avanzadas:
[1] http://algo2.iti.kit.edu/singler/mcstl/parallelmode_se.pdf

[2] http://ls11-www.cs.uni-dortmund.de/people/gutweng/AD08/VO11_parallel_mode_overview.pdf

Hasta luego,
Bruno

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