Algunos entresijos de TopHat2.0.9

Todo software tiene sus pieles y sus tripas. Las pieles suelen venir en forma de parámetros que permiten modificar el comportamiento de un programa. Conocer las tripas suele referirse a sumergirse en el código fuente y conocer en profundidad los detalles del algoritmo en cuestión, las estructuras de datos que utiliza y cómo trabaja con la información.

Sin embargo, en algunos casos es posible quedarse digamos en el tejido conectivo, una mezcla de información proveniente de documentación, listas de correo, foros, logs del programa, diversos tests y publicaciones que explican ciertos entresijos del programa sin llegar al detalle del código. Así avanzaremos un poco en la compresión de su funcionamiento y, lo que es más importante, de los resultados obtenidos. Sin embargo, puede llevar, como veremos, a nuevas preguntas sin resolver, la identificación de posibles bugs y a ser muy cauto cuando uno tiene que interpretar los resultados de un programa que no conoce en detalle. Sin embargo, lo positivo es que podremos estar al tanto de dichas peculiaridades, e incluso podremos colaborar a mejorar los programas indirectamente informando a los meritorios desarrolladores de los resultados de nuestras pruebas. Podremos también ser más críticos ante los resultados generados con uno de estos programas.

En el caso de esta entrada, vamos a hablar un poco de TopHat2.0.9, que es el programa más utilizado para mapear reads cortos de secuenciación masiva de cDNA a una referencia, normalmente un genoma ensamblado. Es una entrada bastante prescindible para quien no vaya a trabajar con TopHat2 o programas de alineamiento similares, aviso.

Digamos que TopHat2 tiene 3 formas básicas de operación, según vayamos añadiendo una nueva capa de filtrado.

La fase común a toda ejecución de TopHat2 consiste en un mapeo contra el genoma. Todos aquellos reads no mapeados serán fragmentados, y dichos segmentos volverán a ser alineados, con el objetivo de encontrar posibles nuevas uniones entre exones, de forma que finalmente podamos mapear dichos reads en toda su longitud.

