Jornadas Bioinformáticas JBI 2012 (XI Edición)

Aprovecho esta última entrada del año para dar difusión a las Jornadas de Bioinformática JBI 2012. Como muchos sabréis, este congreso es el principal punto de encuentro anual de nuestra comunidad en la península ibérica, así que nuestro laboratorio también estará en Barcelona del 23 al 25 de Enero. El programa completo de las jornadas se puede descargar en este enlace.

http://sgu.bioinfo.cipf.es/jbi2012

Este año presentaremos parte de nuestro trabajo reciente:

"Genome-wide clustering of transcription factors by comparison of predicted protein-DNA interfaces"

donde explicamos y evaluamos la anotación de interfaces de reconocimiento de DNA en secuencias de proteínas por medio de diferentes aproximaciones como BLAST, TFmodeller, DP-Bind y DISIS.






El tema principal de las jornadas será "Arquitectura genómica, anotación y diseño", sobre el cual se discutirán los diferentes avances en la integración de los campos de la Biología, Medicina e Informática en el campo de la Genómica. Además se tratarán los siguientes temas:
- Análisis de datos de secuenciación de alto rendimiento (NGS)
- Bioinformática estructural
- Algoritmos de biología computacional y computación de alto rendimiento
- Análisis de sequencias, filogenética y evolución
- Bases de datos, herramientas y tecnologías en biología computacional
- Bioinformática en transcriptómica y proteómica
- Biología de sistemas

ENGLISH:

The XIth Spanish Symposium on Bioinformatics (JBI2012) will take place in January 23-25, 2012 in Barcelona, Spain. Co-organised by the Spanish Institut of Bioinformatics and the Portuguese Bioinformatics Network and hosted by the Barcelona Biomedical Research Park (PRBB). The full program can be downloaded from this link.

This year, the reference topic is “Genome Architecture, Annotation and Design” for which the conference will provide the opportunity to discuss the state of the art for the integration of the fields of biology, medicine and informatics. We invite you to submit your work and share your experiences in the following topics of interest including, but not limited to:

- Analysis of high throughput data  (NGS)
- Structural Bioinformatics
- Algorithms for computational biology and HPC
- Sequence analysis, phylogenetics and evolution
- Databases, Tools and technologies for computational biology
-  Bioinformatics in Transcriptomics and Proteomics
- System and Synthetic Biology

Our contribution to the congress:

Genome-wide clustering of transcription factors by comparison of predicted protein-DNA interfaces

Transcription Factors (TFs) play a central role in gene regulation by binding to DNA target sequences, mostly in promoter regions. However, even for the best annotated genomes, only a fraction of these critical proteins have been experimentally characterized and linked to some of their target sites. The dimension of this problem increases in multicellular organisms, which tend to have large collections of TFs, sometimes with redundant roles, that result of whole-genome duplication events and lineage-specific expansions. In this work we set to study the repertoire of Arabidopsis thaliana TFs from the perspective of their predicted interfaces, to evaluate the degree of DNA-binding redundancy at a genome scale. First, we critically compare the performance of a variety of methods that predict the interface residues of DNA-binding proteins, those responsible for specific recognition, and measure their sensitivity and specificity. Second, we apply the best predictors to the complete A.thaliana repertoire and build clusters of transcription factors with similar interfaces. Finally, we use our in-house footprintDB to benchmark to what extent TFs in the same cluster specifically bind to similar DNA sites. Our results indicate that there is substantial overlap of DNA binding specificities in most TF families. This observation supports the use of interface predictions to construct reduced representation of TF sets with common DNA binding preferences.

Felices Navidades bioinformáticas!!!

Ya que estamos en Navidad os dejo nuestra felicitación navideña especialmente pensada para los bioinformáticos... con algo de inspiración en estructura tridimensional de proteínas y dianas terapéuticas.

Adjunto la imagen a alta resolución para que podáis intentar adivinar con qué estructura del PDB se han generado las bolas (proteínas) y adornos (ligandos) que decoran el árbol, dejar en los comentarios vuestras respuestas...

Y por supuesto... FELIZ NAVIDAD Y PRÓSPERO AÑO 2012!!!

English version:

In this Christmas time we want to share our Christmas card especially designed for bioinformaticians ... with some inspiration from three-dimensional structure of proteins and therapeutic targets.

It is a high resolution image so you can try to guess from which PDB structure have been extracted the balls (proteins) and ornaments (ligands) that decorate the tree, give your answers in the comment section ...

And of course... MERRY CHRISTMAS AND HAPPY NEW YEAR 2012!!!

metaDBSite: un metaservidor para predecir la unión de proteínas a DNA

Los meta-servidores son valiosas herramientas en bioinformática para no perder el tiempo mandando trabajos a diferentes servidores y tener que juntar después manualmente los resultados. Un meta-servidor nos ofrece utilizar varios servicios web a la vez desde una interfaz única. Ejemplos de meta-servidores son: NCBI, UniProt, GeneSilico...

Recientemente he  tenido que usar diferentes programas de predicción de interfaces de unión proteína-DNA, y he encontrado la grata sorpresa de que existía un meta-servidor que permitía realizar la predicción de una sola vez utilizando la mayoría de técnicas conocidas, este valioso meta-servidor se llama metaDBSite. No sólo eso, si no que he podido comprobar que siguen dándole soporte.

Os animo a probarlo con una sencilla secuencia proteica:

>1a02_N  
 wplssqsgsyelrievqpkphhRahYetEgsRgavkaptgghpvvqlhgymenkplglqifigtaderilkphafyqvhritgktvtttsyekivgntkvleiplepknnmratidcagilklrnadielrkgetdigRkntrvrlvfrvhipessgrivslqtasnpiecsQRsahelpmverqdtdsclvyggqqmiltgqnftseskvvftekttdgqqiwemeatvdkdksqpnmlfveipeyrnkhirtpvkvnfyvingkrkrsqpqhftyhpv  

Si comprobáis los resultados con los datos de la estructura original (cadena N de la estructura 1a02) podréis ver cómo el mejor de los métodos (BindN-RF) 'adivina' 7 de los 7 residuos que contactan el DNA entre los 18 predichos, si tenemos en cuenta que la proteína tienen 280 aminoácidos, el resultado no está nada mal...

Sincronizando el Protein Data Bank

Hola,
para los que trabajamos de manera habitual con archivos del Protein Data Bank es muy conveniente tener acceso a la colección de estructuras en formato PDB. Ojo, necesitarás al menos unos 14Gb de disco, pero tendrás al  alcance de tu mano más de 77 mil estructuras de proteínas (y sus ligandos). La sección de descargas del PDB ofrece varias opciones, pero en mi opinión, si necesitas acceder desde tus programas a los archivos PDB sin supervisión, lo más eficiente es tener una copia local, lo que en la jerga se llama mirror o sitio espejo. Como el PDB se va actualizando semanalmente y no sólo se añaden entradas nuevas, si no que se van corrigiendo o completando estructuras antiguas, el proceso de sincronización no es tan sencillo como pudiera parecer. Es aquí donde nos salva la vida el software rsync, parte integrante de cualquier Linux moderno, que permite hacer estas actualizaciones de manera eficiente, y que fácilmente podemos añadir a un cron semanal:

 #!/usr/bin/perl -w   
 use strict;  
   
 my $PDBSERVER = 'rsync.wwpdb.org';  
 my $PDBPORT  = 33444;  
 my $LOCALPDBDIR = '/path/to/my/PDB/mirror/';  
 my $LOGFILE   = '/path/to/PDBmirror.log';  
   
 my $rsync_command = "/usr/bin/rsync -rlpt -v -z --delete --exclude=biounit ".  
      "--exclude=monomers --exclude=status --exclude=obsolete --exclude=models ".  
      "--exclude=mmCIF --exclude=nmr_restraints* --exclude=structure_factors " .  
   "--exclude=XML --exclude=XML-extatom --exclude=XML-noatom ".  
      "--port=$PDBPORT $PDBSERVER\::ftp_data/ $LOCALPDBDIR";  
   
 open(LOG,">$LOGFILE") || die "# $0 : cannot write to $LOGFILE\n";  
 open(RSYNC,"$rsync_command 2>&1 |") || die "# $0 : cannot execute $rsync_command\n";  
 while(<RSYNC>)  
 {  
   print LOG;   
   if(/rsync: failed to connect to/)  
   {  
       die "# $0 : error:\n$_\n";  
   }  
 }  
 close(RSYNC);  
 close(LOG);  

