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

Extraer coordenadas de átomos en un fichero PDB

Los ficheros PDB (Protein Data Bank) contienen las coordenadas espaciales de los átomos de proteínas cuya estructura está resuelta por técnicas de rayos X o resonancia magnética nuclear (RMN). En una entrada de blog de Bruno podéis encontrar una muy buena explicación de estas técnicas realizada por el periódico El País.

Las estructuras de proteínas en formato PDB contienen mucha información que habitualmente no nos interesa para nuestros experimentos e interfiere con muchos programas que únicamente necesitan las coordenadas espaciales de los átomos. Además un fichero PDB normalmente contiene datos de múltiples moléculas y en diferentes posiciones en el cristal. Por ello es muy conveniente extraer de los ficheros PDB únicamente la información estructural que vamos a utilizar y descartar el resto.

En el post anterior vimos como extraer líneas de un texto o archivo. En la entrada de hoy veremos como reusar ese código para extraer los átomos deseados de un fichero PDB.

En el ejemplo se indica el fichero PDB original en la variable $pdb_file, en nuestro caso será la estructura de la famosa insulina que se puede descargar aquí (2INS.pdb). Y en el array @desired_coords le indicaremos el tipo de átomos que queremos extraer en un formato muy similar al usado por PyMOL.

   
 # File with 3D coords of insulin  
 my $pdb_file = "2INS.pdb";  
   
 # Extract coords of the alpha carbons of the chain A, and alpha and beta carbons of histidines of chain B  
 my @desired_coords = ('a/*/ca/', 'b/*/ca|cb/his');  
   
 open(PDBFILE,$pdb_file);  
 my $pdb_content = join('',<PDBFILE>);  
 close PDBFILE;  
   
 my $pdb_coords = join('',extract_pdb_coords($pdb_content,\@desired_coords));  
   
 open(PDBFILE,">$pdb_file.out");  
 print PDBFILE $pdb_coords;  
 close PDBFILE;  
   
 # Extract desired PDB coordinates from PDB entries  
 sub extract_pdb_coords {  
   
      my ($pdb_content, $types) = @_;  
   
      my $pdb_data;  
      my $patterns;  
      my @pattern_parts = ('chain','res_number','res_type','res_name');  
      foreach my $type (@{$types}) {  
           my $count = 0;  
           my %pattern;  
           while ($type =~ /([^\/]+)/g){  
                $type = $'; # Text to the right of the match  
                $pattern{$pattern_parts[$count]} = $1;  
                $count++;  
           }  
           push(@{$patterns},\%pattern);  
      }  
   
      foreach my $pattern (@{$patterns}){  
           my ($chain,$res_number,$res_type,$res_name);  
           if (!defined($pattern->{'chain'}) || !$pattern->{'chain'} || $pattern->{'chain'} eq '*'){  
                $chain = '\w{1}';  
           } else {  
                $chain = uc($pattern->{'chain'});  
           }  
           $pattern->{'res_number'} =~ s/\s+//;  
           if (!defined($pattern->{'res_number'}) || !$pattern->{'res_number'} || $pattern->{'res_number'} eq '*'){  
                $res_number = '\d+';  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)$/) {  
                $res_number = $1;  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)-(\d+)$/) {  
                $res_number = '('.join('|', $1 .. $2).')';  
           } elsif ($pattern->{'res_number'} =~ /\d,/) {  
                $res_number = '('.join('|', split(",",$pattern->{'res_number'})).')';  
           }  
           if (!defined($pattern->{'res_type'}) || !$pattern->{'res_type'} || $pattern->{'res_type'} eq '*'){  
                $res_type = '[\w\d]+';  
           } else {  
                $res_type = uc($pattern->{'res_type'});  
           }  
           if (!defined($pattern->{'res_name'}) || !$pattern->{'res_name'} || $pattern->{'res_name'} eq '*'){  
                $res_name = '\w{3}';  
           } else {  
                $res_name = uc($pattern->{'res_name'});  
           }  
           $pattern = '/(ATOM|HETATM)\s+\d+\s+('.$res_type.')\s+('.$res_name.')\s('.$chain.')\s+('.$res_number.')\s+.+/';  
      }  
   
      my @pdb_data_lines = extract_lines_from_text($pdb_content, $patterns);  
      if (@pdb_data_lines){  
           $pdb_data = join("\n",@pdb_data_lines);  
      }  
   
      return $pdb_data;  
   
 }  
   
 # Extract lines from text with the desired patterns  
 sub extract_lines_from_text {  
   
      my ($text, $patterns) = @_;  
   
      my @data;  
      my @lines = split("\n",$text);  
   
      foreach my $line (@lines){  
           foreach my $pattern (@{$patterns}){  
                if ($pattern =~ /^\/(.+)\/$/){  
                     if ($line =~ /$1/){  
                          push(@data,$line);  
                          last;  
                     }  
                } else {  
                     if ($line =~ /\Q$pattern\E/){  
                          push(@data,$line);  
                          last;  
                     }  
                }  
           }  
      }  
   
      return @data;  
   
 }