27 de septiembre de 2010

Ordenamiento de resultados de BLAST

Un problema con el que me he encontrado recientemente es el de ordenar resultados de BLAST contenidos en diferentes ficheros. Para definir mejor el problema supongamos que tenemos 4 genomas A,B,C y D y queremos comparar, todas contra todas, sus secuencias de proteínas por medio del programa blastp. Una manera de hacerlo sería hacer 4x4 comparaciones por parejas (AA,AB,AC,AD,BA,BB,BC,BD,CA,CB,CC,CD,DA,DB,DC,DD) teniendo en cuenta que la dirección en BLAST normalmente importa.

Una vez completados estos 16 trabajos de BLAST encontraremos que cada uno de los archivos de salida están internamente ordenados en términos de E-value. Sin embargo, dada la secuencia 1 del genoma A, si queremos averiguar en qué genomas se encuentra su secuencia más similar (best hit), deberemos mirar en todos los archivos (7 de 16) donde aparece el genoma A. Otra alternativa, la que da pie a este artículo, es fusionar, intercalar y ordenar (merge-sort) esos 16 archivos internamente ordenados para crear uno sólo, que facilite posteriores consultas.
Sirvan de muestra los archivos AA.blast:
16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8
y AB.blast:
16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103

que contienen, en formato tabular, los alineamientos significativos de la secuencia 16 dentro del genoma A (contra las secuencias 16,18,78) y dentro del genoma B (582,637), con valores esperados de 0.0, 2e-24,1e-19,5e-23,3e-21 respectivamente.
El problema que vamos a resolver es como combinar estos dos archivos (o en general N archivos) en uno sólo hasta obtener una lista ordenada por E-value:
16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8

