30 de diciembre de 2013

Curso NGS en Zaragoza (30-31 Enero 2014)

Hola,
nuestro colega Miguel Pérez-Enciso impartirá los días 30-31 de Enero (30, mañana y tarde; 31, mañana) un curso básico de análisis de datos NGS en el Salón de actos de la facultad de veterinaria de la Universidad de Zaragoza, invitado y organizado por Luis Varona. El curso es gratuito y no requiere inscripción.

El horario será:
jueves   30: 11-14H y 16-18H
viernes  31: 10-13H


Si consigo un enlace con el material del curso lo pondré aquí, de momento no sé mucho más que lo que aparece en el propio material que Miguel usa (http://www.icrea.cat/Web/OtherSectionViewer.aspx?key=1403&titol=2012):

"This is a short introductory course for NGS analyses.
The course assumes no knowledge on NGS data or data analyses, rather to the contrary, if you have experience , you will be disappointed. It is not yet a hands-on course, I will rather present the material theoretically but I am open with help in the analyses. I assume you are a computer illiterate biologist but that you have realized that you need computers in daily life. You are willing to:
Analyze by yourself the data
Willing to learn linux and perl or similar (python ...)

The main topics covered are:
FASTQ format and sequence quality (Phred score and fastqc software)
A simple pipeline: BWA + Samtools
Visualization: IGV
Identifying SNPs: Samtools
Annotating SNPs: ensembl vep.pl
Identifying structural variants
RNAseq: tophat + cufflinks
Population genetics: pools and individuals
Association studies with sequence: IBSeq

All topics are treated at very beginners level and superficially (check accompanying slides)"
Feliz año a todos,
Bruno

3 de diciembre de 2013

ofertas de postdoc en Bioninformática (Dic2013)

Hola,
comparto información de ofertas de trabajo en bioinformática que nos han llegado, hasta luego,
Bruno

1.
"We are looking to recruit a young, post-doctoral researcher for work on: Development and application of bioinformatic analysis methods for identification of biomarkers to assess the efficacy of TGF-β RI inhibitors in vivo"
http://www.genxpro.info/contact_us/Open_Positions/


2.

"[..]The Computational Biology Group of the Luxembourg Centre for Systems Biomedicine seeks a highly-skilled Postdotoral Fellow to work on an exciting project on reconstruction and analysis of an integrated gene regulatory network model to elucidate key mechanisms of cellular reprogramming. The model will rely on the integration and mining of diverse transcriptomics and epigenomics data of different cell types from the Central Nervous System. The
Postdoctoral Fellow is expected to supervise a Ph.D. student working on this project. Both of them are expected to collaborate with other members of the CBG to develop a computational methodology aiming at designing, in-silico, cellular reprogramming events, with a focus on the nervous system. This project will be carried out in collaboration with Prof. Noel Buckley’s lab at Kings College London.
Requirements of the ideal candidate:
• Ph.D. degree in Bioinformatics, Computer Science, Biology or a related discipline
• Strong computational skills
• Prior experience in mathematical modelling of biological networks, especially in network inference and analysis
• Excellent working knowledge in English.

We offer:
• Full contract for three years with possibility of renewal
• Opportunity to supervise a Ph.D. student working on the project
• Opportunity to do applied research to medical problems within a highly dynamic research institution (LCSB) and in collaboration with internationally recognized partners
• An exciting international environment
• A very competitive salary
For further information, please contact:
Prof. Dr. Antonio del Sol
E-mail: antonio.delsol@uni.lu
Applications should contain the following documents:
• A detailed curriculum vitae
• cover letter mentioning the reference number
• description of past research experience and future interests
• name and addresses of three referees
All applications should be sent preferably in electronic version until December 31 st, 2013 to the following address:
Luxembourg Centre for Systems Biomedicine (LCSB)
University of Luxembourg
7, avenue des Hauts-Fourneaux
L-4362 Esch-sur-Alzette
Tel: +352-466644-6982 (Office)
"

http://wwwen.uni.lu/lcsb/people/antonio_del_sol_mesa

21 de noviembre de 2013

GET_HOMOLOGUES for pan-genome analysis

Hola,
en el último número de Applied and Environmental Microbiology mi colega Pablo Vinuesa y yo publicamos un artículo describiendo el software GET_HOMOLOGUES, que tiene como abstract:
GET_HOMOLOGUES is an open source software package that builds upon popular orthology-calling approaches making highly customizable and detailed pan-genome analyses of microorganisms accessible to non-bioinformaticians. It can cluster homologous gene families using the bidirectional best-hit, COGtriangles or OrthoMCL clustering algorithms. Clustering stringency can be adjusted by scanning the domain-composition of proteins using the HMMER3 package, by imposing desired pair-wise alignment coverage cut-offs or by selecting only syntenic genes. Resulting homologous gene families can be made even more robust by computing consensus clusters from those generated by any combination of the clustering algorithms and filtering criteria. Auxiliary scripts make the construction, interrogation and graphical display of core and pan-genome sets easy to perform. Exponential and binomial mixture models can be fitted to the data to estimate theoretical core and pan-genome sizes, and high quality graphics generated. Furthermore, pan-genome trees can be easily computed and basic comparative genomics performed to identify lineage-specific genes or gene family expansions. The software is designed to take advantage of modern multiprocessor personal computers as well as computer clusters to parallelize time-consuming tasks. To demonstrate some of these capabilities, we survey a set of 50 Streptococcus genomes annotated in the Orthologous Matrix Browser as a benchmark case.
El  software  se puede descargar de http://www.eead.csic.es/compbio/soft/gethoms.php y también de http://maya.ccg.unam.mx/soft/gethoms.php y está escrito mayoritariamente en Perl, aunque contiene también trozos en R.
El manual del programa describe en detalle ejemplos de uso y está disponible en http://www.eead.csic.es/compbio/soft/manual.pdf .

Este paquete de programas se diseñó para el estudio de los pan y core-genomas de grupos de microorganismos, que es con lo que trabaja el grupo de Pablo fundamentalmente, y permite generar figuras como éstas:


Un saludo,
Bruno

9 de noviembre de 2013

Trucos para la biología computacional

Buenas,
hoy quiero invitaros a leer lo que nos cuentan dos bioinformáticos (Mick Watson y Nick Loman) sobre el trabajo y el aprendizaje de este oficio, publicado en Nature Biotechnology. Además de tocar temas más relacionados con el desarrollo de software (como el control de versiones) y la construcción de tuberías de análisis, el artículo repasa obviedades que no obstante conviene no olvidar, como
"knowledge of biology is vital in the interpretation of computational results"
 u otra más concreta:
 "Laboratory scientists wouldn’t dream of running experiments without the necessary positive and negative controls... tests are the computational biology equivalent".
El texto, breve, toca temas importantes como la elección apropiada de métodos en bioinformática, la validación de tu propio código y el que te descargaste de otros autores, y la búsqueda de opiniones expertas en foros como SEQanswers. Si te interesa, puedes seguir leyendo en http://www.nature.com/nbt/journal/v31/n11/full/nbt.2740.html?WT.ec_id=NBT-201311,
un saludo,
Bruno

14 de octubre de 2013

Premio Nobel en biología computacional

Buenos días,
han pasado ya unos días desde el anuncio, pero hasta hoy no he caído en que sería buena idea hablar en este blog sobre el reciente premio Nobel de Química, que ha caído sobre los hombros de Michael Levitt, Martin Karplus y Arieh Warshel. He de reconocer que al último no le conocía, pero a los dos primeros les he leído, visto en congresos y escuchado mencionar muchas veces a lo largo de mi aprendizaje, sobre todo en la vertiente estructural de la bioinformática.


Martin Karplus es famoso por sus trabajos sobre dinámica de macromoléculas, a menudo apoyados por datos cristalográficos,  y por el desarrollo de CHARMM. Por tanto, creo que tiene un perfil muy fuerte en química.

Sin embargo, Levitt ha tenido un papel muy relevante en la bioinformática estructural y en el desarrollo de algoritmos, y por tanto creo que es justo decir que ha sido muy influyente en el desarrollo de nuestra disciplina. Prueba de ello es que participa en los consejos editoriales de revistas que han sido clave en nuestra área como Journal of Molecular Biology, PNAS y PLoS Computational Biology, y que ha participado como evaluador en CASP. Además, como recoge una reciente nota de la ISCB, sus propias palabras traen un poco de este premio a la biología computacional:

“It’s sort of nice in more general terms to see that computational science, computational biology is being recognized. [...] It’s become a very large field and it’s always in some ways been the poor sister, or the ugly sister, to experimental biology.”

Así que estamos de enhorabuena, verdad?
Un saludo, Bruno

Añadido el 15/11/2013:
http://www.sciencedirect.com/science/article/pii/S0022283613006943

Añadido el 04/12/2013:
http://www.pnas.org/content/110/49/19656.extract.html



8 de octubre de 2013

What does an SNP look like?

So about to use samtools pileup or mpileup? Would you like to take a look at the file, instead of going straightforward for the BCF and bcftools automation?

Ok! this is what an SNP looks like!

contig_100029   698     C       39      .$.....,.,,,,.,,,.,.,..,......,,,,,,.,,^S. 
contig_100029   699     A       38      .....,.,,,,.,,,.,.,..,......,,,,,,.,,. 
contig_100029   700     A       38      .....,.,,,,.,,,.,.,..,......,,,,,,.,,. 
contig_100029   701     C       39      TTTTTtTttttTtttTtTtTTtTTTTTTttttttTttT^ST       
contig_100029   702     A       40      .....,.,,,,.,,,.,.,..,......,,,,,,.,,..^S,     
contig_100029   703     G       41      .....,.,,,,.,,,.,.,..,......,,,,,,.,,.,.,       
contig_100029   704     G       42      .....,.,,,,.,,,.,.,..,......,,,,,,.,,.,.,^S.   

I have ommited the last column (read base qualities).

And what about a reference skip? Maybe want to see an spliced alignment? Here you are!

contig_100029   516     T       43      ,,,.,,...,.,,.................,.,,,,.,,,.,^S.   
contig_100029   517     A       43      ,,,.,,...,.,,.................,.,,,,.,,,.,.     
contig_100029   518     G       43      ,,,.,,...,.,,.................,.,,,,.,,,.,.  
contig_100029   519     A       43      ,,,.,,...,.,,.................,.,,,,.,,,.,.     
contig_100029   520     G       43      ,,,.,,...,.,,.................,.,,,,.,,,.,.    
contig_100029   521     G       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>     
contig_100029   522     T       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>     
contig_100029   523     G       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>     
contig_100029   524     A       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>    
.    .    .    .    .    .    .    .    .    .    .    .    .    .    .
.    .    .    .    .    .    .    .    .    .    .    .    .    .    .
.    .    .    .    .    .    .    .    .    .    .    .    .    .    .
contig_100029   651     C       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>    
contig_100029   652     A       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>     
contig_100029   653     G       43      <<<><<>>><><<>>>>>>>>>>>>>>>>><><<<<><<<><>     
contig_100029   654     G       44      ,,,.,,...,.,,.................,.,,,,.,,,.,.^S, 
contig_100029   655     G       44      ,,,.,,..$.,.,,............$.....,.,,,,.,,,.,., 
contig_100029   656     C       42      ,,,.,,..,.,,................,.,,,,.,,,.,.,     
contig_100029   657     G       42      ,,,.,,..,.,,................,.,,,,.,,,.,.,     

Who wants an IGV when you can scroll up and down seeing how things pile down and up?
An it is memory efficient, right?
You just have to grep or awk filter your region of interest!

Do you want me to explain about columns?

1st: reference_name
2nd: position in reference
3rd: base on reference
4th: depth, or number of bases pilling-up over this reference position
5th: each symbol comes from a read, with some exceptions. Lets look at it:

"." and ",": matches! one in fwd strand, "," in reverse.
"$": end of read. Followed by mapping symbol of that read base ("." or "," for example). So, 2 symbols.
"^": start of read. Followed by mapping quality of the read (the MAPQ field in SAM format) and the mapping symbol of that base ("." or "," for example).
">" "<": reference skip. That is, the read still maps, but no in these reference bases. So probably an spliced alignment with part of the read aligning before the ">" and part after that. ">" fwd strand, "<" reverse.

6th: (not shown, just working with HQ bases ;)

Looking for more information? Ask please!

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.

20 de agosto de 2013

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.

30 de julio de 2013

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

17 de julio de 2013

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

29 de junio de 2013

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? :)