Un saludo, Bruno

bajar de una sola vez todos los proteomas (ensembl)

Me informan en helpdesk@ensembl.org que, para bajarse la última versión de todos los proteomas a la vez, pueden utilizarse las siguientes lineas (en bash):

---------------------

#!/usr/bin/sh

version='64'

prefix=rsync://ftp.ensembl.org/ensembl/pub/release-${version}

for f in $(rsync --verbose --recursive $prefix/fasta | grep '.*pep.*.fa.gz' | perl -n -e 'print [split(/\s+/)]->[4]."\n"');

do

rsync --partial --verbose $prefix/$f .

done

----------------------

Como debe haber muchas más formas de hacer esto, abro la lata y vais añadiendo.

Hope this helps,

JR

'Fold it' o como el juego agudiza el ingenio

Buenas,
desde hace al menos una década el grupo de David Baker (U.Washington) ha tenido un impacto tremendo en el estudio del plegamiento de las proteínas, yo diría que hasta el punto de ir cocinando poco a poco un premio Nobel, pero el tiempo dirá. Muy brevemente, su trabajo ha consistido en desarrollar algoritmos para la predicción de la estructura de proteínas (por comparación y ab initio) que han ido refinando en sucesivas rondas de CASP (donde siempre quedan entre los mejores) y poniendo a punto con experimentos bastante espectaculares donde comparan sus predicciones con estructuras obtenidas por cristalografía y NMR, obteniendo modelos de una calidad y complejidad sin precedentes, llegando por ejemplo a diseñar una proteína artificial y predecir su estructura terciaria con un error promedio de 1.2A. Una de las aportaciones más extravagantes de este laboratorio tan productivo fue la creación de un juego basado en el problema del plegamiento, que llamaron FoldIt.

El juego permite poner a jugadores de todo el mundo a resolver ejercicios de plegamiento de dificultad variable, tras superar un tutorial donde aprendes las reglas del juego y las herramientas disponibles. La verdad es que hasta hoy no he encontrado una manera más amena de entender qué es un puente de hidrógeno o un choque estérico. Enfin, a lo que voy, lo que parecía simplemente un intento por popularizar el 'problema del plegamiento', que a día de hoy todavía se considera sin resolver, se ha convertido también en una herramienta de innovación, como se publica en el último número de PNAS. Tras estudiar la evolución de las estrategias ganadoras encontradas por la comunidad jugadores, que se pueden guardar como recetas, se observa que se parecen mucho a los propios algoritmos que el grupo de Baker está desarrollando. O dicho de otro modo, parecería que poniendo en forma de juego un problema complejo como esto, con puntuaciones objetivas y reglas sencillas, el ingenio colectivo puede llegar a soluciones parecidas a las de los científicos especialistas, sorprendente?
Un saludo,
Bruno

grafos de De Bruijn

Hola,
un artículo educativo publicado recientemente en Nature Biotechnology me ha vuelto a recordar los grafos de De Bruijn, que ya había mencionado de pasada en una entrada anterior sobre vectores de sufijos. Hoy los veremos un poco más de cerca, por su importancia para la reconstrucción de secuencias genómicas a partir de fragmentos más pequeños que han sido previamente secuenciados.

Los modernos ensambladores de secuencias emplean grafos de De Bruijn para guardar en memoria los prefijos (nodos) y sufijos (aristas) de las lecturas o reads obtenidas en experimentos de secuenciación. Este tipo de representación permite reducir en parte la redundancia natural de este tipo de datos a la vez que facilita su ensamblaje posterior. La figura  compara los métodos de ensamblaje tradicionales (a) con los actuales (d). Los tradicionales se usan para pegar entre si secuencias obtenidas por el método de Sanger, que suelen ser largas, de buena calidad y en números pequeños. Los métodos que acompañan a los secuenciadores actuales se adaptan al nuevo escenario de secuencias cortas y en números varios órdenes de magnitud más elevados.

fuente: http://www.nature.com/nbt/journal/v29/n11/full/nbt.2023.html

Como se muestra en el panel d de la figura, un grafo de De Bruijn facilita la reconstrucción del genoma de partida, que se obtiene buscando un ciclo euleriano que recorra todas las aristas (sufijos) una sola vez. Esta estrategia es computacionalmente mucho más barata que la alternativa de buscar ciclos hamiltonianos, que en la práctica es imposible al ser un problema de tipo NP completo. Sin embargo, como pasa a menudo, la realidad se resiste a ser capturada por una sola teoría, e incluso los grafos de De Bruijn son todavía incapaces de resolver el problema de las secuencias repetidas, que aparecen en un genoma o región genómica no una sinó muchas veces.

En el Russell's blog encontraréis código fuente en Python para aprender cómo se construyen estos grafos. Espero haberos despertado la curiosidad, un saludo,
Bruno

Mapamundi NGS

¿Te has preguntado alguna vez dónde puedes secuenciar ácidos nucleicos?
¿Estás buscando un servicio que ofrezca 'runs' en el Polonator o en Illumina HiSeq?
¿Quieres saber cuántos secuenciadores de DNA hay en tu ciudad y qué modelos son?
¿Dónde está el secuenciador más septentrional del planeta? ¿Y el más meridional?
¿Hay secuenciadores en la Antártida o en Hawai?

Hola,
en esta nueva entrada recomendamos un pequeño pasatiempos, una aplicación web basada en Google Maps que nos permite visualizar la localización de los distintos centros de secuenciación y, más concretamente, de los aparatos de secuenciación (NGS, claro) a lo largo y ancho del planeta.

Next Generation Genomics: World Map of High-throughput Sequencers

 

Y no sólo permite visualizar, ya que pretende que sean los usuarios los que añadan los centros y/o secuenciadores, y corrijan las localizaciones. Así que ya sabes, si conoces un secuenciador que no aparece en el mapa añádelo. Además, cuenta con un sencillo filtro en la parte superior, con el cual visualizar solamente los modelos de interés o hacer una rápida comparativa; y con un desplegable para poder ir directamente a un país concreto.

La aplicación fue creada mediante Google Maps API, JQuery, MarkerClusterer, LabeledMarker, Django y SQLite, en la Universidad de Birmingham.

Un saludo, hasta una nueva entrada,
Carlos

Adiós a ClustalW

Hola,
en esta entrada quería enterrar al ya vetusto ClustalW, el programa de alineamiento múltiple más citado de la historia. En realidad lo acaban de enterrar sus propios autores en este artículo, donde presentan a su sucesor Clustal Omega. Como es habitual en Bioinformática, los autores ponen a prueba el nuevo programa de alineamiento comparándolo con las principales opciones de software disponibles, incluyendo al viejo ClustalW. Las comparaciones son bastante extensas, usando hasta 3 baterías de alineamientos, como son BALiBASE y HomFam (creadas por el entorno de Clustal) y Prefab, creada por el autor de MUSCLE, Robert Edgar. En la siguiente tabla, abreviada de la original, se muestran los resultados promediados sobre el conjunto BALiBASE de 218 familias de proteínas:

