Perl: get specific base from FASTA file.

[njh_TO]

I've written a piece of Perl code that retrieves a base from a FASTA file, given the contig label (e.g. chr1) and base position. However- I think it runs slowly, and am sure is not optimal: mainly because requests are usually made in an ordered manner, which is not reflected by my code.

Is anyone aware of a function that has already been written to do this well?

[Dameon]

samtools faidx /path/to/your/genome.fa chr1:base_position

[njh_TO]

Fast response- thanks

[Dameon]

No problem. I was just incorporating that little bit into one of my perl scripts yesterday. Glad it could help.

[njh_TO]

FYI:
using Benchmark:

Rate mine samtools
mine 111/s -- -92%
samtools 1449/s 1201% --

[krobison]

Are you using Bio:: DB ::Fasta? (Apologies for the spaces; otherwise the I get a smilie version: Bio:B::Fasta)

For your benchmark, are you calling samtools repeatedly or once? If you have multiple requests, better to batch them up and call samtools -- there's an overhead in process creation and probably in firing up samtools.

[njh_TO]

Thanks. No, I'm not- but I'll take a look.

Not sure how you would go about bundling these requests, is it possible to pass samtools an array of positions, returning a corresponding array of bases?

Benchmarking calls the following function X times.

Currently looks something like:

sub getFASTABaseSamtools
{
my $chrID = shift;
my $basepos = shift;
my $fastapath = shift;
my @result = `samtools faidx $fastapath $chrID:$basepos-$basepos`;
chomp($result[1]);
return uc($result[1]);
}