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

curso APIs de Ensembl

Hola,
en línea con la anterior entrada del blog os paso este anuncio de Miguel Pérez-Enciso (miguel.perez@uab.es ) :

Del 4 al 6 de octubre organizo el siguiente curso en el CRAG, más detalles + adelante. Lo darán dos personas del EBI. Se requiere conocimientos sólidos de perl. Máximo 20 personas.

"Ensembl uses MySQL relational databases to store its information. A comprehensive set of Application Programme Interfaces (APIs) serve as a middle-layer between underlying database schemes and more specific application programmes. The APIs aim to encapsulate the database layout by providing efficient high-level access to data tables and isolate applications from data layout changes.

This 2-day workshop is aimed at developers interested in exploring Ensembl beyond the website. Participants will be expected to have experience in writing Perl programs. A background in object oriented programming techniques and familiarity with databases (MySQL) would be an advantage.

The workshop covers the following four databases and APIs. For each of them the database schema and the API design as well as its most important objects and their methods will be presented. This will be followed bypractical sessions in which the participants can put the learned into practice by writing their own Perl scripts.

**Ensembl Core databases and API

The set of species-specific Ensembl Core databases stores genome sequences and most of the annotation information. This includes the gene, transcript and protein models annotated by the Ensembl automated genome analysis and annotation pipeline. Ensembl Core databases also store assembly information, cDNA and protein alignments, external references, markers and repeat regions data sets.

**Ensembl Variation databases and API

The large amount of genetic variation information is organised in a set of species-specific Ensembl Variation databases. This includes basic variation information (alleles, flanking sequences, genotypes), the functional consequence of a variation and more complex data such as LD values from genotype data and tagged SNPs in a population.

**Ensembl Compara database and API

The Ensembl Compara multi-species database stores the results of genome-wide species comparisons re-calculated for each release. The comparative genomics set includes pairwise whole genome alignments and synteny regions. The comparative proteomics data set contains orthologue predictions and protein family clusters. "

Bajar de ensembl todas las secuencias de proteína de los ortólogos de un gen (o de varios)

Dada una secuencia de proteína, puede interesarnos bajar de una sola vez todas las secuencias de proteínas ortólogas conocidas. Si aceptamos el criterio de ensembl para determinar qué es y qué no es ortólogo (una buena opción en una decisión a veces complicada), las siguientes instrucciones pueden servir de guía para esta tarea.

i) pasos previos

Tendréis que tener instalado el ensemble-api (ver instrucciones de instalación) y el módulo DBD::Mysql. En Mac OS X he conseguido instalar este módulo usando macports. En 5 líneas, sería:

~% sudo port selfupdate

~% sudo port install mysql5-server

~% sudo port install p5-dbd-mysql

~% sudo port install p5-dbi

~% sudo port -f activate p5-dbi

ii) bajar los ortólogos

En un fichero de texto, ponéis una lista con los identificadores de genes para los que queréis hacer la búsqueda.

(por ejemplo)

~% cat IDs.txt

ENSMUSG00000031628

ENSMUSG00000027997

ENSG00000003400

Ejecutáis get_orthologues.pl dándole este fichero como $ARGV[0]

~% perl get_orthologues.pl IDs.txt

Si veis algo así en pantalla, todo OK.

.

.

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000031628

main IDs for each orthologue written to ENSMUSG00000031628.mainIDs.txt

all IDs for each orthologue written to ENSMUSG00000031628.allIDs.txt

sequeences for the main IDs dumped to ENSMUSG00000031628.ort.fa

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000027997

main IDs for each orthologue written to ENSMUSG00000027997.mainIDs.txt

all IDs for each orthologue written to ENSMUSG00000027997.allIDs.txt

sequeences for the main IDs dumped to ENSMUSG00000027997.ort.fa

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSG00000003400

main IDs for each orthologue written to ENSG00000003400.mainIDs.txt

all IDs for each orthologue written to ENSG00000003400.allIDs.txt

sequeences for the main IDs dumped to ENSG00000003400.ort.fa

