Curso 'Next Generation Sequence Analysis: Practice and Departure to New Frontiers'

Hola,
en el 11 GABI Status Meeting, que tuvo lugar recientemente en Alemania,
nos anunciaron el siguiente taller de verano, que tendrá lugar en el mes de Agosto:

'Next Generation Sequence Analysis: Practice and Departure to New Frontiers'

organized by: R. Fries (Animal Breeding, TUM), T. Meitinger (Human Genetics, TUM), K.F.X. Mayer (MIPS, Helmholtz Zentrum München), T. Strom (Human Genetics, Helmholtz Zentrum München)

-> OPEN FOR APPLICATION! (see link below)

Dates:       08.08. – 15.08.2011

Location:  Herrsching am Ammersee, Germany

Description:
The Synbreed Summer School provides an introduction to next generation sequence analysis for PhD students and postdoctoral researchers in animal and plant breeding. The course consists of a practical part with hands-on exercises guided by experienced Synbreed scientists and guest lecturers covering the frontiers of next gen sequencing in animal and plant improvement.

Program Overview:

• Practice the pipeline from next gen sequence reads to annotated variants (alignment, visualization, variant calling, imputing and annotation)
• Visiting the sequencing facility of the Helmholtz Zentrum Munich
• Guest lectures about advances in sequencing technology, genome assembly, genomic selection, next gen population genomics, efficient computation and agricultural genome projects
• Download preliminary programme

Target Group and Requirements:

• The course is directed towards PhD students and postdoctoral researchers in animal and plant breeding
• Knowledge of Linux and a scripting language (e.g. Phyton)
• Dual-core processor equipped laptop with a recent Linux distribution (e.g. Ubuntu 10.10) and at least 150 GB free storage. Linux setup assistance will be provided during the introductory session.

Course language: English

Costs:
Course fee                                    500 Euro    
Participant in single bedroom*         826 Euro
Participant in double bedroom*        756 Euro
(*includes: full accommodation and meals)

Services at Course Location: please visit http://www.hdbl-herrsching.de.

-> For your application click here!  

Application deadline:  01.06.2011

Number of participants: 20

Contact:
Project Coordination Synbreed
Natalie Ohl / Wolf-Christian Saul
Chair of Plant Breeding
Center of Life and Food Sciences Weihenstephan
Technische Universität München
Emil-Ramann-Str. 4
D-85354 Freising

Ph.:    +49 (0)8161/71-5226
Fax:    +49 (0)8161/71-4511
Email: synbreed[at]wzw.tum.de

Vectores de sufijos para buscar patrones (suffix arrays)

Buenas,
hace poco tuve la suerte de cursar un taller práctico sobre problemas en el análisis, mapeo y ensamblado de secuencias genómicas, impartido por Paolo Ribeca (autor del software GEM). A parte de volver un poco abrumado por las dificultades que encierran los datos producidos por estas tecnologías de secuenciación de gran capacidad, sobre todo para bichos que no tienen todavía un genoma de referencia de buena calidad, he aprendido qué estructuras de datos soportan el manejo y análisis de todas estas secuencias. Además de tablas hash y grafos de De Bruijn, y la magia de Burrows-Wheeler, que se merecerá su propia entrada, hemos aprendido las ventajas de usar vectores de sufijos (suffix arrays) para el problema de buscar patrones o subcadenas exactas en textos muy largos, como son las secuencias genómicas. Pero, qué es un vector de sufijos? Trataré de explicarlo con el siguiente código en Perl:

 # 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 @suffix = (0 .. length($seq)-1); # en base 0  
   
    # ordena los length($seq) posibles sufijos lexicogr´aficamente 
    @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\n",$suf,substr($seq,$suf));   
       }  
       print "\n";  
    }  
    return @suffix;  
 }  

Si llamamos a la subrutina

my @suffix = make_suffix_array('TTTTAGATCGATCGACTAGACTACGACTCGA',1);

obtendremos:

30 A
22 ACGACTCGA
19 ACTACGACTCGA
14 ACTAGACTACGACTCGA
25 ACTCGA
17 AGACTACGACTCGA
4  AGATCGATCGACTAGACTACGACTCGA
10 ATCGACTAGACTACGACTCGA
6  ATCGATCGACTAGACTACGACTCGA
28 CGA
12 CGACTAGACTACGACTCGA
23 CGACTCGA
8  CGATCGACTAGACTACGACTCGA
20 CTACGACTCGA
15 CTAGACTACGACTCGA
26 CTCGA
29 GA
18 GACTACGACTCGA
13 GACTAGACTACGACTCGA
24 GACTCGA
9  GATCGACTAGACTACGACTCGA
5  GATCGATCGACTAGACTACGACTCGA
21 TACGACTCGA
16 TAGACTACGACTCGA
3  TAGATCGATCGACTAGACTACGACTCGA
27 TCGA
11 TCGACTAGACTACGACTCGA
7  TCGATCGACTAGACTACGACTCGA
2  TTAGATCGATCGACTAGACTACGACTCGA
1  TTTAGATCGATCGACTAGACTACGACTCGA
0  TTTTAGATCGATCGACTAGACTACGACTCGA

que es una lista de los 31 posibles sufijos de la cadena original, en orden lexicográfico. De hecho, si te fijas en el código, el vector realmente contiene sólo las posiciones (en base 0) de los sufijos ordenados, no su secuencia. Obviamente la construcción de este vector es costosa al necesitar de una ordenación (en Perl es por defecto un mergesort con un caso peor O(NlogN) ), pero luego permite consultas  O(logN), es decir, mucho más rápidas que simplemente recorrer la secuencia de principio a fin, al usar implícitamente un árbol binario.

Puedes probarlo mediante el siguiente ejemplo:

 match_pattern('TAG','TTTTAGATCGATCGACTAGACTACGACTCGA'); 

 # aprovecha el orden del vector de sufijos para buscar cualquier subcadena   
 # o patro´n (pattern) por medio de un árbol binario   
 sub match_pattern   
 {  
    my ($pattern,$seq) = @_;  
    print "# looking for pattern $pattern in sequence $seq (base 0)\n";  
    my @suffix = make_suffix_array($seq,1);  
    my $low = 0;  
    my $high = $#suffix;  
    my $patl = length($pattern);  
    my ($med,$submed);  
    while ($low <= $high)   
    {  
       my $med = int (($low+$high)/2); # punto medio de la búsqueda        
       # comparacion lexicográfica en punto medio, mira 'perldoc perlop'  
       my $comp = $pattern cmp substr($seq,$suffix[$med],length($pattern));  
       if($comp < 0){ $high = $med-1 }  # retrocedemos en @suffix   
       elsif($comp > 0){ $low = $med+1 } # avanzamos  
       else   
       {  
          my $submed = $med - 1; # sufijo inmediatamente anterior al punto medio  
          while($submed > 0 && $pattern eq substr($seq,$suffix[$submed],$patl))  
          {   
             print "# match at position $suffix[$submed]\n";  
             $submed--;   
          }  
          while ($med < $#suffix-1 && $pattern eq substr ($seq,$suffix[$med],$patl))  
          {  
             print "# match at position $suffix[$med]\n";  
             $med++;   
          }  
          last;  
       }  
    }  
 }  

Hasta otra,
Bruno

Convertir un archivo PostScript en PDF con Perl

Normalmente convertimos los archivos PostScript en formato PDF y viceversa abriéndolos con nuestro visor de documentos (por ejemplo Okular para los que somos fans de KDE) y usando el menú de Impresión donde tenemos la opción de imprimir nuestro archivo en formato PostScript o PDF. Los usuarios un poco más expertos en linux prefieren usar el comando ps2pdf (también existe pdf2ps):

Si queremos automatizar esta tarea en nuestros scripts de Perl, existe el módulo PostScript::Convert que nos permitirá hacerlo de una forma sencilla:

 use PostScript::Convert;  
 psconvert($infile, filename => $outfile, format => 'pdf');  

Sin embargo, todos los métodos anteriores generarán PDFs recortados cuando el tamaño del PostScript a convertir no es estándar, por ello os propongo una pequeña subrutina que solucionará este problema:

 # Converts a PostScript file into a PDF  
 sub convert_ps_to_pdf{  
      my ($infile,$outfile) = @_;  
      my ($height,$weight);  
      open(IDENTIFY, "identify $infile |")|| die "# $0 : cannot run 'identify $infile'\n";  
      my $ps_properties = join('',<IDENTIFY>);  
      if ($ps_properties =~ /$infile PS (\d+)x(\d+)/){  
           $weight = $1;  
           $height = $2;  
      } else {  
           die "# convert_ps_to_pdf failed to identify PostScript dimensions.\n";  
      }  
      close IDENTIFY;  
      `ps2pdf -dDEVICEWIDTHPOINTS=$weight -dDEVICEHEIGHTPOINTS=$height $infile $outfile`;  
      return $outfile;  
 }  

El Protein Data Bank en El País

Buenos días,
antes de seguir con las estructuras de datos, que dejamos para otro día, hoy me gustaría destacar que recientemente ha salido en el diario El País, uno de los más leídos en español, un artículo donde se resumían los últimos avances en las técnicas de resolución de estructuras moleculares, que poco a poco parece que se van encaminando al estudio de moléculas individuales. El artículo original, muy interesante,  es Destellos brillantes y ultracortos iluminarán la nueva biología estructural.
La verdad me sorprendió gratamente encontrarme en la prensa generalista con información reciente y relevante de acerca del Protein Data Bank (PDB), el recurso sobre el que se construye la Bioinformática Estructural, pero más me sorprendí al descubrir que en realidad el PDB había sido objeto de al menos otros 3 artículos en el mismo diario en años recientes, todos ellos firmados por Cele Abad Zapatero:
2004) La revolución de los rayos X
2007) ¿Morirá de éxito la biología estructural?
2010) Medio siglo de las primeras estructuras de proteínas