Sin embargo, incluyendo la anotación de un transcriptoma, en formato GFF o GTF, TopHat2 mapeará antes de nada nuestros reads sobre las regiones especificadas, de forma que los reads que alineen en el transcriptoma los va a considerar directamente como mapeos buenos y dejará de buscar alineamientos para esos reads ("Only the reads that do not fully map to the transcriptome will then be mapped on the genome", http://tophat.cbcb.umd.edu/manual.shtml). Posteriormente, todo lo que no alineó contra el transcriptoma pasa por la fase descrita en el párrafo anterior: mapeo contra el genoma, realineamiento de segmentos de los reads sin hits, definición de nuevas uniones entre exones, alineamiento final de reads.

Por último, para ciertas aplicaciones, como variant calling, ensamblaje o caracterización estructural del transcriptoma; también según con la naturaleza de los datos y referencias con los que se trabaje; puede convenir filtrar aquellos reads con mapeos múltiples. Especialmente cuando se trabaja con referencias ricas en repeticiones, como es el caso de los genomas de algunas plantas, por ejemplo. Con TopHat2 podemos utilizar la opción -M (o --prefilter-multihits) para realizar un prefiltrado mapeando reads contra el genoma. En esta fase, todos aquellos reads que superen el número de alineamientos indicado por el parámetro -g serán descartados. Así, las fases de alineamiento al transcriptoma y subsiguientes, se realizarán únicamente con los reads con el número máximo de alineamientos que deseemos. Obviamente, si queremos reads que alineen una única vez al genoma deberíamos utilizar la opción "-g 1".

Éste es un pequeño pero importante detalle del funcionamiento de TopHat2 cuando se quiere utilizar la opción -M. (Ojo que vienen curvas a partir de aquí :) ). Hay que ser conscientes de que el parámetro -g (o --max-multihits) se comporta como un threshold en la fase de prefiltrado, de forma que los reads que superen dicho threshold serán descartados ("[...] directs TopHat to first align the reads to the whole genome in order to determine and exclude such multi-mapped reads (according to the value of the -g option)", http://tophat.cbcb.umd.edu/manual.shtml). Sin embargo, el parámetro -g en las fases subsiguientes sirve simplemente para indicar a TopHat2 cuántos alineamientos como máximo queremos que aparezcan en el BAM para cada read. Si hay más alineamientos que el valor de -g, simplemente se elige cuáles incluir en la salida aleatoriamente ("[...] will report the alignments with the best alignment score. If there are more alignments with the same score than this number, TopHat will randomly report only this many alignments", http://tophat.cbcb.umd.edu/manual.shtml).

Esto trae efectos secundarios, quizás inesperados. Sólo podemos definir -g una vez, y será utilizado en todas las fases. En el prefiltrado, utilizando "-g 1", eliminaremos todos los reads con alineamientos múltiples, al genoma. Pero no mapeos múltiples generados en fases posteriores*. Si queremos detectar dichos mapeos múltiples no podremos en principio, porque con la opción "-g 1" se nos informa de un sólo alineamiento por read, y además los valores MAPQ y NH:i:1 del BAM indican todos y cada uno de los reads como de alineamiento simple (!). Esto implicaría que no podríamos diferenciar dichos reads tampoco en programas que utilicen el BAM como entrada (como samtools, por ejemplo).

Hablando de NH:i, esta es otra peculiaridad de TopHat2 que hay que conocer. A pesar de que utiliza Bowtie2 como algoritmo de alineamiento, de alguna forma se ha descartado el uso de la tag XS:i (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#bowtie2-build-opt-fields-xs) y se ha optado por utilizar NH:i. Sin entrar en las tripas del programa no sabemos exactamente en qué momento se realiza dicha sustitución. Ni por supuesto el porqué o para qué. Esto no está recogido en la documentación de TopHat2.

Volviendo al prefiltrado: ¿por qué no obviar el uso de la opción -M y directamente filtrar los mapeos múltiples marcados tras el mapeo al transcriptoma y al genoma? Utilizando la opción del prefiltrado obtendremos muchos más mapeos múltiples, ya que no se da por sentado que el mapeo al transcriptoma es el correcto para un read en particular (como se ha explicado más arriba, al comentar la fase de mapeo al transcriptoma). Véase la comparación de un mapeo con y sin -M:

	contigs mapped	Reads mapped
tophat_s_04	508,193	28,660,380	87.95%
tophat_s_04_M_g2	41,321	16,375,099	50.25%

La diferencia de las dos filas anteriores se debe únicamente al uso de "-M -g 2" (activos en la segunda fila). Como puede verse, el resultado final es muy diferente. Esta diferencia podría ser mayor incluso con la opción "-g 1". La disminución en el número de reads mapeados es importante, pero cómo se refleja esto en la diversidad de hits sobre la referencia es abismal. En el segundo caso estamos mapeando con una confianza mucho mayor a una región más reducida del genoma, quizás indicando que está compuesta por zonas de alta mapabilidad.

Ahí ya depende del criterio de cada uno, del objetivo del análisis, de la mapabilidad del genoma y la completitud de su anotación. Quizás para genomas con un transcriptoma muy bien anotado, en el que todo lo que quede fuera puedan ser pseudogenes o repeticiones que se sabe no son codificantes y si se están análizando únicamente genes que codifiquen para proteínas, por poner un ejemplo, sea más conveniente no usar el prefiltrado. Quizás para genomas sin acabar, con regiones de transcriptoma no del todo delimitadas sea más seguro utilizar la opción -M, especialmente si se están tratando de detectar loci que puedan ser regiones transcritas no identificadas previamente, como posibles parálogos.

Otra cosa interesante que me ha quedado por resolver respecto a los mapeos múltiples es en referencia a la nueva salida que genera TopHat2.0.9: align_summary.txt. Utilizando "-g > 1" los reads con alineamientos múltiples aparecen con MAPQ < 50 y NH:i>1 en el BAM. Sin embargo, no he conseguido ningún alineamiento en el que el número de reads así indicado coincida con el número de reads al que se refiere en el fichero align_summary.txt.

Respecto al pre-filtrado y su compatibilidad con la detección de mapeos múltiples en fases posteriores (debido a lo comentado antes), podríamos de nuevo cuestionarnos ¿Por qué no realizar un mapeo inicial mediante Bowtie2 para detectar dichos mapeos y luego ya utilizar TopHat2 de forma habitual? Esto puede realizarse, sin duda. Pero una de las características más interesantes de TopHat2 es el control que tenemos sobre el número de mismatches y bases en gaps, así como la edit distance final en el alineamiento (como vamos a ir viendo). Estos parámetros además son comunes para todas las fases de alineamiento, lo que nos permite ser consistentes en cierta medida en cuanto a qué consideramos un alineamiento válido para nuestras muestras. Veamos un poco más en detalle qué pasos realiza TopHat2 para los alineamientos, gracias a la información presente en logs/run.log.

Para cada alineamiento, TopHat2 ejecuta Bowtie2 mediante un wrapper (bowtie2_align). En principio, TopHat2 realiza los alineamientos con los reads por separado, sean reads pareados o no ("For paired-end reads, TopHat2 processes the two reads separately through the same mapping stages described above. In the final stage, the independently aligned reads are analyzed together to produce paired alignments, taking into consideration additional factors including fragment length and orientation", http://genomebiology.com/2013/14/4/R36). Para la ejecución de Bowtie2 podemos pasarle ciertos parámetros que serán utilizados siempre por el wrapper, como los presets de sensibilidad (que incluyen -D, -R, -N, -L, -i; todos parámetros relacionados con las seeds del alineamiento de Bowtie2) o los de scoring (mismatches, gaps, etc.). También podemos incluir --b2-score-min, que hará las veces de filtro de alineamientos que no superen dicho score.

Es importante diferenciar el uso que hace el wrapper de -k en las distintas fases de alineamiento, porque éste es un parámetro clave de la configuración de Bowtie2 que no podemos cambiar en TopHat2, al menos no de forma directa. Parece, junto a la actividad de fix-map-ordering (que vemos un poco más abajo), uno de los pilares fundamentales del programa, de lo que lo hace un mapeador especializado en datos de RNAseq.

Para la fase de alineamiento al transcriptoma se utiliza -k 60, mientras que para el alineamiento posterior al genoma se utiliza -k 20. Para los segmentos se pasa a -k 41, al igual que para el mapeo de los reads a las nuevas uniones de exones. Curiosamente, en el paso de prefiltrado (-M) se utiliza -k 2. ¿Qué demonios es -k? Pues el -k de Bowtie2 ("[...] it searches for at most distinct, valid alignments for each read", http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#bowtie2-options-k). 

En principio tengo curiosidad por ver qué pasaría si al prefiltrado le indicáramos "-g 3". Pero no lo he probado, ni tiene pinta de que lo haré. Por otro lado, digamos que la búsqueda más amplia que va a realizar TopHat2 para encontrar el mejor alineamiento posible va a ser sobre el transcriptoma, aprovechando además que es un conjunto de secuencias más reducido que el genoma. Unos cuantos alineamientos menos realizará TopHat2 con los segmentos y los realimientos a uniones de exones. Respecto al alineamiento al genoma, si vamos a jugar con un genoma altamente repetitivo, qué garantías nos da -k 20? ¿Habría que modificarlo al trabajar con un genoma de bacteria, o con el genoma del pino? Preguntas sin resolver ¿tú qué opinas?

Siguiendo con los parámetros que le pasamos al wrapper, es importante no ser demasiado estricto en cuanto al score mínimo permitido (--b2-score-min), o incluso el scoring de los mismatches y gaps, si se plantea uno el alinear reads con 3 mismatches (http://genomebiology.com/2013/14/4/R36, Table S2).

El siguiente paso que TopHat2 realiza sobre cada alineamiento es ejecutar el script fix_map_ordering, al que le pasa los parámetros --read-mismatches, --read-gap-length, --read-edit-dist, --read-realign-edit-dist, --bowtie2-min-score. Éste último en principio no lo podemos cambiar. Con lo cual parece un threshold fijo, aunque no me ha quedado claro aún que relación tiene con el score-min de Bowtie2 que podemos cambiar para el wrapper bowtie2_align. Este --bowtie2-min-score cambia en las distintas fases, siendo de 15 para el transcriptoma y el genoma, de 10 para los segmentos y realineamientos a junctions. Su valor es 35 para el alineamiento al genoma de la opción -M. El resto de parámetros sí los podemos cambiar, y esto puede ser realmente útil si queremos controlar el número de mismatches y bases en gaps que queremos aceptar. Cuando cambiemos estos parámetros, serán utilizados en todas las fases de alineamiento, eso sí.

                     Reads mapped
tophat_s_a 12,806,680 39.30% --read-mismatches 0 --read-gap-length 0 --read-edit-dist 0
tophat_s_b 12,789,305 39.25% --read-mismatches 0 --read-gap-length 0 --read-edit-dist 6
tophat_s_c 15,038,778 46.15% --read-mismatches 1 --read-gap-length 0 --read-edit-dist 6
tophat_s_d 15,056,015 46.20% --read-mismatches 1 --read-gap-length 0 --read-edit-dist 1
tophat_s_e 15,539,057 47.69% --read-mismatches 1 --read-gap-length 3 --read-edit-dist 6
tophat_s_f 16,388,235 50.29% --read-mismatches 2 --read-gap-length 3 --read-edit-dist 6
tophat_s_g 16,478,122 50.57% --read-mismatches 2 --read-gap-length 6 --read-edit-dist 6
tophat_s_h 17,038,705 52.29% --read-mismatches 4 --read-gap-length 6 --read-edit-dist 6
tophat_s_ia 17,046,981 52.31% --read-mismatches 4 --read-gap-length 6 --read-edit-dist 12
tophat_s_ib 17,053,281 52.33% --read-mismatches 4 --read-gap-length 12 --read-edit-dist 12
tophat_s_l  17,188,433 52.75% --read-mismatches 8 --read-gap-length 12 --read-edit-dist 18

--b2-very-sensitive --b2-score-min C,-28,0 -M -g 1 (Parámetros usados en todos los casos)

Esta es una tabla con alineamientos cambiando los parámetros de fix-map-ordering. En éste caso el cambio en el número de mismatches aceptados tiene un impacto mayor en el resultado que el aumento del número de bases en gaps. Especialmente aumenta el número de reads mapeados en el cambio a 1 mismatch, 2 mismatches, 4 mismatches... Ya no aumenta tanto al utilizar 8 mismatches, con 12 bases en gaps, y apenas unos 20,000 alineamientos de los más de 100,000 reads mapeados nuevos, eran de mapeo único (datos no mostrados).

Para el caso de tophat_s_h, se generaron alineamientos modificando --b2-score-min. Como se observa, el número de reads mapeados contra el genoma parece saturarse ya con score entre 8 y 16.

                             mapped

0 13,848,299 42.50%

4 14,205,101 43.59%

8 16,350,857 50.18%

16 17,118,163 52.53%

32 17,486,162 53.66%

64 17,351,895 53.25%

¿Es esta una saturación real, debido a que los reads mapeables independientemente son cerca de los 17 M en esta muestra? O se debe quizás a que está actuando el score mínimo de fix_map_ordering (--bowtie2-min-score 15). En la tabla anterior vemos los alineamientos del prefiltrado al genoma. Sin embargo, a medida que se reduce el número de reads sin mapear, aumenta el número de reads con mapeos múltiples (datos no mostrados).

De esta forma, el número reportado con TopHat2 es bastante robusto frente a --b2-score-min usando la opción -M (con el genoma de éste ejemplo, etcétera.), como se ve en la tabla siguiente (--b2-score-min de 0 a 64).

0	16,796,134	51.54%
4	16,756,691	51.42%
8	16,573,331	50.86%
16	16,649,626	51.09%
32	17,091,930	52.45%
64	17,125,504	52.55%

Otras opciones importantes, en éste caso para el uso con reads pareados, son --no-discordant y --no-mixed (http://tophat.cbcb.umd.edu/manual.shtml). Curiosamente, estas opciones no aparecen en logs/run.log, por lo que no queda claro cómo y en qué momento las usa TopHat2. En principio parecería que lo lógico es que sean las opciones de Bowtie2 (del mismo nombre), pero por otro lado esto parece poco probable si realmente TopHat2 mapea cada pareja por su lado y finalmente acaba realizando un chequeo de las parejas (momento en el cual podría ser que se apliquen estos parámetros). En cualquier caso, su reflejo en los resultados es importante:

            contigs  Reads unmapped filtered multiple Reads mapped
M_1_b 54,537 17,447,158 19.62% 36,442,332 40.99% 35,017,524 39.39%
M_1_c 31,550 24,460,620 27.51% 39,179,540 44.07% 25,266,854 28.42% --no-discordant --no-mixed

Como vemos afecta al número de contigs mapeados, así como al número de reads que de los que no se ha conseguido un alineamiento válido y en pareja. Esto tiene su efecto en el número de reads mapeados final. Pero deberían ser estos reads los que más confianza nos den en el caso de referencias de baja mapabilidad, ya que son alineamientos únicos (en su mayoría), con la pareja alineando al mismo sitio y las distancias y sentidos de los miembros del par dentro de lo esperado. Sorprendentemente, aumenta también el número de reads prefiltrados, hecho que hemos visto consistentemente con nuestros datos. Si --no-discordant y --no-mixed son opciones primordialmente restrictivas ¿qué puede estar pasando realmente para que aumente el número de reads con alineamientos múltiples al genoma? ¿Sugerencias?

Hasta aquí un breve vistazo a partes del tejido conjuntivo de TopHat2.0.9. Hemos visto que podemos utilizar básicamente 3 modos de operación, y que su uso tiene un reflejo importante en los resultados obtenidos. También hemos resuelto algunas dudas respecto al funcionamiento de parámetros como -g y -M, o al uso de las tags NH:i con "-g 1" (probablemente erróneo). Hemos averiguado también que hay ciertos parámetros que no podemos modificar directamente (como -k o --bowtie2-min-score), y parecen ser características intrínsecas de la naturaleza de TopHat2 en unos casos, auténticas incógnitas en otros. Finalmente, se han generado algunas dudas sin resolver en la interpretación de los resultados, como por qué aumenta el número de reads con alineamiento múltiple con la opción --no-mixed, o qué demonios intenta decirnos el fichero align_summary.txt acerca de nuestros reads y su singularidad.

En fin, si alguien se anima a meterse en las tripas y hacer una segunda parte a éste post, sería fantástico. Por el momento, a ver si en la próxima entrada hago algo un poco más light y os cuento un poco sobre www.circos.ca.

Hasta la vista!

* No está muy claro de dónde proceden dichos mapeos múltiples, pero aparecen. Quizás sea el resultado de varios efectos. Por ejemplo, el mapeo a las nuevas junctions generadas durante la fase de re-alineamiento de segmentos; o que no todo lo que mapeaba de forma múltiple fue encontrado en la fase de prefiltrado. La significación cuantitativa tampoco queda del todo clara. Un ejemplo de los que tengo, de ~32,5 M de reads se prefiltraron ~ 9.7 M. Además, ~ 5.8 M no alinearon al genoma en esta primera fase por lo que fueron descartados también. El resto de reads se volvieron a pasar por TopHat2, ahora sin -M y con -g 2. De unos 17 M de reads, 0,8 M (~ 4,7% en el peor de los casos) tenían mapeos múltiples según align_summary.txt; unos 31 K (~ 0,18% en el mejor) según NH:i y MAPQ.

contrato para tesis "Genome annotation, comparative genomics and phylogenomics of the Brachypodium model complex"

El Laboratorio de Biología Computacional (EEAD-CSIC), en colaboración con el grupo de la Dra. Pilar Catalán de la Escuela Politécnica Superior de Huesca (EPS, U.Zaragoza),  ofrece un contrato de 4 años para la realización de la tesis doctoral "Genome annotation, comparative genomics and phylogenomics of the Brachypodium model complex (Poaceae)", ligado al proyecto CGL2012-39953-C02-01. Buscamos preferentemente a personas con formación universitaria en ciencias biológicas, ambientales o agrarias. Valoraremos especialmente la experiencia en áreas como bioinformática, biología molecular y/o botánica.

Se plantea realizar el trabajo en la EPS, con sede en Huesca, con alguna estancia temporal en la EEAD en Zaragoza. Los directores del trabajo serán la propia Dra. Pilar Catalán (biología molecular, filogenómica y evolución de Brachypodium sp.) y el Dr. Bruno Contreras (EEAD-Zaragoza: bioinformática, ensamblaje y análisis genómicos).

 
Referencias relevantes de trabajo previo:
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0051058
http://aob.oxfordjournals.org/content/early/2012/01/01/aob.mcr294.full
http://www.brachy2013.unimore.it/images/stories/file/Brachy2013_Abstract_book.pdf


Se ruega a las personas interesadas se pongan en contacto con los directores del trabajo (bcontreras@eead.csic.es, pcatalan@unizar.es) para enviar una carta de interés y una copia del CV. El plazo para la presentación de solicitudes comprende del 26 de  agosto de 2013 al 10 de septiembre de 2013.

El texto completo de la convocatoria se encuentra en:
http://www.boe.es/diario_boe/txt.php?id=BOE-A-2013-8984 
http://tinyurl.com/k2eymfu


Requisitos:
1. Podrán ser solicitantes todas aquellas personas que se encuentren en disposición de estar matriculado o admitido en un programa de doctorado, en el momento de la formalización del contrato. 
2. Con carácter general, los solicitantes deberán haber finalizado sus estudios, considerándose como fecha de finalización aquella en la que se acredite que se han superado todas las materias y requisitos académicos que dan acceso a un programa dedoctorado, en fecha igual o posterior al 1 de enero de 2010.
3. No podrán ser solicitantes quienes ya estén en posesión del título de Doctor, por cualquier universidad española o extranjera.

English abstract: 
The Laboratory of Computational and Structural Biology (EEAD-CSIC), in collaboration with the group of Dr.Pilar Catalán (EPS,U.Zaragoza),  is currently offering a 4-year PhD position entitled "Genome annotation, comparative genomics and phylogenomics of the Brachypodium model complex (Poaceae)". 
Note: check SCIMAGO (http://www.scimagoir.com/pdf/SIR%20Global%202013%20O.pdf)
for the impact of CSIC among global research institutions.

Comparando proveedores de secuenciación

Hola,
para romper con mis últimas entradas, donde hablaba de código fuente, hoy quisiera discutir un tema recurrente cuando planeamos experimentos de secuenciación masiva.

Cuál es el problema? Pues que a menudo se nos ocurren experimentos pero no sabemos dónde podemos llevarlos a cabo, porque no tenemos en casa un HiSeq, ni cuánto nos costarán. El mapamundi NGS, del que hablaba Carlos P Cantalapiedra en 2011, y que ahora vive en http://omicsmaps.com, nos puede ayudar a buscar proveedores cercanos, lo cual puede ser importante para no mandar muestras difíciles de obtener y facilísimas de estropearse, a la otra punta del munco.

Sin embargo, otra cuestión es el precio y el plazo de entrega, y si ofrecen apoyo bioinformático básico. Nosotros mismos hemos tenido que tirar de teléfono y llamar a proveedores variopintos para tener unos cuantos presupuestos en la mano y decidirnos por el mejor. Hasta la fecha, dentro de España hemos contratado con el Parque Científico de Madrid y con el CNAG. Pero también hemos contactado con proveedores europeos y norteamericanos, y aquí es donde viene a cuento https://genohub.com, al que no quiero dar publicidad gratuita, pero es que es la única herramienta de este tipo que conozco, y que puede ahorrarnos mucho esfuerzo en la elaboración de presupuestos y la preparación de experimentos. Si conocéis otros recursos similares por favor ponedlos como comentarios,
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

Potenciando el entorno Python I: webapps con CherryPy

Buenas!

Todos necesitamos utilizar distintos lenguajes de programación y distintas sintaxis y convenciones para muy diversas tareas. A saber: bash scripting, sed, awk, C, Perl, Python, PHP, Java, XML, HTML, JavaScript, JSON, AJAX, etcetc

Yo soy partidario de tratar de evitar esto lo máximo posible. Así que últimamente estoy tratando de potenciar mi entorno de trabajo Python y moverme a otros lenguajes sólo cuando sea realmente necesario. He hecho algunas nuevas adquisiciones a mi repertorio últimamente, aunque todavía sin mucho dominio de ellas. A pesar de esto, los resultados me han parecido muy satisfactorios. Aquí voy a hablar de un Framework para aplicaciones web en Python, CherryPy (a minimalist python web framework).

Cuales son las claves para que me haya decidido por utilizar CherryPy:
- Python no incorpora las utilidades de aplicación web que vienen incluídas en, por ejemplo, PHP. En mi caso, especialmente, control de sesión HTTP.
- Necesitaba trabajar tras un servidor Apache y tanto el módulo CGI de Python, como el mod_python de Apache son desaconsejados.
- Una aplicación CherryPy es un servidor HTTP en sí mismo, por lo que puede trabajarse muy cómodamente sobre la aplicación sin necesidad de involucrar otros sistemas en un principio.
- Permite trabajar única y exclusivamente en Python. Incluso la configuración se puede realizar con diccionarios de Python y decorators, si se prefiere esto a tener ficheros de configuración externos. La petición HTTP (la request) y sus parámetros se manejan directamente desde métodos Python.
- Es una Framework ligera, que no se inmiscuye en creación de vistas, control de usuarios, etc. Si bien permite la creación y utilización de herramientas y plugins complementarios. Siempre me han atraido las pequeñas cosas modulares, por encima de las grandes aplicaciones todo en uno. Las primeras me sugieren versatilidad y claridad; las segundas dependencia.

Quizás la parte más complicada para tirarse a trabajar con CherryPy es la configuración. Aunque más que por su complejidad, es debido a la errática documentación existente y los ligeros cambios entre versiones de la framework que para un novicio pueden ser difíciles de atribuir a una u otra versión*. En mi caso, lo mejor ha sido obviar a menudo la documentación de la página principal y trabajar directamente con el libro CherryPy Essentials y apoyándome en la lista de usuarios de CherryPy.

Sin embargo, ya digo que no es muy complicado. Una vez uno se tira a escribir sobre la aplicación y va viendo que las cosas funcionan, en seguida se pasa a trabajar con la aplicación en Python y se olvida de lo que hay detrás. Eso hasta que ejecutamos la aplicación y voilá! Ya podemos hacer peticiones desde un navegador directamente a nuestros métodos Python. Pero ¿cómo se hace una aplicación con CherryPy?

Lo primero es descargar e instalar CherryPy, para lo cual hay información detallada y fácil de llevar a cabo en la documentación oficial.
Lo segundo es crear nuestro fichero principal de Python que hará las veces de servidor. En mi caso "server.py".
A éste fichero le añadimos el import de cherrypy:

    import cherrypy

Luego hacemos el diseño de las clases y métodos que serán dianas de las distintas URLs. Por ejemplo, para:

    http://mysite/myapp/exec

con un parámetro "myparam" que viene por ejemplo de un submit en HTML o en un GET:


class NameOfMyClass(Base):
    @cherrypy.expose
    def exec(self, myparam):
        # tratar los datos de la request, por ejemplo añadirlos a sesión
        cherrypy.session['param_01'] = myparam
        # realizar operaciones
        result = mymodel.mymethod(myparam)
        # devolver una respuesta, por ejemplo en formato HTML
        return html_generator(result)

Especialmente fijarse en que el decorator (@cherrypy.expose) hace al método exec visible para el mapeador de URLs, es decir, para el cliente HTTP en definitiva.
El nombre de la clase no tiene que coincidir con el de la aplicación en éste caso, precisamente por ser la clase de la aplicación, cuyo mapeo desde URL definiremos a continuación:

    cherrypy.config.update(conf)

    cherrypy.tree.mount(NameOfMyClass(), script_name='/myapp', app_conf)

    cherrypy.engine.start()

    cherrypy.engine.block()

Faltaría añadir aquí el contenido de las variables "conf" y "app_conf". En ambos casos pueden ser rutas a ficheros de configuración o diccionarios de Python. Sobre la configuración mejor consultar aquí o el libro que ya hemos comentado. Pero un ejemplo con diccionarios podría ser:


conf = {
    'global': {
        'server.socket_host': '127.0.0.1',
        'server.socket_port': 9091,
    },
}

app_conf = {
    '/style.css': {
        'tools.staticfile.on': True,
        'tools.staticfile.filename': os.path.join(_curdir,
        '/mydeploydir/css/style.css'),
    }
}

Y con esto sería suficiente. Ejecutamos:

    python server.py

Y ya deberíamos poder hacer peticiones desde un navegador a:

    http://127.0.0.1:9091/myapp/exec?myparam=eead

Por último, tocaba desplegar la aplicación tras el servidor Apache, de forma que este siguiera sirviendo otras aplicaciones (escritas en PHP, Perl, etc.) y CherryPy sirviera mi aplicación Python. Hay muchos ejemplos sobre como desplegar CherryPy tras Apache en el libro, y también en algunas páginas en caché de Google ;) pero para un profano de Apache no quedaba claro cuál de ellas era mejor para que ambos sirvieran páginas, ni cómo modificar los ejemplos para ello. Finalmente, para el esquema que muestro a continuación, añadiendo 2 líneas a la configuración de Apache fue suficiente.

    ProxyPass /myapp http://127.0.0.1:9091/myapp/
    ProxyPassReverse /myapp http://127.0.0.1:9091/myapp/

Esto es, un proxy reverso para la aplicación creada en CherryPy con script_name en el método cherrypy.tree.mount a "/myapp".

Y esto es todo lo que quería contar sobre CherryPy y servir aplicaciones web escritas en Python. Pythonic uh? :)