Aligner	Av score (218 families)							Tot time (s)
MSAprobs	0.607							12 382.00
Probalign	0.589							10 095.20
MAFFT (auto)	0.588							1475.40
Probcons	0.558							13 086.30
Clustal Ω	0.554							539.91
T-Coffee	0.551							81 041.50
Kalign	0.501							21.88
MUSCLE	0.475							789.57
MAFFT (default)	0.458							68.24
FSA	0.419							53 648.10
Dialign	0.415							3977.44
PRANK	0.376							128 355.00
ClustalW	0.374							766.47

En cuanto a poder decir si efectivamente Clustal Omega es el mejor alineador disponible en la actualidad habrá que esperar el veredicto de la comunidad y ver qué le parecen a Rober Edgar los parámetros de MUSCLE que se usaron en las pruebas, pero no parece superar en precisión a MAFFT en modo automático.
Sin embargo, las tablas del artículo (1 , 2 , 3) muestran que la inclusión de modelos ocultos de Markov (HMMs) y árboles guía mBed produce alineamientos de mucha mayor calidad (y en menor tiempo) que ClustalW, que es más lento para producir peores alineamientos. Clustal Omega puede alinear conjuntos de miles de secuencias, pero por el momento, sólo de aminoácidos, por lo que talvez sigamos usando ClustalW o MAFFT para alinear nucleótidos, no?

Por cierto, al compilar el código fuente de http://www.clustal.org/omega/clustal-omega-1.0.3.tar.gz en mi Ubuntu 10.04 tuve que instalar el paquete libargtable2-dev , un saludo,
Bruno

Perl Moderno?

Hoy me gustaría recomendar 'Modern Perl', que es un librito curioso porque:

1) contiene muchas recetas explicadas a base de ejemplos sobre los diferentes aspectos del lenguaje Perl, reflexionando sobre el famoso principio 'there is more than way to do it' y sobre las ventajas de unas maneras sobre otras de hacer la misma operación. Por ejemplo, se explica que el siguiente bucle:

 my @tripled;  
 my $count = @numbers;  
 for (my $i = 0; $i < $count; $i++)  
 {  
   $tripled[$i] = $numbers[$i] * 3;  
 }  

se puede expresar también como:

 my @tripled = map { $_ * 3 } @numbers;  

También por ejemplo muestra que igual que escribimos a un archivo podemos hacerlo a una variable escalar:

 my $captured_output;  
 open(OUTSTRING,'>',\$captured_output);  

2) es adecuado para principiantes del lenguaje y para gente más experta, ya que reflexiona sobre cuestiones de estilo y eficiencia, y además propone módulos clave de CPAN para facilitarnos la vida y hacer el código más robusto y portable

3) los autores del libro han dispuesto su libre descarga en diferentes formatos, por ejemplo en PDF

Que aproveche,
Bruno

Regulación por microRNAs exógenos de la dieta?

Hola,
uno de los temas principales que estudiamos en nuestro laboratorio es la regulación de la expresión génica, que como sabéis puede darse a varios niveles.
Si hace poco nuestro colega Lorenz Bülow nos contaba en un seminario de la EEAD  la integración de la regulación transcripcional y postranscripcional en la planta modelo Arabidopsis thaliana, organizada en la base de datos relacional AthaMap, hoy descubro un artículo reciente en Cell Research donde los autores publican evidencia de la presencia de cerca de 30 microRNAs de arroz en muestras de sangre de poblaciones humanas y de ganado en China. El artículo, que parte de muestreos masivos de secuencias y luego confirma los resultados por PCR, sugiere de manera convincente que algunos microRNAs expresados en el grano de arroz, dieta fundamental de las poblaciones estudiadas, pueden regular la expresión génica de sus comensales. Seguro que el artículo será puesto a prueba en posteriores análisis y estudios, para validarlo de manera inequívoca, porque estas observaciones desvelan dos hechos que sin duda tendrán mucho impacto:
1) puede haber transferencia de ácidos nucleicos de las plantas a los mamíferos que se las comen, a pesar de la digestión
2) puede haber fenómenos de regulación genética en la naturaleza a través de la dieta sin que medien hormonas, directamente por microRNAs, moléculas en torno a los 22 ribonucleótidos de tamaño, capaces de atravesar los epitelios del tracto digestivo

fuente: http://mcb.berkeley.edu/labs/he/Research.htm

Recomiendo la lectura de la fuente original y si tenéis algo que añadir por favor usad los comentarios, un saludo,
Bruno

Guía de campo de tecnologías de secuenciación

Hola,
ayer me encontré en la Red una revisión, publicada en Mayo de 2011 por TC Glenn, que contiene la siguiente tabla, muy útil  para comparar de un vistazo las plataformas de secuenciación de segunda generación disponibles actualmente:

Tabla original publicada en Molecular Ecology Resources

Esta tabla se complementa con otras disponibles en la 'NGS Field Guide', actualizadas regularmente, incluyendo por ejemplo los costes y los tipos de errores más frecuentes en cada una de ellas. De hecho habrá que esperar para tener datos empíricos de los errores típicos de la plataforma IonTorrent, que por ahora se basan en datos proporcionados por la compañía (previas a su publicación en Nature el pasado mes de Julio).
Hasta otra, Bruno

Leyendo archivos comprimidos .gz

Hola,
tras la entrada 'Compresión de secuencias de ADN', que ya utilizaba el módulo estándar Compress::Zlib, en este ejemplo se muestra como leer línea a línea un archivo comprimido con los algoritmos de la librería zlib. En particular este ejemplo lee un archivo FASTA de gran tamaño sin necesidad de descomprimirlo entero.  Lo he probado con éxito en Linux (con el intérprete Perl que viene instalado en Ubuntu 10.04) y también en Windows (con ActiveState Perl).

 use strict;  
 use Compress::Zlib;  
 # http://perldoc.perl.org/Compress/Zlib.html  
   
 my $filename = '/path/to/swissprot.gz';  
 # ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/  
   
 my $gz = gzopen($filename,'rb') ||   
    die "# Cannot open $filename: $gzerrno\n" ;  
 while($gz->gzreadline($_) > 0)   
 {  
   # haz algo con esta linea, guardada en $_, por ejemplo   
   # imprime los identificadores GenBank encontrados  
    if($_ =~ /^>/ && /(gi\|\d+)/){ print "$1\n" }  
 }  
   
 if($gzerrno != Z_STREAM_END)  
 {  
    die "# Error reading from $filename: $gzerrno\n"  
 }  
   
 $gz->gzclose();  

Hasta pronto, Bruno

Decodificando la Gene Ontology (C/C++)

Hola,
como continuación de la entrada anterior os pongo hoy código en C/C++ para hacer exactamente la misma tarea, aprovechando las estructuras de datos e iteradores de la Standard Template Library (STL). En mis pruebas, esta implementación es aproximandamente el doble de rápida que la de Perl. Por cierto, se compila con:  
$ gcc -o get_go_annotation get_go_annotation.cpp myGOparser.cpp -lstdc++

El código fuente incluye 3 ficheros:


