[Bioperl-l] Determining strand orientation from bl2seq report

[Jason Stajich]

Dave - a little bit of detective work.  This stems from the fact that
bl2seq reports don't include information about which blast program was
used (blastn, blastp, etc) and that determines whether or not we set
stranded information or set it to 0 (for proteins).  Surprised this was
not caught before, I've added tests and testfiles to make sure this
doesn't happen again.

I've committed changes to CVS allow you to specify the program used and
added tests for blastn bl2seq reports.  You can see how it works by
looking at t/BPbl2seq.t .

Basically when you initialize a BPbl2seq object you add in the flag

-report_type => 'blastn'

You'll need to get the changes from CVS and either drop that file in your
distro or do a full CVS checkout.

We may support bl2seq in SearchIO in the future if someone wants to try
and do this I'm more than happy to look over the code / patch to the
'blast' SearchIO parser.

-jason
On Fri, 24 Jan 2003, Dave Arenillas wrote:

> Hi Jason,
>
> Thanks for the info. I forgot to mention in my original post that I did try
> using $hsp->query->strand and $hsp->hit->strand and for some reason
> $hsp->hit->strand always returns 0. I also tried the alternate $hsp->hstrand
> to no avail. Is this a bug in bioperl? It is critical to my application that
> I am able to determine weather the blast HSPs are Plus / Plus or Plus /
> Minus. I guess my next step is to try to run blast in such a way as to be
> able to parse with SearchIO and see if that correctly returns strand info.
> This is a bit of a pain though.
>
> Dave
>
>
> On Friday 24 January 2003 08:27 am, Jason Stajich wrote:
> > Searchio can't parse bl2seq (yet), you have to use Bio::Tools::BPbl2seq,
> > so StandAloneBlast returns a Bio::Tools::BPbl2seq report.
> >
> > start is always less than in end in the bioperl world to get strand you
> > call $hsp->query->strand or $hsp->hit->strand.
> > The HSP object you get back is a FeaturePair object so you can call all
> > the methods that are available there, in addition the $hsp->hit and
> > $hsp->query methods return objects which are Bio::SeqFeatureI compliant
> > so you can call any method available in Bio::SeqFeatureI on them.
> >
> > perldoc doesn't do a good job with the fact that methods inherit from a
> > superclass so you don't see all those methods in the documentation, if you
> > use doc.bioperl.org and click through the inheritance that is available
> > in the ISA part of a module documentation you will be able to see what
> > are all the methods.
> >
> > -jason
> >
> > On Thu, 23 Jan 2003, Dave Arenillas wrote:
> > > I'm trying to figure out if there is any way to determine the strand
> > > orientation (Plus / Plus or Plus / Minus) of an HSP from the blast report
> > > created by running bl2seq.
> > >
> > > ex.
> > > my @params = ('program' => 'blastn',
> > >                 'outfile' => "$BlastFile");
> > > my $factory = Bio::Tools::Run::StandAloneBlast->new(@params);
> > > my $report = $factory->bl2seq($seq1, $seq2);
> > >
> > > while (my $hsp = $report->next_feature) {
> > > 	#
> > > 	# get orientation info from $hsp???
> > > 	#
> > > }
> > >
> > > Apparently Bio::Tools::BPlite::HSP objects have no methods to return this
> > > information. I (incorrectly) assumed that I could test
> > > $hsp->hit->start > $hsp->hit->end
> > > but start is always less than end even for Plus / Minus orientation.
> > >
> > > I then attempted to create a SearchIO object and use this to parse the
> > > blast report, for ex.
> > > my $searchio = new Bio::SearchIO( -format => 'blast',
> > >                         -file   => "$BlastFile" );
> > >
> > > while (my $result = $searchio->next_result) {
> > > 	#
> > > 	# Results in the error below
> > > 	#
> > >         while (my $hit = $result->next_hit) {
> > >                 while (my $hsp = $hit->next_hsp) {
> > > 		}
> > > 	}
> > > }
> > >
> > > However, this produced the following error:
> > >
> > > ------------ EXCEPTION  -------------
> > > MSG: no data for midline Lambda     K      H
> > > STACK Bio::SearchIO::blast::next_result
> > > /usr/lib/perl5/site_perl/5.8.0/Bio/SearchIO/blast.pm:566
> > >
> > > Is the report produced by bl2seq some kind of simplified blast report
> > > which SearchIO cannot parse? Is there any other way of determining HSP
> > > strand orientation?
> > >
> > > Thanks,
> > > Dave
> > >
> > >
> > > _______________________________________________
> > > Bioperl-l mailing list
> > > Bioperl-l at bioperl.org
> > > http://bioperl.org/mailman/listinfo/bioperl-l
>
>

--
Jason Stajich
Duke University
jason at cgt.mc.duke.edu