[Bioperl-l] Bio::SearchIO::hmmer

[Jason Stajich]

It really wasn't designed to parse A0 format.  I think we've slowly  
tried to plug the gaps.

Okay I've fixed it, give the latest code from CVS a whirl.

And I don't see how your code is going to work wrt printing out the  
hit_string if you don't have any alignments in the file.

-jason
On Jun 1, 2005, at 9:08 AM, Sean O'Keeffe wrote:

> Hi,
> I was wondering how Bio::SearchIO::hmmer parses hmmpfam/hmmsearch
> result files. I have a set of hmmpam result hits without alignments in
> a file which I generated using hmmpfam locally (-A0 option on the
> command line). Is this considered a valid result by
> Bio::SearchIO::hmmer? If so, what might be wrong with the code below
> (which gets as far as the first while loop and doesn't enter any
> other):
>
> use strict;
> use Bio::SearchIO;
>
> my $inhmmfile = "test-hmm.smart";
> my $outputfilename = "HMM-test.hmmer.parsed";
> my $fastafilename = "$outputfilename".".fasta";
> my $inevalue =1;
> my $inlength =20;
> my ($myresult,$myhit,$myhsp,$mysignificance,$mylength,
> $mynohit,$mylasthit,$mylastresult,$mypercent) = 0;
>
> unless (open(PARSEDFILE, ">$outputfilename")) {
>                print "Could not open file $outputfilename !\n";
>                exit;
>        }
>
> unless (open(FASTAFILE, ">$fastafilename")) {
>                print "Could not open file $fastafilename !\n";
>                exit;
>        }
>
> my $in = new Bio::SearchIO(-format => 'hmmer', -file => $inhmmfile);
>
> while(my $result = $in->next_result ) {
>        $myresult++;
>        while (my $hit = $result->next_hit ) {
>                $myhit++;
>                while (my $hsp = $hit->next_hsp ) {
>                        $myhsp++;
>                        if( $hsp->length('total') >= $inlength ) {
>                                $mylength++;
>                                if ( $hit->significance <=  
> $inevalue ) {
>                                        $mysignificance++;
>
>                                        print PARSEDFILE
> $result->query_name,"\t",
>                                        $result- 
> >query_description,"\t",
>                                        $result->query_length, "\t",
>                                        $hit->description, "\t",
>                                        $hit->accession, "\t",
>                                        $hit->bits, "\t",
>                                        $hit->significance, "\t",
>                                        $hsp->num_identical, "\t",
>                                        $hsp->num_conserved,"\t",
>                                        $hsp->start('query'),"\t",
>                                        $hsp->end('query'),"\t",
>                                        $hsp->start('hit'),"\t",
>                                        $hsp->end('hit'),"\n";
>
>                                        print FASTAFILE  "> ",
> $hit->description,"\n",
>                                        $hsp->hit_string,"\n";
>                                }
>                        }
>                }
>        }
>
>        if ($myhit == 0) {
>               $mynohit++;
>        }
>
>        $myhit = 0;
> }
>
> $mypercent = $mynohit*100 / $myresult;
>
> print "\n\n", $myresult, " query sequence(s)\n";
> print "\n", $myhsp, " HSP sequence(s) \n";
> print "\n", $mylength, " hit(s) presenting the minimum requested  
> length\n";
> print "\n", $mysignificance, " hit(s) presenting the minimum requested
> E-value\n";
> print "\n", $mynohit, " query sequence(s) presenting < NO HITS > ";
> close (PARSEDFILE);
> close (FASTAFILE);
> exit;
>
>
> Output is:
> Use of uninitialized value in numeric eq (==) at ../hmm-test-parse.pl
> line 58, <GEN1> line 40126.
>
>
> 1 query sequence(s)
>
> Use of uninitialized value in print at ../hmm-test-parse.pl line 68,
> <GEN1> line 40126.
>  HSP sequence(s)
>
> Use of uninitialized value in print at ../hmm-test-parse.pl line 69,
> <GEN1> line 40126.
>  hit(s) presenting the minimum requested length
>
> Use of uninitialized value in print at ../hmm-test-parse.pl line 70,
> <GEN1> line 40126.
>  hit(s) presenting the minimum requested E-value
>
> 1 query sequence(s) presenting < NO HITS > (100.00 %)
>
>
> Thanks very much,
> Sean.
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>

--
Jason Stajich
Duke University
http://www.duke.edu/~jes12/