11 de enero de 2013

Indexando archivos FASTA

Hola y feliz año nuevo con retraso!
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";  

No hay comentarios:

Publicar un comentario