1) get_go_annotation.cpp

 #include <stdio.h>  
 #include <stdlib.h>  
 #include <string>  
 #include <map>  
 #include "myGOparser.h"  
 using namespace std;  
   
 #define DEFFLATFILE "/path/to/gene_ontology.1_2.obo";  
   
 int main(int argc, char *argv[])  
 {   
    if(argc == 1)  
    {  
       printf("# usage: ./get_go_annotation <GO:0046983>\n\n");  
       exit(-1);  
    }  
      
    string input_term = argv[1];  
    string flatfile_name = DEFFLATFILE;   
    int n_of_records = 0;  
    map <string, GOnode> parsedGO;  
      
    printf("# parsing GO flat file ... ");  
    n_of_records = parse_GO_flat_file(flatfile_name, parsedGO);  
    printf("done (%d records)\n\n",n_of_records);  
      
    string annot = get_full_annotation_for_term(input_term,parsedGO);  
    printf("%s\n",annot.c_str());  
      
    exit(0);  
 }  

2) myGOparser.h

 /* Bruno Contreras-Moreira, 2011 EEAD/CSIC */   
 #include <stdio.h>  
 #include <stdlib.h>  
 #include <string>  
 #include <vector>  
 #include <map>  
   
 using namespace std;  
   
 #define LINE 400  
 #define FIELD 200  
   
 struct GOnode   
 {  
    string name;  
    vector <string> parents;     
 };  
   
 int parse_GO_flat_file(string &filename, map <string, GOnode> &GO);  
   
 string get_full_annotation_for_term(string &term, map <string, GOnode> &GO);  

3) myGOparser.cpp

 // Bruno Contreras-Moreira, 2011 EEAD/CSIC  
 #include <string.h>  
 #include <map>  
 #include "myGOparser.h"  
   
 using namespace std;  
   
 int parse_GO_flat_file(string &filename, map <string, GOnode> &GO)  
 {  
    // 1) parse flat GO file  
    FILE * gofile = fopen(filename.c_str(),"r");  
    if(gofile == NULL)  
    {  
       printf("# parse_GO_flat_file : cannot read %s, exit ...\n",  
          filename.c_str());  
       return 0;  
    }  
      
    int n_of_records = 0;  
    char line[LINE], field[FIELD];  
    string term,name,alt_id,parent;  
    map <string, string> synonym;  
    while( fgets(line, LINE, gofile) != NULL )  
    {  
      /*[Term]  
       id: GO:0006355  
       name: regulation of transcription, DNA-dependent  
       namespace: biological_process  
       alt_id: GO:0032583  
       is_a: GO:0010468 ! regulation of gene expression */  
         
       if(strncmp(line,"id:",3) == 0)  
       {   
          sscanf(line,"id: %s", (char *) &field);  
          term = field;   
          n_of_records++;  
       }  
       else if(strncmp(line,"name:",5) == 0)   
       {  
          strncpy(field,line+6,FIELD-1);  
          field[strcspn (field,"\n")] = '\0'; // chomp  
          name = field;   
          GO[term].name = name;  
       }  
       else if(strncmp(line,"alt_id:",7) == 0)  
       {  
          sscanf(line,"alt_id: %s",(char *) &field);  
          alt_id = field;  
          synonym[alt_id] = term;   
       }  
       else if(strncmp(line,"is_a:",4) == 0)  
       {  
          sscanf(line,"is_a: %s",(char *) &field);  
          parent = field;     
          GO[term].parents.push_back(parent);  
       }  
    }  
    fclose(gofile);  
      
    // 2) link synonims  
    map <string, string>::iterator syn;  
    for(syn = synonym.begin(); syn != synonym.end(); ++syn )   
       GO[syn->first] = GO[syn->second];  
      
    return n_of_records;  
 }  
   
 string get_full_annotation_for_term(string &term, map <string, GOnode> &GO)  
 {  
    string annot, pterm;  
    vector <string>::iterator pit;  
      
    if(GO.find(term) == GO.end())  
    {  
       annot = "[cannot find GO term ";  
       annot += term;  
       annot += "]";  
    }  
    else  
    {  
       if(GO[term].parents.size() == 0)  
       {  
          annot += GO[term].name;  
          annot += "(";  
          annot += term;  
          annot += ") | ";  
       }  
       else  
       {  
          for(pit=GO[term].parents.begin();pit != GO[term].parents.end();++pit)  
          {      
             annot += GO[term].name;  
             annot += "(";  
             annot += term;  
             annot += "), ";  
             annot += get_full_annotation_for_term(*pit,GO);  
            
          }  
       }     
    }  
      
    return annot;  
 }  

Un saludo, Bruno

Decodificando la Gene Ontology

Hola,
a estas alturas es difícil no haber oído hablar de la Gene Ontology (GO), el vocabulario controlado más utilizado en Biología para describir secuencias y sus funciones. Diría yo que la principal aplicación de la GO es facilitar la anotación, manipulación y comparación computacional de grandes volúmenes de secuencias. De hecho, los repositorios de referencia, como por ejemplo UniProt, llevan tiempo anotando sus secuencias usando la jerarquía de la GO. Por lo tanto, actualmente es fácil encontrarse con un montón de secuencias, por ejemplo procedentes de un experimento o generadas por un colaborador, que en vez de descripciones legibles tienen asociados códigos de la GO. Por ejemplo, puedes encontrarte con una proteína de A.thaliana,

> AT3G04220  GO:0005524  

ATGGATTCTT CTTTTTTACT CGAAACTGTT GCTGCTGCAA CAGGCTTCTT
CACACTTTTG GGTACAATAC TTTTTATGGT TTACAGAAAA TTCAAAGACC
ATCGAGAAAA CAAAGAAAAT GATTCTTCTT CTTCAACACA ...

y preguntarte, qué será? Cuál es su función? En el Taller de (bio)Perl ya apuntamos una manera de averiguar la genealogía del término GO:0005524, pero dadas las limitaciones del módulo GO::Parser aquí os muestro una manera más eficiente de extraer el significado y la genealogía de cualquier identificador de la GO, que puede pertenecer a una de las 3 ramas fundamentales de la jerarquía  (proceso biológico, función molecular y localización celular). Para ello es necesario descargar la versión más reciente de la GO en un archivo de texto en formato OBO, y guardar el siguiente módulo escrito en Perl:

 package myGOparser;  
 require Exporter;  
   
 @ISA = qw(Exporter);  
 @EXPORT = qw( $DEFGOFLATFILE parse_GO_flat_file get_full_annotation_for_id );   
 use strict;  
   
 our $DEFGOFLATFILE = "/path/to/GO/gene_ontology.1_2.obo";  
   
 sub parse_GO_flat_file  
 {  
    my $flatfile = $_[0] || $DEFGOFLATFILE;  
    my ($id,$alt_id,$name,$namespace,%synonim,%GO);  
      
    open(GOFILE,$flatfile) ||   
       die "# parse_GO_flat_file : cannot read $flatfile\n";  
    while(<GOFILE>)  
    {  
       #[Term]  
       #id: GO:0006355  
       #name: regulation of transcription, DNA-dependent  
       #namespace: biological_process  
       #alt_id: GO:0032583  
       #is_a: GO:0010468 ! regulation of gene expression  
       if(/^id: (GO:\d+)/){ $id = $1; }  
       elsif(/^name: (.+?)\n/){ $GO{$id}{'name'} = $1; }  
       elsif(/^alt_id: (GO:\d+)/){ push(@{$synonim{$id}},$1); }  
       elsif(/^is_a: (GO:\d+)/){ push(@{$GO{$id}{'parents'}},$1); }  
    }  
    close(GOFILE);  
      
    # link synonims  
    foreach $id (keys(%synonim))  
    {  
       foreach my $syn (@{$synonim{$id}}){ $GO{$syn} = $GO{$id} }  
    }  
      
    return \%GO;  
 }  
   
 sub get_full_annotation_for_id  
 {  
    my ($id,$ref_parsed_GO) = @_;  
    my $annot;  
      
    if(!$ref_parsed_GO->{$id})  
    {  
       return "[cannot find GO term $id]";  
    }  
    else  
    {  
       if(!$ref_parsed_GO->{$id}{'parents'})  
       {  
          $annot .= $ref_parsed_GO->{$id}{'name'}."($id) | ";  
       }  
       else  
       {  
          foreach my $pid (@{$ref_parsed_GO->{$id}{'parents'}})  
          {  
             $annot .= $ref_parsed_GO->{$id}{'name'}."($id)".', '.  
               get_full_annotation_for_id($pid,$ref_parsed_GO);  
          }  
       }     
    }  
      
    return $annot;  
 }  
 1;  