13 de junio de 2013

oferta de empleo en Bioinformática

Puesto de Postdoc en laboratorio de Genética de Plantas del CRAG. Se trata del grupo de genómica de melón del Dr. Jordi García-Mas:

Funciones: análisis bioinformático de datos de secuenciación de RNA de melón, en el marco de un proyecto PLANT-KBBE (SAFQIM: Sugar and Fruit Quality in Melon). Además, participará en la mejora de una base de datos de genómica de melón.

·         Requisitos: PhD en Biología o similar. Gran experiencia en bioinformática y Linux.

·         Se ofrece: contrato de 6 meses, renovable por 6 meses más. 37’5h semanales.

·         Inscripción: CRAG, Ref 18/2013 (antes del 30 de junio de 2013)



5 de junio de 2013

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



15 de mayo de 2013

Multicore con threads, threads::shared y Thread::Queue para controlar la RAM utilizada

Multicore con threads, threads::shared y Thread::Queue para controlar la RAM utilizada


Hace unos dias tanto Alvaro como Bruno publicaron unos blogs incluyendo scripts, tanto en Perl como en Python, en los cuales discutían sobre el uso y los problemas operativos al utilizar los modulos threads y threads::shared. Basicamente el problema principal es que, a veces es imposible utilizarlos ya que los problemas de duplicación de memoria hace que el programador no se pueda beneficiar de utilizar todos los procesadores de un ordenador, ya que se ve limitado por la cantidad de RAM que tiene disponible. Esa es una de las premisas que comentaba Bruno en su post. Yo me he encontrado multiples veces con este problema, por ejemplo procesando ficheros de genomas de Uniprot, que pueden tener ~ 30GB en sus casos mas extremos. En estos casos, al utilizar los modulos threads y threads::shared, normalmente lo que ocurre es que la memoria se va llenando mientras avanza el trabajo de las threads, y como resultado final el sistema operativo termina matando al proceso porque ocupa toda la RAM. He encontrado, mezclando código de varias personas, una forma de resolver este asunto utilizando estos modulos, y sin depender de la RAM disponible, combinando el uso de threads y threads::shared con Thread::Queue. De hecho, lo que se hace es generar una cola (en este caso una FIFO) en la cual un número predeterminado de 'worker threads' vayan sacando trabajo mientras vayan terminando. La idea es simultanear los procesos de llenado de la cola y trabajo, lo cual tambien permite controlar el tamaño de la cola, y por lo tanto la memoria utilizada. La idea sería, que esta aproximación podría combinarse con cualquier tarea que se necesite hacer, de hecho el código que compartiré a continuación, se podría adaptar fácilmente al compartido por Bruno para multi Blast.
Básicamente el programa a continuación hace lo que comentaba arriba, leer un fichero en formato Swissprot, y escanear las secuencias buscando un bias de glutamina y asparagina. El tamaño de la cola y por tanto la memoria utilizada se fija con la constante (Q_limit), en este caso definida para guardar 10000 lineas, aunque esto se puede modificar a conveniencia. Luego se generan 8 'worker threads', definidos por la constante (CORES), que van sacando trabajo mientras se llena la FIFO. Luego las 'worker threads' devuelven predicciones de secuencias con este tipo de bias y el programa principal las imprime en un fichero texto. Basado en este ejemplo, es muy facil adaptar esto a cualquier tarea que implique escanear secuencias en la búsqueda de cualquier propiedad, en ficheros enormes. Por ejemplo, un programa similar a este procesa en aproximadamente 2 horas un fichero de 35 GB con 64 CPUs y una utilización máxima de RAM de 5 GB.
Bueno aqui va el código, y espero sea útil la idea:
#!/usr/bin/perl -w

