[Bioperl-l] Blast features added to wrong strand???

[michael watson (IAH-C)]

No, see something definitely weird is going on, but I am happy to accept
it may be my misuse.

To the bl2seq hits I have added the code "-report_type => 'blastn'" and
I get the same results.

The *really* weird thing is that after I have created these features, if
I write them to an EMBL file, they are annotated as being on the -1
strand e.g. complement(1..230) etc etc.  However, when I pass those very
same features to $panel->add_track, they are drawn on the + strand.

If I iterate through them and "print $feat->strand", they all say -1.
Write them out as EMBL, they say "complement(1..23) etc.  Draw them, and
they point --------------------> that way.

If I write them out as EMBL, then read them back in using Bio::SeqIO,
then pass them to $panel->add_track, they point in the right direction.
So something is getting set wrong - could it be "frame"?

????

-----Original Message-----
From: Jason Stajich [mailto:jason.stajich at duke.edu] 
Sent: 14 July 2005 13:41
To: michael watson (IAH-C)
Cc: bioperl-l at portal.open-bio.org
Subject: Re: [Bioperl-l] Blast features added to wrong strand???


Why use bl2seq when you can use fasta....

 From the Bio::SearchIO::blast documentation

=head2 bl2seq parsing

Since I cannot differentiate between BLASTX and TBLASTN since bl2seq
doesn't report the algorithm used - I assume it is BLASTX by default -
you can supply the program type with -report_type in the SearchIO
constructor i.e.

   my $parser = new Bio::SearchIO(-format => 'blast',
                                  -file => 'bl2seq.tblastn.report',
                                  -report_type => 'tblastn');

This only really affects where the frame and strand information are put
- they will always be on the $hsp-E<gt>query instead of on the
$hsp-E<gt>hit part of the feature pair for blastx and tblastn bl2seq
produced reports.  Hope that's clear...


On Jul 14, 2005, at 6:53 AM, michael watson ((IAH-C)) wrote:

> Hi
>
> I'm using bioperl-1.4.
>
> I have a genomic region.  I am first using bl2seq and blastn to
> align it
> with some custom sequences, then blastall and blastx to blast it  
> agains
> uniprot.
>
> I'm using SearchIO and pretty standard code to add the blast hits as 
> features to the sequence e.g.:
>
> my $feature = Bio::SeqFeature::Generic->new(-primary_tag  => 'CDS',
>                               -score
> => $hit->raw_score,
>                               -display_name
> => $hit->name,
>                                -tag
> => {
>                                                              
> locus_tag =>
> $name,
>
> note      => $note,
>                                                               }
>                                 );
>
> # @hsps is a filtered list of HSPs obtained from $hit->next_hsp 
> foreach $hsp (@hsps) {
>     $feature->add_sub_SeqFeature($hsp,'EXPAND');
> }
>
> $genome->add_SeqFeature($feature); # $genome is a Bio::Seq feature
>
> Now, the bl2seq hits all have the strand reported as "Plus / Minus"
> and
> the blastx hits all have the strand reported as -1 i.e. there is a  
> gene
> on the other strand of my sequence.
>
> HOWEVER, using the above code for both the bl2seq results and the
> blastx
> results, ONLY the blastx results get annotated on the reverse strand -
> the bl2seq results, which report the strand as "Plus / Minus", get
> annotated on the forward strand and hence point the wrong way when I
> draw them :-(
>
> So my question is, what am I doing wrong in the above code (which is 
> pretty much ripped off from the bioperl HOWTOs) that makes the bl2seq 
> "Plus / Minus" hits get annotated on the plus strand on my sequence??
>
> Many thanks
> Mick
>
> _______________________________________________
> 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/