Sin más preámbulos, el siguiente código Perl implementa dos posibles soluciones:

 package BlastSort;   
   
 use strict;  
 use Symbol qw(gensym);  
 use sort 'stable';   
   
 require Exporter;  
 use vars qw(@ISA @EXPORT);  
 @ISA = 'Exporter';  
 @EXPORT = qw(merge_BLAST_files merge_BLAST_files_GNUsort);  
   
 sub merge_BLAST_files   
 {  
    # Adapted from File-Sort-1.01(http://search.cpan.org/perldoc?File::Sort)  
    # Assumes infiles are BLAST output files in tabular format with sequences  
    # identified by natural numbers, such as 11,12,1439,1440 in the sample:  
    #11   1439   78.24   625   136   0   4   628   6   630   0.0    993  
    #12   1440   80.88   272   52   0   1   272   1   272   7e-125    446  
    # Order of infiles is RELEVANT as merging is stable, so sequences from  
    # the first files will be given sorting priority  
      
    my ($outfile,@infiles) = @_;  
      
    my (%fh,%order,@fhorder,$filein,$id,$first,$n_of_fhs,$curr,$line);  
   
    # for Schwartzian transform (ST) see  
    # http://www.hidemail.de/blog/perl_tutor.shtml#sort_orcish_schwartzian  
    sub blastsort { $a->[2] <=> $b->[2] || $a->[3] <=> $b->[3] }  
    sub blastmap   
    {   
       my @tmp = split(/\s+/,$fh{$_});   
       [$_,$fh{$_},$tmp[0],$tmp[10]]   
       # returns anonymous array[4]: filehandle, line,query_id,E-value  
    }  
   
    ## open all input BLAST files and keep filehandles, ordered 0..N  
    $n_of_fhs = 0;  
   foreach $filein (@infiles)   
    {  
      $id = gensym(); # get valid id for filehandle  
         
       open($id,$filein) ||   
          die "# merge_BLAST_files : cannot read $filein: $!";  
   
       $order{$n_of_fhs} = $id;  
       $n_of_fhs++;  
    }  
    @fhorder = (0 .. $n_of_fhs-1);  
      
     
    ## open outfile and ensure IO buffer is used  
    $| = 0;   
   unlink($outfile) if(-s $outfile);  
    open(OUT,">$outfile") ||   
       die "# merge_BLAST_files : cannot create $outfile: $!";  
   
   ## get first BLAST line from all filehandles    
   %fh = map {  
     my $fh = $order{$_};  
     ($_ => scalar <$fh>);  
   } @fhorder;  
   
   ## start merging BLAST lines   
    while(scalar(@fhorder)>1)   
    {        
       ($first) = (map {$_->[0]} sort blastsort map &blastmap, @fhorder); #ST  
         
       print OUT $fh{$first};  
   
     $curr = $order{$first};  
     $line = scalar <$curr>;  
     if(defined($line)) # update current filehandle  
       {   
          $fh{$first} = $line;   
       }  
       else # exhausted filehandle  
       {   
          @fhorder = grep { $_ ne $first } @fhorder;  
       }  
   }  
      
    ## take care of last filehandle left and close file  
    print OUT $fh{$fhorder[0]};  
    $curr = $order{$fhorder[0]};  
    while(<$curr>){ print OUT }  
   close(OUT);    
      
    ## close all input files  
    foreach $id (0 .. $n_of_fhs-1){ close($order{$id}); }  
      
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
   
 sub merge_BLAST_files_GNUsort  
 {  
    my ($outfile,$tmpdirpath,$maxbufferMb,@infiles) = @_;  
   
    # local sort -k11g fails with: 0.00,0.006 and 1e-23 (coreutils v6.10)  
    # probably as LC_ALL is not set to 'POSIX'  
    # http://www.gnu.org/software/coreutils/faq/coreutils-faq.html   
   # #Sort-does-not-sort-in-normal-order_0021  
    $ENV{'LC_ALL'} = 'POSIX';  
     
   unlink($outfile) if(-s $outfile);  
   
   my $sort_command = "sort --temporary-directory=$tmpdirpath " .  
     "--buffer-size=$maxbufferMb -s -k1g -k11g -m " .  
     join(' ',@infiles)." > $outfile ";  
   
   system("$sort_command");    
   
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
 __END__  

El módulo contiene dos subrutinas,  merge_BLAST_files y merge_BLAST_files_GNUsort; mientras la primera muestra explícitamente cómo se desarrolla el ordenamiento externo en disco, manteniendo los N ficheros de entrada abiertos y sin guardar nada en memoria, la segunda subrutina es realmente una llamada a la utilidad sort del shell (GNU coreutils 6.10 en mi sistema), donde sí usamos la memoria para acelerar el ordenamiento, 500Mb en el siguiente ejemplo:

 #!/usr/bin/perl  
 use strict;  
 use Benchmark;  
 use BlastSort;  
   
 my @infiles = ( 'A-A.blast','A-B.blast' );    
   
 my ($outfile,$gnu_outfile,$start_time,$end_time,$sortOK) =   
    ('out.merge-sort.txt','out.gnu-sort.txt');  
   
 print "# number of BLAST files to merge-sort = ".scalar(@infiles)."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files($outfile,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime (BlastSort): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files_GNUsort($gnu_outfile,'./',500,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime    (GNU): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  

Como se observa en la gráfica, aunque la versión Perl es autoexplicativa, es exponencialmente más lenta que GNUsort. Por tanto, en general convendrá usar la segunda opción y la primera sólo tendrá interés si queremos reservar la RAM, a costa de tiempo de ejecución. La única posible complicación del GNUsort es que depende de la variable de ambiente LC_ALL, a la que deberemos dar el valor 'POSIX' para tener resultados correctos. De no ser así no ordena bien los E-valores.


17 de septiembre de 2010

¿Cuál es el mínimo valor del E-value de Blast?

Esta extraña pregunta ha surgido hoy en mi laboratorio, y analizando varios ficheros de resultados de Blast he llegado a la conclusión que el menor E-value de Blast es 1e-180, por debajo de este valor Blast devuelve 0.0 como resultado.


¿What is the minimum Blast E-value?
This strange question has an easy solution: 1e-180, it is the minimum Blast E-value, lower values are scored by Blast as 0.0.

13 de septiembre de 2010

Jornadas Bioinformáticas JBI 2010 (X Edición), nuestro laboratorio estará allí...

Las Jornadas Bioinformáticas son la cita anual obligada para los bioinformáticos españoles. Este año se celebrará su décima edición del 27 al 29 de Octubre en Torremolinos (Málaga). La organización de las mismas corre a cargo de la Universidad de Málaga, el Instituto Nacional de Bioinformática y la Red Portuguesa de Bioinformática. Este año el tema central es "La bioinformática aplicada a la medicina personalizada", sobre el cual se discutirá la integración de los campos de la biología, medicina e informática para el desarrollo de terapias más específicas y efectivas. Sin embargo, éste no será el único tema a tratar, también se compartirán resultados y experiencias en otros campos:
- Análisis de datos en técnicas de alto rendimiento como la secuenciación de nueva generación.
- Bioinformática estructural
- Algoritmos de biología computacional y técnicas de computación de alto rendimiento
- Análisis de secuencias, filogenética y evolución
- Bases de datos, herramientas y tecnologías de biología computacional
- Bioinformática en transcriptómica y proteómica
- Biología sintética y de sistemas

IN ENGLISH:


The Xth Spanish Symposium on Bioinformatics (JBI2010) will take place in October 27-29, 2010 in Torremolinos-Málaga, Spain. Co-organised by the National Institute of Bioinformatics-Spain and the Portuguese Bioinformatics Network and hosted by the University of Malaga (Spain).


This year, the reference topic is “Bioinformatics for personalized medicine” for which the conference will provide the opportunity to discuss the state of the art for the integration of the fields of biology, medicine and informatics. We invite you to submit your work and share your experiences in the following topics of interest including, but not limited to:
- Analysis of high throughput data    (NGS)
- Structural Bioinformatics
 - Algorithms for computational biology and HPC
 - Sequence analysis, phylogenetics and evolution
 - Databases, Tools and technologies for computational biology
-  Bioinformatics in  Transcriptomics and Proteomics
- System and Synthetic Biology



Nuestras aportaciones

Nuestro laboratorio va a participar en las Jornadas Bioinformáticas con tres contribuciones que presentaré a continuación:


3D-footprint: a database for the structural analysis of protein–DNA complexes
3D-footprint is a living database, updated and curated on a weekly basis, which provides estimates of binding specificity for all protein–DNA complexes available at the Protein Data Bank. The web interface allows the user to: (i) browse DNA-binding proteins by keyword; (ii) find proteins that recognize a similar DNA motif and (iii) BLAST similar DNA-binding proteins, highlighting interface residues in the resulting alignments. Each complex in the database is dissected to draw interface graphs and footprint logos, and two complementary algorithms are employed to characterize binding specificity. Moreover, oligonucleotide sequences extracted from literature abstracts are reported in order to show the range of variant sites bound by each protein and other related proteins. Benchmark experiments, including comparisons with expert-curated databases RegulonDB and TRANSFAC, support the quality of structure-based estimates of specificity. The relevant content of the database is available for download as flat files and it is also possible to use the 3D-footprint pipeline to analyze protein coordinates input by the user. 3D-footprint is available at http://floresta.eead.csic.es/3dfootprint with demo buttons and a comprehensive tutorial that illustrates the main uses of this resource.


The relation between amino-acid substitutions in the interface of transcription factors and their recognized DNA motifs

Transcription Factors (TFs) play a key role in gene regulation by binding to DNA target sequences. While there is a vast literature describing computational methods to define patterns and match DNA regulatory motifs within genomic sequences, the prediction of DNA binding motifs (DBMs) that might be recognized by a particular TF is a relatively unexplored field. Numerous DNA-binding proteins are annotated as TFs in databases; however, for many of these orphan TFs the corresponding DBMs remain uncharacterized. Standard annotation practice transfer DBMs of well known TFs to those orphan protein sequences which can be confidently aligned to them, usually by means of local alignment tools such as BLAST, but these predictions are known to be error-prone. With the aim of improving these predictions, we test whether the knowledge of protein-DNA interface architectures and existing TF-DNA binding experimental data can be used to generate family-wise interface substitution matrices (ISUMs). An experiment with 85 Drosophila melanogaster homeobox proteins demonstrate that ISUMs: i) capture information about the correlation between the substitution of a TF interface residue and the conservation of the DBM; ii) are valuable to evaluate TFs alignments and iii) are better classifiers than generic amino-acid substitution matrices and that BLAST E-value when deciding whether two aligned homeobox proteins bind to the same DNA motif.


