Bioperl: Flipped out on Bio::Tools::Blast
Steve A. Chervitz
sac@alberich.Stanford.EDU
Tue, 25 Aug 1998 12:32:29 -0700 (PDT)
Lincoln,
This was an important oversight which has already been corrected in
version 0.06. There is now an HSP::strand() method for retrieval of strand
information. You can call the strand() method on HSP objects as in:
print $hsp->strand('query');
print $hsp->strand('sbjct');
These will return 'Minus' or 'Plus'.
Or, to get strand information for both query and sbjct (in that order):
@strands = $hsp->strand();
As you noted, the _set_seq() method will flip start <-> end if
start > end. This was done to facilitate HSP processing. Commenting this
out will have consequences only for TBLAST[NX] reports and only for methods
that rely on HSP::range() or HSP::matches(). The main effect is that the HSP
tiling code in Sbjct.pm could break for these reports, leading to
potentially incorrect values reported by methods of Sbjct.pm that rely on
HSP::_tile_hsp() such as length_aln(), gaps(), frac_identical(), etc. which
collect data across all HSPs. Tables produced by Blast::table_tiled() would
therefore be affected as well. (Since there is now a strand() method,
you shouldn't need to comment out any code :).
Thanks for reporting the bug in Sbjct.pm regarding sequence length
handling. I'll include it in the next release.
Steve Chervitz
sac@genome.stanford.edu
http://genome-www.stanford.edu/perlOOP/bioperl/blast/
On Tue, 25 Aug 1998, Lincoln Stein wrote:
> When I retrieve the list of HSPs for a hit, how do I tell whether
> I'm dealing with a plus or a minus strand HSP? The range() operator
> always sorts the endpoints of the subject and query sequences from
> lower to higher!
> ...
>
> Here's the offending code in Steve's HSP.pm. As far as I can tell, it
> silently flips the sequence endpoints without setting any flag that
> this has happened.
>
> # Check for flipped sequence endpoints (from reverse complement)
> if($self->{$seqType.'Start'} > $self->{$seqType.'Stop'}) {
> my $tmp = $self->{$seqType.'Start'};
> $self->{$seqType.'Start'} = $self->{$seqType.'Stop'};
> $self->{$seqType.'Stop'} = $tmp;
> }
>
> If I comment out this section, what will break?
>
> Lincoln
> --
> ========================================================================
> Lincoln D. Stein Cold Spring Harbor Laboratory
> lstein@cshl.org Cold Spring Harbor, NY
> ========================================================================
> =========== Bioperl Project Mailing List Message Footer =======
> Project URL: http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/
> For info about how to (un)subscribe, where messages are archived, etc:
> http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
> ====================================================================
>
=========== Bioperl Project Mailing List Message Footer =======
Project URL: http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/
For info about how to (un)subscribe, where messages are archived, etc:
http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
====================================================================