#Loading all the libraries to run in a multithreading environment
use warnings;
use threads;
use threads::shared;
use Thread::Queue;

my %aa_scores = (
   'A' => -0.149,
   'C' => -1.000,
   'D' => -0.395,
   'E' => -0.726,
   'F' => -0.125,
   'G' => 0.010,
   'H' => -0.034,
   'I' => -0.397,
   'K' => -0.496,
   'L' => -0.408,
   'M' => 0.044,
   'N' => 0.659,
   'P' => 0.059,
   'Q' => 0.536,
   'R' => -0.314,
   'S' => 0.192,
   'T' => -0.070,
   'V' => -0.450,
   'W' => -0.908,
   'Y' => 0.206,
   'O' => 0,
   'U' => 0,
   'B' => 0,
   'Z' => 0,
   'X' => 0,
);

use constant WINDOW => 60;
use constant CUTOFF => 12;
use constant CORES => 8; #Set the number of cores to use
use constant Q_limit => 10000; #Set the number of elements in the Queue

my $Q = new Thread::Queue (); #Initializing the Queue obect
my @threads;

#Loading the Swissprot file from command line argument
my $dat_file = shift;
$/ = "//\n"; #Reading the multiple lines of a Swissprot entry in a single line. The (//) line corresponds to the entry end makr 
#As this huge files are ussually distributed compressed here I pipe with zcat to avoid decompresing prior to processing
open (DAT_FILE, "zcat $dat_file |") or die "Sequence file couldn't be opened for reading\n";