Si no habéis usado nunca el PDB, el enlace principal es http://www.rcsb.org,  y una entrada típica, como la 1le8, contiene imágenes como ésta:

Por supuesto podéis acceder a él desde vuestros programas Perl, por ejemplo con el módulo WWW::PDB. En el curso de Algoritmos en Bioinformática Estructural hay varios ejemplos de usos de archivos de coordenadas en formato PDB.

Un saludo y que vaya bien la semana,
Bruno

Rendimiento de estructuras de datos en Perl

Buenas,
recientemente, tratando de resolver las dudas de un principiante en Perl, he podido volver a constatar la dificultad que suponen para los novatos los vectores (@arrays) y las tablas asociativas (%hashes). Sin embargo, estas estructuras de datos son muy empleadas por usuarios experimentados, puesto que permiten guardar y modelar datos complejos de manera relativamente sencilla. No obstante, después de tanto tiempo usando estas estructuras, creo que nunca me he puesto a medir el coste de un vector y de un hash para guardar los mismos datos en un programita Perl. Este coste es especialmente relevante para manejar volúmenes de datos grandes, como averiguaremos con ayuda de los módulos Benchmark y Devel::Size.

Si además invocamos el módulo GraphViz::Data::Grapher, podremos dibujar ambos tipos de estructuras de datos, lo cual puede ayudar a los más principiantes:

 #!/usr/bin/perl -w  
   
 use strict;  
 use Benchmark;    
 use Devel::Size;   
 use Scalar::Util;  
 #use GraphViz::Data::Grapher;   
   
   
 my $MUESTRAS = 100;   
   
 # 1) mide tiempo de creacion de estructura de datos  
 print "# construyendo estructuras de datos:\n";  
 print "> hash :\n";  
 timethis( $MUESTRAS, "crea_estructura_datos('hash')" );  
 print "> array:\n";  
 timethis( $MUESTRAS, "crea_estructura_datos('array')" );  
   
 # 2) mide memoria que necesita cada estructura  
 print "\n# memoria RAM usada (KB):\n";  
 my $ref_hash = crea_estructura_datos('hash');  
 printf("\n> hash: %1.1f\n",  
    Devel::Size::total_size($ref_hash)/1024);  
   
 my $ref_array = crea_estructura_datos('array');  
 printf("> array: %1.1f\n",  
    Devel::Size::total_size($ref_array)/1024);  
      
 # 3) mide tiempo de consulta de estructura de datos  
 print "\n# consultando estructuras de datos:\n";  
 print "> hash:\n";  
 timethis( $MUESTRAS, sub{ consulta_estructura_datos($ref_hash) } );  
 print "> array:\n";  
 timethis( $MUESTRAS, sub{ consulta_estructura_datos($ref_array) } );     
   
 ############################################  
   
 sub crea_estructura_datos  
 {  
    my ($hash_o_array) = @_;   
   
     my $referencia;  
      
    if($hash_o_array eq 'hash')  
    {  
        foreach my $n (1..100_000)   
        {  
          # las llaves ocupan menos como cadenas de caracteres  
          # http://codenode.com/perl-memory-usage  
          $referencia->{"$n"} = $n * 10;         #7165.6KB en CPU 64bits  
         #$referencia->{sprintf("%01.3f",$n)} = $n * 10; #7556.2KB "  
          #$referencia->{$n/0.3} = $n * 10;        #7914.3KB "  
        }  
         
       #descomenta para generar grafo de estructura como el del blog  
       # OJO: mejor que sea una estructura no muy grande  
       #my $grafo = GraphViz::Data::Grapher->new($referencia);  
       #print $grafo->as_png("hash.png");  
     }  
    else  
    {  
       foreach my $n (1..100_000)   
        {  
         push(@{$referencia}, $n * 10);  # 3367.9KB           
          #$referencia->[$n-1] = $n * 10; # lo mismo  
        }  
         
       #my $grafo = GraphViz::Data::Grapher->new($referencia);  
       #print $grafo->as_png("array.png");  
    }  
   
    return $referencia;  
 }   
      
 sub consulta_estructura_datos  
 {  
    my ($referencia) = @_;   
      
    my $index;  
      
    if(Scalar::Util::reftype($referencia) eq "HASH")  
    {  
        foreach my $n (1..100_000)   
        {  
          $index = int(rand(100_000));  
          $referencia->{$index} += 1;   
        }  
     }  
    else  
    {  
       foreach my $n (1..100_000)   
        {  
         $index = int(rand(100_000));  
          $referencia->[$index] += 1;   
        }  
    }  
   
    return $referencia;  
 }   
   

Mediante este código, ejecutado en mi máquina de 64 bits, podemos ver que un hash sencillo ocupa un poco más del doble que un array para guardar una lista de enteros, y que su tamaño varía según el tipo de llave que usemos. Además, tardamos el triple de tiempo en llenar el hash que el array equivalente. Finalmente, este pequeño experimento muestra que el acceso y modificación de datos dentro un hash es al menos dos veces más lento que en un array:

# construyendo estructuras de datos:
> hash :
timethis 100: 10 wallclock secs (10.45 usr +  0.01 sys = 10.46 CPU) @  9.56/s (n=100)
> array:
timethis 100:  3 wallclock secs ( 2.17 usr +  0.00 sys =  2.17 CPU) @ 46.08/s (n=100)

# memoria RAM usada (KB):

> hash: 7165.6
> array: 3367.9

# consultando estructuras de datos:
> hash:
timethis 100:  6 wallclock secs ( 6.08 usr +  0.01 sys =  6.09 CPU) @ 16.42/s (n=100)
> array:
timethis 100:  3 wallclock secs ( 2.61 usr +  0.00 sys =  2.61 CPU) @ 38.31/s (n=100)

Un saludo,
Bruno