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.


No hay comentarios:

Publicar un comentario