Ahora podemos crear un programa como éste, que invocado desde el terminal nos devuelve la descripción jerárquica completa de un termino GO:

 use strict;  
 use myGOparser;  
   
 my $input_id = $ARGV[0] || die "# usage: ./go_parser.pl \n\n";  
   
 print "# parsing GO flat file ... ";  
 my $ref_GO = parse_GO_flat_file();  
 print "done\n\n";  
   
 print get_full_annotation_for_id($input_id,$ref_GO)."\n";  

Para el identificador GO:0046983 obtendremos la siguiente descripción:

protein dimerization activity(GO:0046983), protein binding(GO:0005515), binding(GO:0005488), molecular_function(GO:0003674) |

Transformada de Burrows-Wheeler (para comprimir secuencias y buscar patrones)

Hola,
como posiblemente haya vacaciones en el blog estas próximas semanas, os dejo una entrada algo más compleja, donde se muestran los fundamentos de la transformación de Burrows-Wheeler aplicada a la búsqueda de patrones en secuencias biológicas, como continuación natural de los vectores de sufijos y los árboles binarios, que vimos en una entrada anterior.

(fuente: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.12.1158)

Vaya por delante que las implementaciones reales de este algoritmo (como Bowtie o BWA) son bastante más complejas (porque permiten búsquedas inexactas) y eficientes (porque comprimen el corpus para que se pueda cargar en la memoria RAM) que el código de esta entrada, pero su filosofía es similar:

1) se precalcula la transformada de un corpus de secuencias (y se comprime)
2) se buscan patrones exactos/inexactos sobre el corpus precalculado (comprimido)

Para aprender más, además de la lectura del artículo original de Burrows y Wheeler, un material recomendable para iniciarse en estos temas puede ser el libro de Gabriel Valiente  'Combinatorial pattern matching algorithms in computational biology using Perl and R'.

