[Bioperl-l] Parsing "hmmsearch" results - help

[Alberto Davila]

Indeed Jason,

I am listing my full code below.... could anybody kindly provide a
working code as an example for me ?

Thanks, Alberto

*****

#!/usr/bin/perl -w
 
printf("Name of the input Hmmer file ?\n");
my $inblastfile = <STDIN>;
chomp ($inblastfile);
 
printf("Name of the ouput file ?\n");
my $outputfilename = <STDIN>;
chomp ($outputfilename);
$outputfilenamenew = "$outputfilename".".hmmer.parsed";
 
my $logfilename = "$outputfilename".".hmmer.log";
my $fastafilename = "$outputfilename".".hmmer.fastalike";
 
printf("Minimum Hit E-value ?\n");
my $inevalue =<STDIN>;
chomp ($inevalue);
 
printf ("Minimum Hit Length ?\n");
my $inlength = <STDIN>;
chomp ($inlength);
 
unless (open(PARSEDFILE, ">$outputfilenamenew")) {
                print "Could not open file $outputfilenamenew !\n";
                exit;
        }
 
unless (open(LOGFILE, ">$logfilename")) {
                print "Could not open file $logfilename !\n";
                exit;
        }
 
unless (open(FASTAFILE, ">$fastafilename")) {
                print "Could not open file $fastafilename !\n";
                exit;
        }
 
 
my $albresult = 0;
my $albhit = 0;
my $albhsp = 0;
my $albsignificance = 0;
my $alblength = 0;
my $albnohit = 0;
my $alblasthit = 0;
my $alblastresult = 0;
my $albpercent = 0;
 
print PARSEDFILE  "Query @ Q-description @ Q-length @ Hit-Description @
H-Accession-Number @ H-Bit-Score @ E-value @ Identical @ Conserved @
Query_Start @ Query_End @ Target_Start @ Target_End\n";
 
    use strict;
    use Bio::SearchIO;
    use Bio::Search::Result::HMMERResult;
    
    my $in = new Bio::SearchIO(-format => 'hmmer',
                               -file   => $inblastfile);
 
    while( my $result = $in->next_result ) {
      $albresult ++;
 
 
                while( my $hit = $result->next_hit ) {
                $albhit++;
 
                                                                                
 
                        while( my $hsp = $hit->next_hsp ) {
                        $albhsp ++;
 
                                if( $hsp->length('total') >= $inlength )
{
           
                                $alblength++;
                   
                                if ( $hit->significance <= $inevalue ) {
           
                                $albsignificance++;
 
                        print PARSEDFILE $result->query_name,"@",
                                         $result->query_description,"@",
                                         $result->query_length, "@",
                                         $hit->description, "@",
                                         $hit->accession, "@",
                                         $hit->bits, "@",
                                         $hit->significance, "@",
                                         $hsp->num_identical, "@",
                                         $hsp->num_conserved,"@",
                                         $hsp->start('query'),"@",
                                         $hsp->end('query'),"@",
                                         $hsp->start('hit'),"@",
                                         $hsp->end('hit'),"@\n";
 
                        print FASTAFILE  "> ", $hit->description,"\n",
                                               $hsp->hit_string,"\n";
                                          
                                                                   }
                                                            }
                                                         }
 
                                                    }
      
 
         
    print LOGFILE "\nQuery < ", $result->query_name, " > : ", $albhit, "
hit(s)";
    if ($albhit == 0) {
               $albnohit++;
                      }
     
                                                                                
    $albhit = 0;
                                          }
       
    #$albpercent = $albnohit*100 / $albresult;
    print LOGFILE "\n\n", $albresult, " query sequence(s)\n";
    print LOGFILE "\n", $albhsp, " HSP sequence(s) \n";
    print LOGFILE "\n", $alblength, " hit(s) presenting the minimum
requested length\n";
    print LOGFILE "\n", $albsignificance, " hit(s) presenting the
minimum requested E-value\n";
    print LOGFILE "\n", $albnohit, " query sequence(s) presenting < NO
HITS > (",$albnohit*100/$albresult," %)\n\n";
     
    print "\n\n", $albresult, " query sequence(s)\n";
    print "\n", $albhsp, " HSP sequence(s) \n";
    print "\n", $alblength, " hit(s) presenting the minimum requested
length\n";    print "\n", $albsignificance, " hit(s) presenting the
minimum requested E-value\n";
    print "\n", $albnohit, " query sequence(s) presenting < NO HITS > ";
    printf "(%.2f %%)\n\n", $albnohit*100/$albresult;
    close (PARSEDFILE);
    close (LOGFILE);
    close (FASTAFILE);
    exit;





On Tue, 2004-04-27 at 05:34, Jason Stajich wrote:
> did you provide 'hmmer' as the format type?
> 
> 
> On Mon, 26 Apr 2004, Alberto Davila wrote:
> 
> > Hi all,
> >
> > I am trying to write a script to parse the "hmmsearch" results but not
> > sure I am using the appropriate code (eg "$hit->bits"). The code I am
> > listing below work ok parsing blast results but not "hmmsearch"
> > results... any tips ?
> >
> > Thanks, Alberto
> >
> >
> >     use strict;
> >     use Bio::SearchIO;
> >
> >
> >
> >              		while( my $hsp = $hit->next_hsp ) {
> > 	        	$albhsp ++;
> >
> >                    		if( $hsp->length('total') >= $inlength ) {
> >
> > 	             		$alblength++;
> >
> > 		  		if ( $hit->significance <= $inevalue ) {
> >
> > 	                	$albsignificance++;
> >
> >                      	print PARSEDFILE $result->query_name,"@",
> >                                          $result->query_description,"@",
> > 	     		       		 $result->query_length, "@",
> > 			       		 $hit->description, "@",
> > 			       		 $hit->accession, "@",
> >                                          $hit->bits, "@",
> >  			       		 $hit->significance, "@",
> > 			       		 $hsp->num_identical, "@",
> >                                          $hsp->num_conserved,"@",
> >                                          $hsp->start('query'),"@",
> >                                          $hsp->end('query'),"@",
> >                                          $hsp->start('hit'),"@",
> >                                          $hsp->end('hit'),"@\n";
> >
> >                         print FASTAFILE  "> ", $hit->description,"\n",
> > 			                       $hsp->hit_string,"\n";
> >
> >                                                                    }
> >                                                             }
> >                                                          }
> >
> >
> >
> > _______________________________________________
> > Bioperl-l mailing list
> > Bioperl-l at portal.open-bio.org
> > http://portal.open-bio.org/mailman/listinfo/bioperl-l
> >
> 
> --
> Jason Stajich
> Duke University
> jason at cgt.mc.duke.edu