Go Back   SEQanswers > Bioinformatics > Bioinformatics

Similar Threads
Thread Thread Starter Forum Replies Last Post
Hello, and thanks in advance Champi Introductions 0 05-11-2013 08:02 PM
Bioperl installation chandrugcg Bioinformatics 2 01-28-2013 03:45 PM
Bioperl Script to convert fasta to fq Lizex Bioinformatics 0 01-26-2012 11:21 PM
Reminder: X-Gen Congress Advance Registration Deadline and Exclusive Savings LifeScienceMarketing Events / Conferences 0 01-23-2012 05:02 AM
Bioperl jsun529 Bioinformatics 7 09-01-2009 08:27 AM

Thread Tools
Old 06-05-2013, 09:26 PM   #1
Lyn Hsiong
Location: China

Join Date: Sep 2011
Posts: 14
Default Errors in my bioperl script. Thanks in advance!

Dear friends:

I have a blast report, and want to extract following information for each result(query): the Query_name, hit_number, name and description of the hit(HSP) with the highest identity.

my bioperl script is:
#!/usr/bin/perl -w
use strict;
use warnings;
use Bio::SearchIO;

if (@ARGV != 2){die "Usage: $0 <BLAST-report-file> <output-file>\n"};

my ($infile,$outfile) = @ARGV;
print "Parsing the BLAST result ...";
my $blast = new Bio::SearchIO(
-format => 'blast',
-file   => $infile);
open (OUT,">$outfile") or die "Cannot open $outfile: $!";

print OUT "query\tNumber of hits\tGreatest identity %\tAccession (identity %)\tDescription (identity %)\n";

while (my $result = $blast->next_result){
   print OUT $result->query_name . "\t";
   print OUT $result->num_hits. "\t";
        if (my $hit = &sort_hit($result)){
                         if (my $hsp=$hit->hsp){
                                        print OUT $hsp->percent_identity. "\t"; 
                                        print OUT $hit->accession. "\t";
                                        print OUT $hit->description. "\n";
 close OUT;
 print " DONE!!!\n";

 sub sort_hit{
                  my $result = shift;
                  my @hits;
                  while (my $hit = $result->next_hit) {
                        push @hits, $hit;
                  my @hit_sort = sort { $b-> hsp ->percent_identity * 100 <=> $a-> hsp ->percent_identity * 100 } @hits;
                  my $flag=0;
                  my $temp_hit;
                  while (my $hit=$result->next_hit){
                  return $temp_hit
error information:
-------------EXCEPTION: Bio::Root::Exception -------------
MSG: Can't get HSPs: data not collected.
STACK: Error::throw
STACK: Bio::Root::Root::throw /home/xiongtl/Bioperl/dag/lib/perl5/Bio/Root/
STACK: Bio::Search::Hit::GenericHit::hsp /home/xiongtl/Bioperl/dag/lib/perl5/Bio/Search/Hit/
STACK: main::sort_hit
Parsing the BLAST result ...
I don't know what does it mean. Please give me some clues, all your words are welcome!

BTW, my bioperl version is 1.6.1, I use blastx 2.2.27+.
My blast report(-outfmt 0) have many queries and each query have many hits(against nr database), so each hit maybe also have several HSPs.

Last edited by Lyn Hsiong; 06-06-2013 at 06:50 PM. Reason: extra information
Lyn Hsiong is offline   Reply With Quote
Old 06-06-2013, 04:16 PM   #2
Senior Member
Location: The University of Melbourne, AUSTRALIA

Join Date: Apr 2008
Posts: 275

Do you mean $hit->next_hsp instead of $hit->hsp ?
Torst is offline   Reply With Quote
Old 06-06-2013, 05:45 PM   #3
Lyn Hsiong
Location: China

Join Date: Sep 2011
Posts: 14

Originally Posted by Torst View Post
Do you mean $hit->next_hsp instead of $hit->hsp ?
Thanks for your reply!
I tried the $hit->next_hsp before, but running the script still informed me the same error. However, today I have just reruned the script, both $hit->next_hsp and $hit->hsp are ok! They both outputed the first hit for each result. I don't know why. So I tried different sort rules as follows:

# my @hit_sort = sort { $b-> hsp ->hsp_length <=> $a-> hsp ->hsp_length } @hits;
# my @hit_sort = sort { $b-> hsp ->percent_identity * 100 <=> $a-> hsp ->percent_identity * 100 } @hits;
# my @hit_sort = sort { $b-> significance <=> $a-> significance } @hits;
# my @hit_sort = sort { $a-> accession cmp $b-> accession } @hits;
but different sort functions still generate the same output (still the first hit)! It seems that the sub-function does not work at all!
oh, it really troubles me...
Lyn Hsiong is offline   Reply With Quote

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

All times are GMT -8. The time now is 09:28 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO