SEQanswers

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 07:02 PM
Bioperl installation chandrugcg Bioinformatics 2 01-28-2013 02:45 PM
Bioperl Script to convert fasta to fq Lizex Bioinformatics 0 01-26-2012 10:21 PM
Reminder: X-Gen Congress Advance Registration Deadline and Exclusive Savings LifeScienceMarketing Events / Conferences 0 01-23-2012 04:02 AM
Bioperl jsun529 Bioinformatics 7 09-01-2009 07:27 AM

Reply
 
Thread Tools
Old 06-05-2013, 08:26 PM   #1
Lyn Hsiong
Member
 
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:
Code:
#!/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;
                  $result->rewind;
                  my $flag=0;
                  my $temp_hit;
                  while (my $hit=$result->next_hit){
                        if($flag==0){
                           $temp_hit=$hit;
                     $flag++;
                           }
                  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/Root.pm:368
STACK: Bio::Search::Hit::GenericHit::hsp /home/xiongtl/Bioperl/dag/lib/perl5/Bio/Search/Hit/GenericHit.pm:688
STACK: main::sort_hit xtl_new.pl:37
STACK: xtl_new.pl:20
-----------------------------------------------------------
Parsing the BLAST result ...
I don't know what does it mean. Please give me some clues, all your words are welcome!


NOTE:
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 05:50 PM. Reason: extra information
Lyn Hsiong is offline   Reply With Quote
Old 06-06-2013, 03:16 PM   #2
Torst
Senior Member
 
Location: The University of Melbourne, AUSTRALIA

Join Date: Apr 2008
Posts: 275
Default

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

Join Date: Sep 2011
Posts: 14
Default

Quote:
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:

Quote:
# 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
Reply

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 01:01 PM.


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