#First generating the Queue
my $builder = threads->create (\&Queue_generator);

#And then setting up the worker threads
push @threads, threads->create (\&parse_swissprot) for (1 .. CORES);

#Wait until all the threads finish their work
$builder->join;

#Now that the work is done, collect the returned info extracted by the worker threads
#It is neccesary to fill a hash with all the partial results to be able to print the
#results with the scanning of the complete file in an organized fashion
my %domains;
foreach my $thread (@threads)
{
   my %partial_results = %{$thread->join};
 
   while (my ($organism, $predictions_ref) = each (%partial_results))
   {
      my @predictions;
      while (my ($seq_id, $info_ref) = each (%{$predictions_ref}))
      {
         my $score = $info_ref->{'Score'};
         my $domain = $info_ref->{'Seq'};
         my $window_pos = $info_ref->{'Window'};
   
         my $protein_id = $1 if ($seq_id =~ m/^>(\w+);/);
         my $protein_ac = $1 if ($seq_id =~ m/^>$protein_id;\s(\w+);\s/);
   
         push @predictions, "\t$protein_id|$protein_ac\tWindow Position=$window_pos; Score=$score | Predicted Domain: $domain\n";
      }
      push @{$domains{$organism}}, @predictions;#Loading all the predictions returned by each thread for a given organism
   }
}

  
#And finally printing the predictions of domains with a specific amino acid bias in the complete file
#organized by organisms
open (OUTFILE, ">predictions.txt");
while (my ($organism, $predictions_ref) = each (%domains))
{
   my @predictions = @{$predictions_ref};
 
   my $total_predictions = scalar (@predictions);
   print OUTFILE ">$organism: Total=$total_predictions\n@predictions\n";
}
close (OUTFILE);
close (DAT_FILE);