101DNA: a set of tools for Protein-DNA interface analysis

Analysis of protein-DNA interfaces has shown a great structural dependency. Despite the observation that related proteins tend to use the same pattern of amino acid and base contacting positions, no simple recognition code has been found. While protein contacts with the sugar-phosphate backbone of DNA provide stability and yield very little specificity information, contacts between amino acid side-chains and DNA bases (direct readout) apparently define specificity, in addition to some constrains defined by DNA sequence-dependent features, namely indirect readout.
Recent approaches have proposed bipartite graphs as an structural way of analysing interfaces from a protein-DNA-centric viewpoint. With this perspective in mind, we have developed a set of tools for the dissection and comparison of protein-DNA interfaces. Taking a protein-DNA complex file in PDB format as input, the software generates a 2D matrix that represents a bipartite graph of residue contacts  obtained after applying a simple distance threshold that captures all non-covalent interactions. The generated 2D matrices allow a fast and simple visual inspection of the interface and have been successfully produced for the current non-redundant set of protein-DNA complexes in the 3D-footprint database.
As a second utility to compare 2 interfaces, the 101DNA software includes an aligment tool where a dynamic programming matrix is created with the Local Affine Gap algorithm and traced back as a finite state automata. The scores between pairs of  interface amino acid residues are calculated as a function of the observed contacts with DNA nitrogen bases. This tool produces local interface alignments which are independent of the underlying protein sequence, but that faithfully represent the binding architecture. Preliminary tests show that these local alignments successfully identify binding interfaces that share striking similarity despite belonging to different protein superfamilies, and these observations support this graph-theory approach.

