Alineamiento global con el algoritmo Wavefront (WFA)

Hola, en esta entrada vuelvo a una de las piedras angulares de la biología computacional, el alineamiento de secuencias. Este problema tiene múltiples caras, algunas ya discutidas en este blog, pero desde el algoritmo de Needleman-Wunsch su formulación fundamental es el alineamiento global de dos secuencias q y t de longitudes n y m, desde el principio hasta el final. 

La última vuelta de tuerca la acaban de publicar Santiago Marco Sola y colaboradores en Bioinformatics, donde describen su algoritmo de alineamiento llamado wavefront (WFA). En el artículo los autores muestran como WFA y su variante heurística WFA-Adapt son más eficientes tanto en tiempo de cálculo como en consumo de memoria que las alternativas del estado del arte, incluyendo los algoritmos de la librería KSW2, empleada por  minimap2 (visto en este blog).

Recomiendo la lectura del artículo porque es muy didáctico y está escrito en un estilo sencillo, dada la naturaleza del problema. A continuación resumo aquí las ideas principales. 

El problema que resuelve WFA es un alineamiento global entre dos secuencias usando el modelo affine-gap como función coste. Esto significa que para cada par de posiciones alineadas el coste del alineamiento aumenta en 0 puntos en caso de ser idénticas (a), en x puntos en caso de ser diferentes (x=6 en BWA-MEM) o con un coste lineal para las inserciones que se calcula como g(n) = o + e n
donde o es el coste de abrir un indel y e el de extenderlo n bases. WFA es un algoritmo exacto que calcula el alineamiento óptimo como aquel que termina en la celda (n,m) de la matriz de programación dinámica (Figura 1, izquierda) con un coste total más pequeño. Hasta aquí nada nuevo, es un algoritmo que busca diagonales que alinean las dos secuencias.

La novedad de WFA es que define furthest-reaching points (fr), vectores  Fs,k que indican para una diagonal k el punto más lejano donde se alcanza un coste s (ver Figura 1 izquierda, vectores M0, M4 y M8 desde el origen, donde M=matches, de la misms manera que I=indel y D=deletion). En su algoritmo reformulan el alineamiento por programación dinámica calculando vectores fr para un coste s en base a los vectores fr calculados para costes menores, pero de manera que sólo una fracción de los fr se llegan a calcular. El algoritmo descrito en la Figura 2 recibe su nombre porque para cada coste s se define el frente de onda WFs como el conjunto de todos los vectores fr con coste total s. El alineamiento optimo se corresponde a la secuencia de frentes de onda desde WF0 a WFs que alcanzan la coordenada (n, m) con el menor coste s. En el ejemplo de la Figura 1 solamente es necesario guardar en memoria 3 WF (0, 4 y 8) para calcular el alineamiento óptimo y luego reconstruirlo. A diferencia de otros algoritmos de alineamiento, WFA es más eficiente cuánto más se parecen las secuencias a alinear y sus operaciones son fácilmente paralelizables de manera portable con instrucciones SIMD.

 

Figura 1. Alineamiento global de dos secuencias en una matriz de programación dinámica (izq) y su representación en forma de vectores fr (furthest-reaching points, dcha). En la matriz las celdas (0,0) y (6,6) marcan respectivamente el principio y el final del alineamiento. Tomada de Bioinformatics. 

Figura 2. Pseudocódigo para el algoritmo WFA, tomado de Bioinformatics.

Para terminar esta entrada, el código de WFA está escrito en C y se compila fácilmente con gcc si lo clonas desde el repositorio https://github.com/smarco/WFA . Allí encontrarás ejemplos sencillos de cómo llamar a las funciones de alineamiento, un saludo,

Bruno

 

Bioinformática en C (y perl6 async)

Hola, y feliz año!
tras la última entrada, y antes de que escriba algo más completo sobre perl6,
parece que la gestión de tareas asíncronas tiene mucho potencial...

Pero vamos a lo que vamos. Hoy quería compartir un documento sobre cómo programar en C con los compiladores actuales, porque algunos todavía seguimos usando sintaxis anticuada y esto es relevante en nuestro campo, dado que hay muchas aplicaciones bioinformáticas escritas en C, normalmente por cuestiones de eficiencia.

El documento en cuestión es: https://matt.sh/howto-c

y ofrece un montón de trucos y recomendaciones para modernizar y simplificar la escritura de código en C, con ejemplos sencillos como éste para vectores dinámicos:

 uintmax_t arrayLength = strtoumax(argv[1], NULL, 10);
    void *array[arrayLength];

 /* Hace innecesario liberar el puntero al final.

    Ojo: debe ser un tamaño razonable para el sistema donde se ejecute. */

