parsing BLAST output

Anonymous Monk has asked for the wisdom of the Perl Monks concerning the following question:

Hi. i am trying to write a program that will parse the output of a blastx file so that the output looks likr this;


Gene: blahblah
Query: FDVKHGAHJKHBHXYJTJGHCJHZGJHJJFVJKFG
[download]

mb, a BLAST output file looks like this;

    Query= 8
         (545 letters)

Database: /ebi/cgg/one/nrdb-filtered/nrdb.filtered
           905,591 sequences; 283,075,510 total letters

Searching..................................................done

                                                                   Sco
+re     E
Sequences producing significant alignments:                        (bi
+ts)  Value

/:swiss|P75505|YC71_MYCPN Hypothetical lipoprotein MPN271 precu...   1
+36  2e-31
>/:swiss|P75505|YC71_MYCPN Hypothetical lipoprotein MPN271 precursor
           (A65_orf251a).//:pironly|S73889|S73889 probable
           lipoprotein A65_orf251a - Mycoplasma pneumoniae (ATCC
           29342) (SGC3)//:gp|AE000055|1674262 conserved
           hypothetical protein, see: MPN643 [Mycoplasma
           pneumoniae]
           Length = 251

 Score =  136 bits (340), Expect = 2e-31
 Identities = 72/121 (59%), Positives = 88/121 (72%), Gaps = 1/121 (0%
+)
 Frame = -3

Query: 363 KRSFYQELPLLLWFLLLGTVLSACSKVESTKSVHPVKSFKSDLKSLLKETIKQDIDLSK
+T 184
           K+ F         FLL GT+LSAC+ +++            DL++L+KET  +DIDLSK
Sbjct: 3   KKKFLSRFSFSSLFLLCGTLLSACTGIQA------------DLRNLIKETTGKDIDLSK
+A 50
[download]

Here is my attempt below which just doesn't run (p.s i realise that this could be easily achieved using bioperl but i dont want to use that). Error messages especially don't like it when i mention $end_annotation and $start_annotation!! thanks.


#! /usr/local/bin/perl -w

use strict;

my $start_annotation;
my $end_annotation;
my %alignments = ();
my $filename = 'strept_blastx.output';

parse_blast (\$start_annotation, \$end_annotation, \%alignments, $file
+name);

print $start_annotation;

sub get_file_data
{
my ($filename) = @_;
my @filedata = ();

unless(open (GET_FILE_DATA, $filename))
{
print STDERR "cant open file \"$filename\"\n\n";
exit;
}   
 @filedata = <GET_FILE_DATA>;
close GET_FILE_DATA;
return @filedata;
}


foreach my $key (keys %alignments)
{
print "$key\nxxxxxxxxx\n", $alignments{key}, "\nxxxxxxxxxx\n";
}

print $end_annotation;

exit;

sub parse_blast
{
my ($start_annotation, $end_annotation, $alignments, $filename) = @_;

my $blast_output_file = '';
my $alignment_section = '';

$blast_output_file = join('', get_file_data ($filename));   

   ($$start_annotation, $alignment_section, $$end_annotation) = ($blas
+t_output_file
=~/(.*^ALIGNMENTS\n)(.*)(^ Database:.*)/ms);

%alignments = parse_blast_alignment ($alignment_section);
}

sub parse_blast_alignment
{
my($alignment_section) = @_;
my (%alignment_hash) =();

while ($alignment_section =~ /^>.*\n(^(?!>).*\n)+/gm)
{
my ($value) = $&;
my ($key) = (split(/\|/, $value)) [1];
$alignment_hash{$key} = $value;
}

return %alignment_hash;
}
[download]

Comment on parsing BLAST output Select or Download Code

Replies are listed 'Best First'.
Re: parsing BLAST output by Hofmator (Curate) on May 31, 2002 at 13:10 UTC
Well, let me just rewrite your program a bit ... but be cautious, this is untested code so I couldn't check your regexes, ... (and I probably introduced some small bugs of my own :) #! /usr/local/bin/perl -w use strict; my $filename = 'strept_blastx.output'; my ($start_annotation, $end_annotation, $alignments) = parse_blast ($f +ilename); print $start_annotation; print map "$_\nxxxxxxxxx\n$alignments->{$_}\nxxxxxxxxxx\n", keys %$ali +gnments; print $end_annotation; sub parse_blast { my ($filename) = @_; my $blast_output_file; my ($start_anno, $alignment_sec, $end_anno); open( my $data_file, '<', $filename ) or die "Couldn't open file $fi +lename: $!"; $blast_output_file = do {local $/; <$data_file> }; close $data_file; ($start_anno, $alignment_sec, $end_anno) = ($blast_output_file =~/(.^ALIGNMENTS\n)(.)(^ Database:.)/ms); my $align_hashref = parse_blast_alignment($alignment_sec); return ($start_anno, $end_anno, $align_hashref); } sub parse_blast_alignment { my ($alignment_section) = @_; my $alignment_hashref; while ($alignment_section =~ /^>.\n(^(?!>).*\n)+/gm) { my $value = $&; my ($key) = (split(/\\|/, $value)) [1]; $alignment_hashref->{$key} = $value; } return $alignment_hashref; } [download] So, this should run and then you can tell us if you still have some problems with your regexes and the actual extraction. Some more explanations for us non-biologists about the format would then be helpful as well. -- Hofmator	[reply] [d/l]
Re: parsing BLAST output by marvell (Pilgrim) on May 31, 2002 at 11:05 UTC
I'd start with http://bio.perl.org. -- Steve Marvell	[reply]
Re: parsing BLAST output by indapa (Monk) on May 31, 2002 at 18:01 UTC
check out BPlite.pm The documentation is here	[reply]
Re: parsing BLAST output by mrbbking (Hermit) on May 31, 2002 at 12:34 UTC
Where in the BLAST output do you find the Gene? I don't see '[Gg]ene' anywhere in the sample you posted (and I'm not sure what else to look for...) As for the Query: `Query: 363 KRSFYQELPLLLWFLLLGTVLSACSKVESTKSVHPVKSFKSDLKSLLKETIKQDIDLSK +T 184 K+ F FLL GT+LSAC+ +++ DL++L+KET +DIDLSK` [download] - where does it end? Is it both lines, or just the first? Is the three-digit number included?	[reply] [d/l]