#Here I construct the Queue. The idea is that while the generator is feeding the queue other working threads
#are retrieving (emptying) from the queue and the memory used can be controlled during execution
sub Queue_generator
{
   while (my $line = <DAT_FILE>)
   {
      chomp ($line);
      while (1)
      {
         #The line is only send to processing if the number of items in the queue is less than a fixed size defined by Q_limit
         #thus a line is trapped in this infinite loop until the queue has space to receive it
         if ($Q->pending () < Q_limit) 
         {
            $Q->enqueue ($line);
            last; #Once the line has been loaded in the queue we need to break out of the infinite loop to read the next line from the file
         }
      }
   }
 
   #Once the file is exhausted it is needed to alert the worker threads that the work is over
   #When a thread catch and UNDEF variable it stops working
   $Q->enqueue (undef) for (1 .. CORES);
}

#This is the code of the worker threads
sub parse_swissprot
{
   my %predictions; #The variable to return to the main program
 
   #Here the thread DEQUEUE one line from the queue stack
   while (my $line = $Q->dequeue)
   {
      my @patterns = ();
      #Getting the Entry Name
      my $seq_id = $1 if ($line =~ m/^ID\s+(\w+)\s+/);
      
      # And now getting the Accession Number lines
      @patterns = $line =~ m/^AC\s+(.*)/mg;
      # and finally just getting the primary AC -i.e. the first of all the possible multiple ACs-
      my $seq_ac = $1 if ($patterns[0] =~ m/^(\w+);/);
      
      #And now getting the rest of the Entry info
      @patterns = $line =~ m/^DE\s+(.*)/mg;
      my $description = join ' ', @patterns;
      @patterns = $line =~ m/^SQ\s+(.*)/mg;
      my $seq_info = join ' ', @patterns;
      @patterns = $line =~ m/^OS\s+(.*)/mg;
      my $organism = join ' ', @patterns;
      @patterns = $line =~ m/^\s+(.*)/mg;
      my $seq = join '', @patterns;
      $seq =~ s/\s+//g;
      
      my $id = '>' . $seq_id . '; ' . $seq_ac . '; ' . $description . ' ' . $seq_info . "  [$organism]";
      
      next unless (length ($seq) >= WINDOW);
      my %results;
      
      #Loop to scan the complete sequence of the corresponding entry
      #After scanning the sequence only the highest scoring stretch is returned
      for my $i (0 .. (length ($seq) - WINDOW))
      {
         my $domain = substr ($seq, $i, WINDOW);
         my @domain = split (//, $domain);
   
         my $domain_score = 0;
         foreach my $aa (@domain)
         {
            $domain_score += $aa_scores{$aa};
         }
         $domain_score = sprintf ("%.3f", $domain_score);
   
         #Updating the $results variable to store only the highest scoring stretch
         if (exists ($results{$id}))
         {
            $results{$id} = {
               'Score' => $domain_score,
               'Seq' => $domain,
               'Window' => $i,
               'CompleteSeq' => $seq,
            } if ($domain_score > $results{$id}{'Score'});
         }
         else
         {
            $results{$id} = {
               'Score' => $domain_score,
               'Seq' => $domain,
               'Window' => $i,
               'CompleteSeq' => $seq,
            };
         }
      }
  
      #Loading the entry in the variable returned to the main program only if the highest scoring stretch passes the cutoff
      if ($results{$id}{'Score'} >= CUTOFF)
      {
         $predictions{$organism}{$id} = $results{$id};
      }
      $seq_id = $seq_ac = $description = $seq_info = $organism = $seq = undef;
   }
   return (\%predictions);
}



23 de abril de 2013

Multicore Blast con Python

El otro día Bruno publicó en una entrada del blog un script en Perl para optimizar el tiempo de cálculo de Blast en ordenadores con procesadores multicore. Un gran script, pero creo que se puede hacer más sencillo usando el módulo subprocess de Python.

Bruno usó el módulo Parallel::ForkManager de Perl. en el pasado ya habíamos mostrado en otra entrada las ventajas y peligros de paralelizar procesos con Perl con el módulo threads.

La solución con Python creo que es la más sencilla, consiste en usar el módulo subprocess. Otro módulo de Python que podríamos usar sería multiprocessing, pero no es muy sencillo retornar el standard output de programas externos como Blast para un procesamiento posterior de los resultados.

En el siguiente ejemplo se paraleliza la ejecución del comando 'blastp -version' que retorna un texto con la versión instalada de blastp. Dicho comando se puede cambiar por cualquier otro, así como añadir código al script para procesar los resultados. En el ejemplo se ejecuta el comando 10 veces con 4 procesos simultáneos, chequeando si han terminado cada 2 segundos.

 #!/usr/bin/python  
 # -*- coding: utf-8 -*-   
 #  
   
 # Importar módulos  
 import time  
 import os, sys  
 import subprocess  
 from pprint import pprint  
   
 maximum_process_number = 4 # Número máximo de procesos simultáneos  
 checking_processes_interval = 2 # (segundos) Tiempo de espera para checkear si han terminado los procesos  
 data_to_process = range(10) # Datos a procesar, como ejemplo una lista de números del 0 al 9  
   
 # El proceso a ejecutar (como ejemplo se muestran los datos de la version de Blastp)  
 # Se puede cambiar por cualquier otro programa y comando que no sea Blast  
 def run_process(processes) :  
      print "Running process"  
      processes.append(subprocess.Popen("blastp -version", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE))  
   
 # Programa principal  
 if __name__ == '__main__':  

      # Definir variables  
      count_processes = 0  
      finished_processes = []  
      processes = []  
      results = {}

      # Procesar múltiples datos  
      for i in data_to_process:

           # Chequear si los procesos han terminado (cuando hay tantos procesos como el máximo permitido)  
           while count_processes == maximum_process_number:  
                for p in processes:  
                     p.wait()  
                     if p.returncode == 0:  
                          results[p.pid] = "\n".join(p.stdout)  
                          processes.remove(p)  
                          count_processes -= 1  
                          print "Process finished"  
                     else:   
                          results[p.pid] = "\n".join(p.stdout)  
                          processes.remove(p)  
                          count_processes -= 1  
                          print "Process failed"  
                # Si no han terminado, esperar el tiempo establecido  
                if count_processes == maximum_process_number:  
                     print "Waiting for finished processes"  
                     time.sleep(checking_processes_interval)  

           # Añadir procesos a la cola de ejecución  
           run_process(processes)  
           count_processes += 1  
           print "Process started", i, count_processes  

      # Chequear si los procesos han terminado los últimos procesos  
      while count_processes != 0:  
           for p in processes:  
                p.wait()  
                if p.returncode == 0:  
                     results[p.pid] = "\n".join(p.stdout)  
                     processes.remove(p)  
                     count_processes -= 1  
                     print "Last processes finishing"  
                else:   
                     results[p.pid] = "\n".join(p.stdout)  
                     processes.remove(p)  
                     count_processes -= 1  
                     print "Process failed"  
           if count_processes == maximum_process_number:  
                print count_processes,maximum_process_number  
                print "Waiting for last finished processes"  
                time.sleep(checking_processes_interval) 
 
      # Imprimir resultados  
      print "Output from the processes:"  
      pprint(results)