Recientemente me he visto en la necesidad de acceder de maneara aleatoria a secuencias de un archivo FASTA de gran tamaño, en mi caso cercano a los 2 Gb,
con las secuencias formateadas en una sola línea.
Primero probé con las herramientas de BLAST+, en concreto:
$ makeblastdb -in archivo.fasta -parse_seqid
Con la idea de posteriormente consultar las secuencias con algo como (más ejemplos aquí):
$ blastdbcmd -query -db archivo.fasta
Sin embargo, la generación del índice se hizo eterna y de hecho la aborté, por problemas que desconozco y no he seguido mirando.
Posteriormente, tras buscar en Google encontré dos caminos en Perl:
1) con ayuda de scripts de Bioperl: bp_index.pl y bp_fetch.pl
2) con ayuda del módulo core Tie::File, como muestro en el ejemplo a continuación, mi solución preferida, un saludo:
#!/usr/bin/perl -w
use strict;
use Tie::File;
my $fasta_file = '/path/to/file.fasta';
my (%index_fasta,@fasta_array);
open(FASTA,$fasta_file) ||
die "# cannot open $fasta_file\n";
while(<FASTA>)
{
if(/^> any typical marker in header, such as gi (\S+)/)
{
$index_fasta{$1} = $.;
}
}
close(FASTA);
printf("# indexed %d sequences in file %s\n\n",
scalar(keys(%index_fasta)),$fasta_file);
tie(@fasta_array,'Tie::File',$fasta_file) ||
die "# cannot tie file $fasta_file\n";
print "ejemplo: secuencia con identificador 'id12':\n";
print "$fasta_array[$index_fasta{'id12'}-1]\n".
"$fasta_array[$index_fasta{'id12'}]\n";