.

.

(en mi ordenador, para estos 3 IDs tarda unos 5 minutos)

El resultado, para cada identificador ENSG.., son como véis 3 ficheros:

ENSG*.mainIDs.txt --> listado con los protein-IDs del transcript principal de cada ortólogo

ENSG*.allIDs.txt --> listado con los protein-IDs de todas las secuencias que ensembl guarda para cada ortólogo

ENSG*.ort.fa --> secuencias de proteina (formato fasta) correspondientes al transcript principal de cada ortólogo

Aclaración por si alguien se ha liado con esto de transcript_principal / todas_las_secuencias. Como sabéis, para cada gen (y para cada uno de sus ortólogos) hay descritos varios transcripts normalmente. ¿Cuál utiliza ensembl como transcript principal? El más largo. La secuencia de proteína de ese se graba en OUT3, su ID en OUT1 y las IDs de todos en OUT2.

El código de get_orthologues.pl es el siguiente:

#!/usr/bin/perl

# Script to retrieve protein sequences for 1-to-1(and 1-to-many)-orthologues for a given gene.
# (adaptado del que me pasó Bert Overduin de helpdesk@enseml.org)

# Definicion de librerias y modulos

use strict;
use warnings;

use lib '/Users/jramon/Soft/ensembl_API/ensembl/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-compara/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-variation/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-functgenomics/modules/';
use lib '/Users/jramon/Soft/BioPerl-1.6.1/';
use Bio::EnsEMBL::Registry;

my $reg = "Bio::EnsEMBL::Registry";
$reg->load_registry_from_url('mysql://anonymous@ensembldb.ensembl.org');

my $member_adaptor = $reg->get_adaptor("Multi", "compara", "Member");
my $homology_adaptor = $reg->get_adaptor("Multi", "compara", "Homology");


# Leer identificadores de genes y buscar ortologos para cada uno
# guardando IDs en OUT1,2 y secuencias en OUT3

open (INPUT_IDs,"$ARGV[0]");
while () {
my $gene_stable_id = $_;
chop $gene_stable_id;

my $member = $member_adaptor->fetch_by_source_stable_id("ENSEMBLGENE", $gene_stable_id);

open (OUT1,">$gene_stable_id.mainIDs.txt");
open (OUT2,">$gene_stable_id.allIDs.txt");
open (OUT3,">$gene_stable_id.ort.fa");
print "--\nObtaining 1-to-1 and 1-to-many orthologues for $gene_stable_id\n";
print "main IDs for each orthologue written to $gene_stable_id.mainIDs.txt\n";
print "all IDs for each orthologue written to $gene_stable_id.allIDs.txt\n";
print "sequences for the main IDs dumped to $gene_stable_id.ort.fa\n";

my $all_homologies = $homology_adaptor->fetch_all_by_Member($member);

foreach my $homology (@$all_homologies) {
if($homology->description eq (('ortholog_one2one')||('ortholog_one2many'))){
my $members = $homology->get_all_Members();
foreach my $member (@$members) {
 if($member->stable_id ne $gene_stable_id){
  my $canonical_peptide = $member->get_canonical_peptide_Member;
  print OUT1
   $canonical_peptide-> stable_id, "\n";

  my $gene_adaptor = $reg->get_adaptor($member->genome_db->name, "core", "Gene");
  my $gene = $gene_adaptor->fetch_by_stable_id($member->stable_id);
  my @transcripts = @{$gene->get_all_Transcripts};
  foreach my $transcript(@transcripts){
   if($transcript->translation){
     print OUT2 $transcript->translation->stable_id, "\n";
   }
  }

  print OUT3
   ">", $canonical_peptide->stable_id, "\n",
   $canonical_peptide->sequence, "\n";
  }
}
}
}
close OUT1;
close OUT2;
close OUT3;
}
close INPUT_IDs;

De todas formas, para el volumen de datos que ha de bajar, se me hace un poco lento el código. Se agradecen mejoras en este sentido.