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";