Ha sido traducido al español aquí y comentado/criticado aquí y en menéame,
un saludo,
Bruno

kseq::klib to parse FASTQ/FASTA files

Buenas,
éste tiene pinta de ser el último artículo del año, no? Desde hace tiempo venimos usando en el laboratorio la librería kseq.h , integrante de [ https://github.com/attractivechaos/klib ], para las desagradecidas tareas de manipular grandes archivos en formato FASTQ. Esta librería es un único archivo .h en lenguaje C que de manera muy eficiente permite leer archivos FASTQ, y también FASTA, incluso comprimidos con GZIP, con zlib. La eficiencia es fundamental para estas aplicaciones, porque modificar archivos .fastq de varios GB puede consumir mucho tiempo.

Fuente: http://drive5.com/usearch/manual/fastq_files.html

Como resultado de esta experiencia hace unos meses colgamos de la web del labo un programita, que llamamos split_pairs. Este software facilita trabajos tediosos como pueden ser modificar las cabeceras de una archivo FASTQ con expresiones regulares para que un programa posterior no se queje, y de paso emparejar o interlazar los reads según convenga.  Está escrito en C++ y tiene asociado un script Perl.

De lo que quería hablar hoy es de que es posible utilizar kseq.h directamente desde nuestro lenguaje preferido. Para ellos es necesario emplear herramientas que permitan incluir, y compilar, código C mezclado con el de nuestro lenguaje interpretado. En el caso de Python estaríamos hablando por ejemplo de cython , aquí os muestro cómo hacerlo en Perl con el módulo Inline::CPP, dotando a kseq de una interfaz orientada a objetos para evitar pasar tipos de datos con tipos no estándar:

 package klib::kseq;  
 use Inline CPP => config => namespace => 'klib';    
 use Inline CPP => config => libs => '-lz';  
 use Inline CPP => <<'END_CPP';  
 /*  
 The Klib library [ https://github.com/attractivechaos/klib ], is a standalone and   
 lightweight C library distributed under MIT/X11 license.   
 Klib strives for efficiency and a small memory footprint.  
 This package contains C++ wrappers for components of Klib to make them easy to use from Perl scripts.  
 */  
 using namespace std;  
 #include <iostream>  
 #include <string>  
 #include <zlib.h>  
 #include "kseq.h"   
 /*   
 kseq is generic stream buffer and a FASTA/FASTQ format parser from Klib.  
 Wrapper written by Bruno Contreras Moreira.  
 */  
 KSEQ_INIT(gzFile, gzread)  
 class kseqreader
 {  
 public:  
    kseqreader()  
    {  
       fp = NULL;  
       fromstdin=false;  
       openfile=false;  
       filename = "";  
    }  
    ~kseqreader()  
    {  
       close_stream();  
    }  
    void rewind_stream()  
    {  
       if(openfile == true)  
       {  
          if(fromstdin == false) kseq_rewind(seq);  
          else  
             cerr<<"# reader: cannot rewind STDIN stream"<<endl;  
       }  
    }           
    void close_stream()  
    {  
       if(openfile == true)  
       {  
          kseq_destroy(seq);  
          if(fromstdin == false) gzclose(fp);  
          fp = NULL;  
          fromstdin=false;  
          openfile=false;  
          filename = "";  
       }     
    }  
    // accepts both FASTA and FASTQ files, which can be GZIP-compressed  
    // returns 1 if file was successfully opened, 0 otherwise  
    int open_stream(char *infilename)  
    {  
       filename = infilename;  
       if(filename.empty())  
       {   
          cerr<<"# reader: need a valid file name"<<endl;  
          filename = "";  
          return 0;  
       }     
       else  
       {  
          FILE *ifp = NULL;  
          if(filename == "stdin")  
          {  
             ifp = stdin;  
             fromstdin = true;  
          }  
          else ifp = fopen(filename.c_str(),"r");  
          if(ifp == NULL)     
          {  
             cerr<<"# reader: cannot open input file: "<<filename<<endl;  
             filename = "";  
             return 0;  
          }  
          else // open file and initialize stream  
          {  
             fp = gzdopen(fileno(ifp),"r");   
             seq = kseq_init(fp);   
             openfile=true;  
             return 1;  
          }  
       }     
    }  
     // returns length of current sequence, which can be 0,  
    // -1 in case of EOF and -2 in case of truncated quality string  
    int next_seq()  
    {   
       return kseq_read(seq);   
    }  
     const char *get_sequence()  
    {  
       return seq->seq.s;  
    }  
    // returns FASTA/FASTQ header up to first found space  
    const char *get_name()  
    {  
       return seq->name.s;  
    }  
    const char *get_quality()  
    {  
       if(seq->qual.l) return seq->qual.s;  
       else return "";   
    }  
    const char *get_comment()  
    {  
       if(seq->comment.l) return seq->comment.s;  
       else return "";   
    }  
    const char *get_full_header()  
    {  
       _header = seq->name.s;  
       if(seq->comment.l)  
       {   
          _header += " ";  
          _header + seq->comment.s;  
       }     
      return _header.c_str();  
    }  
 private:  
    gzFile fp;  
    kseq_t *seq;  
    bool fromstdin;  
    bool openfile;  
    string filename;  
    string _header;  
 };  
 END_CPP  
 1;  

Ahora estamos en disposición de probar a leer un archivo con este módulo:

 use strict;  
 use klib;  
 my ($length,$header,$sequence,$quality);  
 my $infile = "sample.fq.gz";  
 # my $infile = "stdin";  
 my $kseq = new klib::kseq::kseqreader();  
 if($kseq->open_stream($infile))  
 {  
   while( ($length = $kseq->next_seq()) != -1 )  
   {  
     ($header,$sequence,$quality) =  
       ($kseq->get_full_header(),$kseq->get_sequence(),$kseq->get_quality());   
   }  
 }  
 $kseq->close_stream()  

Qué ganancia en tiempo de procesamiento de archivos FASTQ se obtiene? Una prueba realizada con el módulo Benchmark , con diez repeticiones y un archivo .fastq.gz con 250K secuencias sugiere que kseq::klib es el doble de rápido:

            s/iter  perlio   kseq
perlio    1.32      --     -50%
kseq    0.660  100%     --

Hasta el próximo año, un saludo,
Bruno

C++ STL en paralelo

Hola,
como contábamos en una entrada anterior, últimamente hemos estado liados trabajando con archivos FASTQ de varios GBs, con decenas de millones de lecturas o reads. Para algunas de las tareas que hemos tenido que realizar, como ordenar o renombrar, nos hemos dado cuenta que lenguajes interpretados como Perl o Python son varias veces más lentos que caballos de carreras compilados como C/C++, como también comenta el autor de readfq, por ejemplo.
Por eso he estado refrescando la Standard Template Library (STL) de C++, que para alguien acostumbrado a programar en lenguajes de alto nivel es esencial, como recurso para crear estructuras de datos flexibles y eficientes en programas escritos en C/C++. Como dice Scott Meyers en su libro Effective STL, probablemente los contenedores más populares de la STL son vectores y strings. Sólo por estas dos estructuras dinámicas, a las que yo añado también las tablas asociativas (maps), vale la pena aprender a usar la STL. Sin embargo, hay mucho más que esto, y por eso escribo esta entrada, porque la STL implementada en el compilador GNU gcc (libstdc++) incluye una biblioteca de algoritmos en paralelo de propósito general, que podemos usar fácilmente en nuestros programas para sacarle el jugo a los procesadores multicore y resolver de manera más eficiente múltiples problemas. Entre otros, la versión actual de libstdc++ paraleliza los siguientes algoritmos de la STL, todos ellos útiles para tareas habituales en programación:

• std::count
• std::find
• std::search
• std::replace
• std::max_element
• std::merge
• std::min_element
• std::nth_element
• std::set_union
• std::set_intersection
• std::sort

Esquema de sort parelelo y merge posterior,
tomado de http://javaero.org/tag/parallel-merge-sort

El siguiente ejemplo de sort en paralelo, que depende de la librería OpenMP, se compila con $ g++ -O3 -fopenmp -o testP testP.cc para generar un ejecutable que empleará 4 cores para ordenar un millón de palabras de 25 caracteres de longitud:

 
#include <stdlib.h>  
#include <omp.h>  
#include <vector>  
#include <parallel/algorithm>  
using namespace std;  

// g++ -O3 -fopenmp -o testP testP.cc  
// http://gcc.gnu.org/onlinedocs/libstdc++/manual/parallel_mode.html  
#define MAXTHREADS 4;  
 
string randomStrGen(int length) 
{
   //http://stackoverflow.com/questions/2146792/how-do-you-generate-random-strings-in-c    
   static string charset = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:|.@";  
   string result;  
   result.resize(length);  
   for (int i = 0; i < length; i++) result[i] = charset[rand() % charset.length()];  
   return result;  
}  

int main()  
{  
   unsigned threads = MAXTHREADS;  
   omp_set_dynamic(false);  
   omp_set_num_threads(threads);  
   srand(time(NULL));  
   std::vector<string> data(1000000);  
   for(int i=0;i<data.size();i++) data[i] = randomStrGen(25);  
   __gnu_parallel::sort(data.begin(), data.end()); //std::less<string>() , std::smaller<string>() );  
   printf("%s\n%s\n%s\n%s\n",data[0].c_str(),data[1].c_str(),data[2].c_str(),data[data.size()-1].c_str());  
   return 0;  
}  

Referencias más avanzadas:
[1] http://algo2.iti.kit.edu/singler/mcstl/parallelmode_se.pdf

[2] http://ls11-www.cs.uni-dortmund.de/people/gutweng/AD08/VO11_parallel_mode_overview.pdf

Hasta luego,
Bruno

comparando secuencias de igual longitud

Hola,
hoy quería compartir recetas en Perl para la comparación eficiente de secuencias de igual longitud. Las dos funciones que vamos a ver devuelven como resultado listas con las posiciones (contando desde 0) en las que difieren dos secuencias.

Para qué sirve esto? En general, para calcular distancias de edición (edit distances); en bioinformática es una manera, por ejemplo, de comparar oligos de igual longitud o secuencias extraídas de un alineamiento múltilple, incluyendo gaps. El siguiente código fuente incluye una subrutina en Perl (con el operador binario ^ y la función pos) y otra en C embebido para hacer la misma operación por dos caminos:

 use strict;  
 use warnings;  
 my $seq1 = 'ACTGGA';  
 my $seq2 = 'AGTG-A';  
   
 my @Pdiffs = find_diffs($seq1,$seq2);  
 my @Cdiffs = find_diffs_C($seq1,$seq2);  
   
 printf("# version Perl\ndiferencias = %d\n%s\n",scalar(@Pdiffs),join(',',@Pdiffs));  
 printf("# version C\ndiferencias = %d\n%s\n",scalar(@Cdiffs),join(',',@Cdiffs));  
   
 sub find_diffs  
 {  
    my ($seqA,$seqB) = @_;  
    my $XOR = $seqA ^ $seqB; 
    my @diff;   
    while($XOR =~ /([^\0])/g)  
    {   
       push(@diff,pos($XOR)-1);  
    }  
    return @diff;  
 }  
   
 use Inline C => <<'END_C';  
 void find_diffs_C(char* x, char* y)   
 {                        
   int i;   
    Inline_Stack_Vars;                               
   Inline_Stack_Reset;                                                                
   for(i=0; x[i] && y[i]; ++i) {                          
    if(x[i] != y[i]) {                              
     Inline_Stack_Push(sv_2mortal(newSViv(i)));                 
    }                                       
   }   
    Inline_Stack_Done;                                                                  
 }     
 END_C  

Podéis ver en stackoverflow otras versiones y su comparación en base a su tiempo de cálculo,
un saludo,
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

Código en C dentro de un programa Perl

Hola, siguiendo el hilo de los comentarios a la entrada anterior de Álvaro, hoy pongo un ejemplo de como incluir una o más subrutinas escritas en C dentro de un texto en Perl, recurriendo al módulo Inline:

 #!/usr/bin/perl -w 
 # programa 7 
 use Inline C; 
 use strict; 
 my $BASE = 0.987654321; 
 my $EXP = 1000; 
 print "# $BASE^$EXP = \n"; 
 print "# ".potencia_perl($BASE,$EXP)." (perl)\n"; 
 print "# ".potencia_C($BASE,$EXP)." (C)\n"; 
 sub potencia_perl 
 { 
    my ($base,$exp) = @_; 
    my $result = $base; 
    for(my $p=1; $p<$exp; $p++) { $result *= $base; } 
    return $result; 
 } 
 __END__  
 ### alternativa en C 
 __C__ 
 #include <stdio.h> 
 float potencia_C( double base, double exp )  
 { 
    double result = base; 
    int p; 
    for(p=1; p<exp; p++) 
    {  
       result *= base; 
    }    
    return result; 
 } 

La salida obtenida pone en evidencia que diferentes lenguajes suelen tener diferentes precisiones a la hora de representar números y operar con ellos:

$ perl prog7.pl 
# 0.987654321^1000 = 
#  4.02687472213071e-06 (perl)
#  4.02687464884366e-06 (C)