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.
Cualquiera que lleve cierto tiempo trabajando en Bioinformática seguro que tiene programitas y scripts que llevan tiempo languideciendo en directorios o incluso en copias de seguridad condenadas al olvido. Lo que propone el autor es que muchos de esos programas se podrían publicar por ejemplo como material suplementario junto con los artículos para los que se escribieron. Yo añadiría que son todavía más visibles, y por tanto potencialmente más útiles para otros usuarios, si los colgamos de la web en blogs como éste.