10 de septiembre de 2010

Cómo saber si un valor es numérico mediante una simple expresión regular

Hoy no estoy muy inspirado para escribir en el blog, así que he repasado mis librerías de subrutinas de perl para ver si encontraba algo curioso para publicar y me he encontrado con esto...
¿Quién no ha tenido nunca el problema de comprobar si una expresión es numérica o no? 
Perl no posee una función que lo haga automáticamente (por lo menos que yo conozca), sin embargo con una simple línea de código podemos salir de dudas:

 # Return TRUE if an expression is numerical  
 sub is_numeric {  
      my ($exp) = @_;  
      if ($exp =~ /^(-?[\d.\-]*e[\d.\-\+]+|-?[\d.\-]\^[\d.\-]+|-?[\d.\-]+)$/){  
      # Corregida siguiendo las indicaciones de Joaquín Ferrero:
      if ($exp =~ /^-?(?:[\d.-]+*e[\d.+-]+|\d[\d.-]*\^[\d.-]+|[\d.-]+)$/){
           return 1;  
      } else {  
           return 0;  
      }  
 }  

1 de septiembre de 2010

Adaptando scripts antiguos para Blast+

Hola,
los que uséis programas de la familia BLAST habitualmente habréis notado que las últimas versiones instalables, que se pueden descargar por FTP, son de la rama nueva de desarrollo Blast+, que tiene algunas nuevas funcionalidades muy interesantes, pero que cambia los nombres de los ejecutables y la forma de pasar argumentos que sabíamos usar.


Sin embargo, los desarrolladores del NCBI ya habían previsto esta transición y acompañan a los nuevos binarios un script Perl, que se llama legacy_blast.pl, que puede ayudarnos a reconvertir código que invoca versiones antiguas de BLAST. Con este programa podemos por ejemplo traducir este comando de la versión 2.2.18, que tiene opciones de filtrado un tanto especiales:

$ blastall -F 'm S' -z 10000 -p blastp -i problema.faa -d bases_datos/seqs.faa

a esta llamada al binario de la versión 2.2.24+:

$ blastp -db bases_datos/seqs.faa -query problema.faa -dbsize 10000 -seg yes -soft_masking true

Por cierto, en mis pruebas veo que el nuevo binario puede aprovechar una base de secuencias formateada con el formatdb antiguo, que ahora se llama makeblastdb.

Un saludo,
Bruno