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