El siguiente código en Perl muestra un prototipo con estas ideas, empleando para la búsqueda el árbol binario implícito en el vector de sufijos:

 use strict;  
   
 my $VERBOSE = 1;  
 my @ALPHABET = ( '$', 'a' .. 'z' ); # orden lexicografico, $ marca el final   
   
 my ($text,$pattern) = ('mississippi','ssi');           
            
 printf("# original text: '%s'\n\n",$text);  
   
 my ($bwt,$start,$F) = BWT($text,$VERBOSE);   
   
 printf("# BWT('%s') = ('%s',%d) | F: %s\n",$text,$bwt,$start,$F);   
   
 # ahora se pueden comprimir $bwt y $F, por ejemplo, con run-length coder:  
 # 'mmmmmsssiiii' (12 bytes) quedaría como 'm5s3i4' (6 bytes)  
   
 my ($orig_text,$ref_array_W) = revBWT($bwt,$start,$VERBOSE);  
   
 printf("# BWTinv('%s',%d) = '%s'\n\n",$bwt,$start,$orig_text);   
   
 print "# searching for '$pattern' matches in '$text' [base 0]:\n";  
 match_pattern_BWT($pattern,$F,$start,$ref_array_W);  
   
   
   
 # Calcula un vector de sufijos de manera que las diferentes subcadenas   
 # de una secuencia queden preordenadas lexicograficamente  
 sub make_suffix_array  
 {  
   my ($seq,$verbose) = @_;       
   
   my $l = length($seq);  
   my @suffix = (0 .. $l-1); # en base 0  
   
   # ordena los length($seq) posibles sufijos lexicográficamente (cmp),  
   # este paso consume RAM y es uno de los pasos que se puede optimizar  
   @suffix = sort {substr($seq,$a) cmp substr($seq,$b)} (@suffix);  
   if($verbose)  
   {  
     print "# suffix array for $seq :\n";  
     foreach my $suf (@suffix)  
     {   
       printf("%3d %s%s\n",  
         $suf,substr($seq,$suf),substr($seq,0,$l-($l-$suf)));   
     }  
     print "\n";  
   }  
   return @suffix;  
 }  
   
 # Devuelve: 1) string BWT, 2) índice $start, que indica la fila del vector  
 # de sufijos que contiene la cadena original, y 3) string $F, la columna   
 # F1..Fn concatenada.  
 # BWT realmente es la columna de caracteres que preceden a la primera (F)  
 # de un vector de sufijos, o la última (L) si en vez de sufijos la cadena  
 # se permuta circularmente  
 # F1 .............$ L1  
 # F2 ............$. L2  
 # Fn ......$....... Ln  
 sub BWT  
 {  
   my ($seq,$verbose) = @_;  
       
     $seq = $seq . '$'; # marca el final   
       
   my @suffix_rotation = make_suffix_array($seq,$verbose);  
   my ($bwt,$start,$F) = ('',0);  
   foreach my $c (0 .. length($seq)-1)  
   {  
     $bwt .= substr($seq,$suffix_rotation[$c]-1,1);  
         $F .= substr($seq,$suffix_rotation[$c],1);  
     if($suffix_rotation[$c] == 0)  
     {  
       # fila del vector de sufijos que contiene cadena original    
       $start = $c;   
     }  
   }  
   return ($bwt,$start,$F);  
 }  
   
 # Devuelve: 1) la cadena original correspondiente a una cadena  
 # previamente transformada con BWT y 2) una referencia a un array @W  
 # que contiene, para cada posicion del vector L (cadena transformada), la que  
 # le sigue en la secuencia original; ojo, el orden de @W se conserva en @F  
 sub revBWT  
 {  
   my ($bwt,$start,$verbose) = @_;  
   my (@L,@C,@K,@M,@W,@T,$original_string);  
   my ($c,$l,$r,$total,$last_char);  
     my $last_alphabet = scalar(@ALPHABET)-1;  
   
   # convierte cadena $bwt en array @L  
   @L = split(//,$bwt);  
     $last_char = scalar(@L)-1;  
   
   # calcula vectores auxiliares C y K, contienen las observaciones   
   # de cada letra del alfabeto en $bwt (@L)  
   foreach $l (0 .. $last_alphabet){ $K[$l] = 0 }  
       
   foreach $c (0 .. $last_char)  
   {   
     foreach $l (0 .. $last_alphabet)  
     {  
       if($L[$c] eq $ALPHABET[$l])  
       {  
         $C[$c] = $K[$l];    
         $K[$l]++;  
       }   
     }   
   }  
       
     foreach $l (0 .. $last_alphabet)  
   {  
     if($verbose && $bwt =~ /$ALPHABET[$l]/)  
     {  
       print "# $ALPHABET[$l] K=$K[$l]\n";  
     }               
   }  
   
   # calcula vector auxiliar M, que contiene la primera aparición  
   # de cada letra del alfabeto en vector implícito F, que es   
   # la primera columna del suffix array subyacente  
   $total = 0;  
   foreach $l (0 .. $last_alphabet)  
   {  
     $M[$l] = $total;  
     $total += $K[$l];    
     if($verbose && $bwt =~ /$ALPHABET[$l]/)      
     {  
       print "# $ALPHABET[$l] M=$M[$l]\n" if($verbose);  
     }  
   }  
   
   # calcula vector decodificador @W (forward), destruye @M  
   foreach $c (0 .. $last_char)  
   {  
     foreach $l (0 .. $last_alphabet)  
     {  
       if($L[$c] eq $ALPHABET[$l])  
       {  
         $W[$M[$l]] = $c;     
         $M[$l]++;   
       }    
     }    
   }  
   
   # decodifica cadena de texto original (@T)  
   $original_string = '';  
   $l = $start;  
   foreach $c (0 .. $last_char)  
   {  
     $original_string .= $L[$l] if($L[$l] ne '$');  
     $l = $W[$l]  
   }  
   
   return ($original_string,\@W);  
 }  
   
 # Funcion equivalente al operador cmp para comparar una cadena $pattern   
 # con la BWT de un texto, dada una fila $row del suffix_array.  
 # Devuelve < 0 si $pattern < $bwt (orden lexicográfico), > 0 en caso contrario,  
 # y 0 si la fila $r comienza con $pattern   
 sub BWTstrcmp  
 {  
   my ($pattern,$F,$row,$ref_W) = @_;  
     
   my @F = split(//,$F);  
   my $m = length($pattern);   
   my ($c,$match) = (0,'');    
       
   while($m >= 0 and $F[$row] eq substr($pattern,$c,1))  
   {  
         $match .= $F[$row];  
     $row = $ref_W->[$row];  
     $m--;  
     $c++;  
   }    
   
   if($m == 0){return 0 } # match  
   else{ return substr($pattern,$c,1) cmp $F[$row] }  
 }  
   
 # Procedimiento que busca ocurrencias de un patron recorriendo el vector de sufijos  
 # implicito en la BWT, por medio del vector @F. Imprime las posiciones encontradas  
 sub match_pattern_BWT  
 {  
   my ($pattern,$F,$start,$ref_W) = @_;  
     
     my ($med,$submed,$low,$high,$comp);  
     my $last = length($F);  
   $low = 0;  
   $high = $last;  
   
   while ($low <= $high)   
   {  
     $med = int (($low+$high)/2);     
     $comp = BWTstrcmp($pattern,$F,$med,$ref_W);  
   
     if($comp < 0){ $high = $med-1 } # retrocedemos   
     elsif($comp > 0){ $low = $med+1 } # avanzamos  
     else   
     {  
       $submed = $med; # filas anteriores  
             while($submed > 0 &&   
         BWTstrcmp($pattern,$F,$submed,$ref_W) == 0)  
       {   
         printf("# match at position %d\n",  
           get_real_position($submed,$start,$ref_W));  
         $submed--;   
       }  
         
             $submed = $med + 1; # filas posteriores  
       while ($submed < $last &&   
         BWTstrcmp($pattern,$F,$submed,$ref_W) == 0)  
       {  
         printf("# match at position %d\n",  
           get_real_position($submed,$start,$ref_W));  
         $submed++;   
       }  
       last;  
     }  
   }  
   
 }  
   
 # Devuelve la posicion en el texto original de una fila de @F  
 sub get_real_position  
 {  
   my ($Fpos,$start,$ref_W) = @_;   
       
   my $real_pos = 0;  
   while($start != $Fpos)  
   {  
     $start = $ref_W->[$start];  
     $real_pos++;       
   }  
   return $real_pos;  
 }  
   

conferencia 'Brazilian genomics in a nutshell'

Hola,
hoy aprovecho para anunciar la conferencia que la Dra. Ana Tereza de Vasconcelos dará el jueves en Zaragoza con el título de 'Brazilian genomics in a nutshell'. Ana Tereza es la directora del grupo de Bioinformática del brasileño laboratorio nacional de computación científica y vendrá a contarnos detalles sobre diversos proyectos de secuenciación en marcha en Brasil y sobre el proceso de anotación en el que su laboratorio está directamente implicado.

(fuente: http://genometools.org)

En qué consiste anotar genomas? Pues es, a grandes rasgos, el proceso por el cual se asigna información biológica a las secuencias de nucleótidos que conforman un genoma, como explica la wikipedia. Por ejemplo, la anotación incluye la tarea de reconocer donde empiezan y acaban cada uno de los genes de un genoma, una tarea nada sencilla, sobre todo en eucariotas.

La charla tendrá lugar el 30 de Junio a las 18H00 en Ibercaja Zentrum, situado en la plaza de los Sitios, 2, en el centro de Zaragoza, y será en portuñol.

Os dejo unos enlaces a artículos representativos del trabajo de Ana Tereza:

 genómica bacteriana y metagenómica

• Swine and poultry pathogens: the complete genome sequences of two strains of Mycoplasma hyopneumoniae and a strain of Mycoplasma synoviae 
•  A System for Automated Bacterial (genome) Integrated Annotation—SABIA
•  Impact of Transcription Units rearrangement on the evolution of the regulatory network of gamma-proteobacteria
•  Genetic characterization of Chromobacterium isolates from black water environments in the Brazilian Amazon 

taxonomía molecular bacteriana

•  Genomic and evolutionary comparisons of diazotrophic and pathogenic bacteria of the order Rhizobiales
•  Genomic taxonomy of Vibrios 

genómica de eucariotas y cáncer

•  Distinct patterns of somatic alterations in a lymphoblastoid and a tumor genome derived from the same individual
•  Combination of measures distinguishes pre-miRNAs from other stem-loops in the genome of the newly sequenced Anopheles darlingi 

bioinformática estructural

•  Structural descriptor database: a new tool for sequence-based functional site prediction
•  Prediction of TF target sites based on atomistic models of protein-DNA complexes

Un saludo,
Bruno

PD: Esta charla se anuncia también en el blog "Te ayudamos a crecer"

Creando un cluster Rocks virtual

Hola,
mientras actualizaba el material didáctico de Programación en clusters Rocks a la versión actual de Linux Rocks (5.4) me pareció buena idea ir probando el código en un cluster virtual, creado como máquinas virtuales corriendo bajo VirtualBox.

(fuente: http://www.rocksclusters.org/rocks-documentation/4.1/getting-started.html)

La ventaja que tiene hacerlo virtual es que cualquiera con unos cuantos Gb de disco libres y más de 2 Gb de RAM puede probar a programar en un cluster y así aprender a manejarse en ese entorno antes de tener acceso a uno real.

Éstos son los pasos necesarios, probados en mi Ubuntu 10.04 LTS y con la versión 4.0.8 de VirtualBox:

• Instala VirtualBox para tu sistema Linux (puedes descargarlo de
  http://www.virtualbox.org/wiki/Downloads)
  
  
• Descarga e instala el 'VirtualBox Extension Pack', que puedes obtener en el mismo lugar
  
  
• Descarga ISO del DVD Jumbo de Rocks, que ya incluye los rolls esenciales
  (puedes descargarlo de http://www.rocksclusters.org,
  comprobando tras la descarga la suma MD5). En este tutorial usaremos la versión 5.4, que se llama area51+base+bio+ganglia+hpc+kernel+os+sge+web-server+xen-16.12.2009-15.16.53.i386.disk1.iso,
  que guardamos en la carpeta /home/pepe/soft/
  
  
• Elige un directorio donde ubicar los discos duros virtuales, como por ejemplo /home/pepe/
  
  
• Crea el nodo principal tecleando (copia y pega) en el terminal estos comandos:
  

   VBoxManage createvm --name "vRocks54-Frontend" --ostype RedHat --register  
     
   VBoxManage modifyvm "vRocks54-Frontend" --memory 1000 --vram 32 --nic1 intnet --nic2 nat --audio none --nictype1 82540EM --nictype2 82540EM --boot1 dvd --boot2 disk --boot3 none --boot4 none  
     
   VBoxManage createhd --filename "/home/pepe/vRocks54-Frontend.vdi" --size 50000 --variant Standard   
     
   VBoxManage storagectl "vRocks54-Frontend" --name "SATA Controller" --add sata --controller IntelAhci  
     
   VBoxManage storageattach "vRocks54-Frontend" --storagectl "SATA Controller" --port 0 --device 0 --type hdd --medium "/home/pepe/vRocks54-Frontend.vdi"  
     
   VBoxManage storagectl "vRocks54-Frontend" --name "IDEcontrol" --add ide  
     
   VBoxManage storageattach "vRocks54-Frontend" --storagectl "IDEcontrol" --port 0 --device 0 --type dvddrive --medium /home/pepe/soft/area51+base+bio+ganglia+hpc+kernel+os+sge+web-server+xen-16.12.2009-15.16.53.i386.disk1.iso  
     
   VBoxManage startvm "vRocks54-Frontend"  
  

  
  Tras estos comandos da comienzo el proceso normal de instalación del nodo maestro o frontend , documentado en la web de Rocks. El proceso es casi automático, pero debes recordar dos cosas:
  
  
  
  • escribe frontend cuando aparezca el terminal de arranque boot:
    
  • cuando llegue el momento elige como mínimo estos rolls : kernel, base, web server, os
    
  
  
• Una vez instalado el frontend puedes añadir el paquete
  VirtualBox guest additions para mayor comodidad
  
  
• En el terminal del nodo maestro teclea como superusuario $ insert-ethers
  
  
• Para cada nodo compute que quieras crear, empezando por el 0-0, hasta el 0-N teclea en el terminal de tu Linux:
  
  

   VBoxManage createvm --name "vRocks54-Compute-0-0" --ostype RedHat --register  
     
   VBoxManage modifyvm "vRocks54-Compute-0-0" --memory 1000 --vram 32 --nic1 intnet --audio none --nictype1 82540EM --boot1 net --boot2 disk --boot3 none --boot4 none  
     
   VBoxManage createhd --filename "/home/pepe/vRocks54-Compute-0-0.vdi" --size 50000 --variant Standard   
     
   VBoxManage storagectl "vRocks54-Compute-0-0" --name "SATA Controller" --add sata --controller IntelAhci  
     
   VBoxManage storageattach "vRocks54-Compute-0-0" --storagectl "SATA Controller" --port 0 --device 0 --type hdd --medium "/home/pepe/vRocks54-Compute-0-0.vdi"  
     
   VBoxManage storagectl "vRocks54-Compute-0-0" --name "IDE Controller" --add ide   
     
   VBoxManage startvm "vRocks54-Compute-0-0"  
  

  
  
• En el terminal del nodo maestro teclea como superusuario la tecla F8 para terminar de añadir nodos
  

Y ésto es todo, con esto tienes ya un cluster en miniatura para probar tus programas, que puedes borrar con dos golpes de ratón desde el VirtualBox,
un saludo,
Bruno

Aprender a divulgar la ciencia y no morir en el intento

La Universidad de la Rioja ha editado online el "Manual de Comunicación para Investigadores" que consiste en un compendio de consejos y ejemplos de comunicación y divulgación científica.


A la hora de divulgar la ciencia que se hace día a día en los laboratorios, los investigadores nos hacemos muchas preguntas: ¿Por qué? ¿Dónde? ¿Cómo? ¿Para quién?... estas preguntas y muchas otras son respondidas en este manual con breves explicaciones y ejemplos sobre cada uno de los problemas y situaciones más habituales en la divulgación de la investigación.

Respecto al primer interrogante que surge: ¿Por qué divulgar la investigación? Me ha gustado mucho una de las contestaciones que se ofrecen en el manual:
"Son los contribuyentes quienes pagan las facturas de los investigadores y de la actividad científica. Por ello, los científicos, además de corresponder haciendo ciencia de la mejor calidad, deben promover la difusión de sus resultados a todos los colectivos sociales." Luis Martínez Saez (profesor del Máster en Comunicación Científica, Médica y Ambiental de la Universidad Pompeu Fabra).

Ejemplos cercanos de divulgación en prensa y televisión: 

Secuencia de Escherichia coli O104:H4

Hola,
en las noticias llevamos una semana escuchando hablar de la que parece una nueva cepa patógena de Escherichia coli, la cepa  enterohemorrágica O104:H4, que ha causado ya varias muertes en Alemania. En los medios ha habido cierta confusión a la hora de divulgar sobre este tema, pero en Internet ya es posible encontrar buenas referencias, como por ejemplo este blog al que me redirigió mi colega Pablo Vinuesa.

Árbol pangenómico de E.coli y especies relacionadas, basado en su contenido génico. Tomado de http://www.springerlink.com/content/kk08p747348hq301/fulltext.html.

Como ejemplo de lo rápido que han ido la cosas me gustaría recordar que las primeras noticias alertando sobre esta nueva cepa aparecieron en la última semana de Mayo y que el día 2 de Junio ya estaba disponible para descarga una primera versión de su secuencia genómica ensamblada, obtenida en el BGI (antiguo Beijing Genomics Institute). Ahora mismo hay disponibles un archivo con los contigs encontrados y las lecturas originales en formato FASTQ (aquí tienes un ejemplo), pero mr imagino que se irán actualizando los contenidos en los próximos días. Tal como se relata aquí y aquí, la secuencia ensamblada y experimentos de multilocus sequence typing sugieren que efectivamente se trata de una cepa nueva,
un saludo,
Bruno

PD Iré anadiendo aquí más enlaces sobre este tema:
1) http://scienceblogs.com/mikethemadbiologist/2011/06/i_dont_think_the_german_e_coli.php
2)  https://github.com/ehec-outbreak-crowdsourced/BGI-data-analysis/wiki
3) http://eterno-bucle.blogspot.com/2011/06/ensamble-de-cepa-de-e-coli-que-esta.html
4)  http://enews.patricbrc.org/1172/e-coli-outbreak-new-comprehensive-comparisons

Bio::Phylo, una librería para análisis filoinformático

Hola,
hace unos días me enteré por mi colega Pablo Vinuesa de la publicación de Bio::Phylo, una colección de 59 módulos escritos en Perl para análisis filogenético (filoinformático, como dicen sus autores). El artículo original aparece en  la revista BMC Bioinformatics e incluye ejemplos de como utilizar esta librería para tareas típicas pero complejas, como por ejemplo el análisis automático (y el almacenamiento persistente como NeXML) de las anotaciones de árboles en formato Newick.
 

Lo más llamativo de Bio::Phylo es su facilidad de instalación, ya que por diseño no tiene dependencias externas al core de Perl, y su extensa documentación en CPAN, lo que lo convierte en una buena opción para este tipo de aplicaciones, superando con creces las prestaciones de módulos como Bio::TreeIO (ver por ejemplo http://www.eead.csic.es/compbio/material/bioperl/node45.html ). En Ubuntu se instala en el terminal en menos de un minuto con la instrucción:
$ sudo cpan -i Bio::Phylo

Un saludo,
Bruno

Guardar en un archivo los contenidos de una variable de perl para reusarlos

Cuando obtenemos resultados de costosos cálculos computacionales los guardamos para poder volver a usarlos sin tener que calcularlos de nuevo.

Normalmente se guardan en archivos de texto con el formato deseado. Pero hay veces que nos interesa guardar los resultados tal y como los guardamos temporalmente en una variable de perl mientras son procesados.

Guardar una variable intermedia de un script de perl en un archivo puede tener las siguientes utilidades:

• No perder ninguna información del procesado de los datos.
• Poder importar la variable rápidamente al volver a correr el script.
• Revisar los datos contenidos por la variable de una forma sencilla.

El código de perl escrito a continuación usa el módulo Data::Dumper de perl para guardar los contenidos de una variable compleja en el archivo 'backup.txt' mediante la subrutina 'store_data' y posteriormente recupera los contenidos de la variable con la subrutina 'recover_data'.

   
 use Data::Dumper;  
   
 my $data = {  
        'Research Centers' => {  
               'EEAD-CSIC' => [  
                       {  
                        'Group' => 'Laboratory of Computational Biology',  
                        'Address' => 'Av. Montañana 1.005, 50059 Zaragoza (Spain)',  
                        'Phone' => '+34 976716089'  
                       },  
                       {  
                        'Group' => 'Genética y Desarrollo de Materiales Vegetales',  
                        'Address' => 'Av. Montañana 1.005, 50059 Zaragoza (Spain)',  
                        'Phone' => '+34 976716079'  
                       }  
                      ]  
              }  
       };  
   
 my $file = 'backup.txt';  
   
 store_data($data, $file);  
   
 my $recovered_data = recover_data($file);  
   
 print $recovered_data->{'Research Centers'}{'EEAD-CSIC'}[0]{'Group'}."\n";  
 print $recovered_data->{'Research Centers'}{'EEAD-CSIC'}[0]{'Address'}."\n";  
 print $recovered_data->{'Research Centers'}{'EEAD-CSIC'}[0]{'Phone'}."\n";  
   
   
 # Stores data from a variable into a file  
 sub store_data {  
   
      my ($data, $file) = @_;  
   
      open(FILE, ">$file");  
      print FILE Dumper($data);  
      close FILE;  
   
 }  
   
 # Recovers data from a file into a variable  
 sub recover_data {  
   
      my ($file) = @_;  
   
      return do $file;  
   
 }  

Algoritmos en Bioinformática Estructural, versión 2011

Buenas,
con motivo de las clases que daré la semana que viene en la Licenciatura en Ciencias Genómicas de la UNAM, en México, he actualizado el material sobre 'Algoritmos en Bioinformática Estructural', disponible en la URL:
http://www.eead.csic.es/compbio/material/algoritmos3D

Los cambios incluyen:

1. actualización de la sección sobre estabilidad de un dúplex de DNA
2. renovación de la teoría sobre diseño de primers
3. he añadido una nueva sección sobre modelado de interfaces por homología
4. un montón de nuevas referencias a lo largo de todo el material, con sus URLs para facilitar su consulta
5. nuevos ejercicios

Espero que el material sea de provecho y, por favor, si encontráis errores, hacedmelos llegar,
un saludo,
Bruno

Paralelizar procesos en perl con threads

Con los modernos ordenadores que tenemos con procesadores de 2, 4, 8 o más núcleos podemos pensar en paralelizar procesos fácilmente en Perl. ¿Qué significa esto? Básicamente que el tiempo que tardará en correr nuestro programa se reducirá proporcionalmente al número de núcleos, si a cada núcleo le mandamos tareas diferentes a realizar (threads)

Pero hay que ser cuidadoso, no todos los programas son adecuados para ser paralelizados. Un programa paralelizable en threads requiere que sus procesos puedan ser independientes y una vez finalizados podamos recoger y unificar sus resultados sin problemas. Además paralelizar significa multiplicar las copias en memoria de las variables, por lo cual sólo podremos paralelizar procesos que no consuman grandes cantidades de memoria o el proceso de copia ralentizará los cálculos.

El código mostrado a continuación realiza una suma de números de tres formas diferentes:

1. De forma tradicional, con 3 procesos de suma en serie.
2. Usando 3 threads que calculan simultáneamente.
3. Usando un bucle con 3 threads que calculan simultáneamente.

La conclusión que podemos obtener de los resultados es que el proceso se realiza casi 3 veces más rápido que de la forma tradicional mediante el método de 3 threads. Pero que cuando se realiza la paralelización para cálculos más sencillos (bucle de threads), no sólo no se mejora el tiempo de cálculo, sino que se ralentiza enormemente el cálculo (5 veces más que el tradicional), esto se debe a que el coste de lanzar los threads no se compensa con el tiempo de cálculo de cada uno.

Resultado:

 Total = 300000000  
 Time taken by the traditional way (3 serial processes) was 26 wallclock secs (26.48 usr 0.00 sys + 0.00 cusr 0.00 csys = 26.48 CPU) seconds  
   
 Total = 300000000  
 Time taken by 3 parallel threads was 9 wallclock secs (25.94 usr 0.00 sys + 0.00 cusr 0.00 csys = 25.94 CPU) seconds  
   
 Total = 300000000  
 Time taken by a loop of 10000 times 3 threads was 104 wallclock secs (103.41 usr 1.01 sys + 0.00 cusr 0.00 csys = 104.42 CPU) seconds  
   

Código:

 use threads;  
 use Benchmark;  
 use Config;  
 $Config{useithreads} or die('Recompile Perl with threads to run this program.');  
   
 # Subroutine to test multithreading  
 sub calc {  
      my ($number) = @_;  
      my $i=0;  
      while ($i<$number) {  
           $i++;  
      }  
      return $i;  
 }  
   
 # Define a number to calculate  
 my $number = 100000000;  
   
 # Start timer  
 my $start = new Benchmark;  
   
 # Run subroutines in the traditional way  
 my $res1 = calc($number);  
 my $res2 = calc($number);  
 my $res3 = calc($number);  
 my $total = $res1 + $res2 + $res3;  
   
 # End timer  
 my $end = new Benchmark;  
   
 # Calculate difference of times  
 my $diff = timediff($end, $start);  
   
 # Print final result  
 print "\nTotal = $total\n";  
   
 # Report benchmark  
 print "Time taken by the traditional way (3 serial processes) was ", timestr($diff, 'all'), " seconds\n\n";  
   
 # Start timer  
 $start = new Benchmark;  
   
 # Create multiple threads running each one the subroutine  
 my $thr1 = threads->create(\&calc, $number);  
 my $thr2 = threads->create(\&calc, $number);  
 my $thr3 = threads->create(\&calc, $number);  
   
 # Check subroutines ending and retrieve results   
 $res1 = $thr1->join();  
 $res2 = $thr2->join();  
 $res3 = $thr3->join();  
 $total = $res1 + $res2 + $res3;  
   
 # End timer  
 $end = new Benchmark;  
   
 # Calculate difference of times  
 $diff = timediff($end, $start);  
   
 # Print final result  
 print "Total = $total\n";  
   
 # Report benchmark  
 print "Time taken by 3 parallel threads was ", timestr($diff, 'all'), " seconds\n\n";  
   
 # Start timer  
 $start = new Benchmark;  
   
 # Divide the process  
 $total = 0;  
 my $divide = 10000;  
 $number = $number / $divide;  
   
 # Create multiple threads running each one the subroutine  
 for (my $i=0; $i<$divide; $i++){  
      $thr1 = threads->create(\&calc, $number);  
      $thr2 = threads->create(\&calc, $number);  
      $thr3 = threads->create(\&calc, $number);  
      # Check subroutines ending and retrieve results   
      $res1 = $thr1->join();  
      $res2 = $thr2->join();  
      $res3 = $thr3->join();  
      $total += $res1 + $res2 + $res3;  
 }  
   
 # End timer  
 $end = new Benchmark;  
   
 # Calculate difference of times  
 $diff = timediff($end, $start);  
   
 # Print final result  
 print "Total = $total\n";  
   
 # Report benchmark  
 print "Time taken by a loop of 10000 times 3 threads was ", timestr($diff, 'all'), " seconds\n\n";