[Bioperl-l] Bio::Search::HSP::GenericHSP::seq_inds

Marc Logghe Marc.Logghe at ablynx.com
Thu Mar 27 15:09:56 UTC 2008


Hi Sendu, Chris

> > At least that is what one would expect looking at the example usage
in
> > the HOWTOs
(http://www.bioperl.org/wiki/HOWTO:SearchIO#Using_the_methods
> > second code block).
> >
> > Am I correct in believing you can only do this if your hsp query
stretch
> > starts at position 1 of the query sequence?
> 
> No...

<HOWTO>
# put all the conserved matches in query strand into an array
my @str_array = split "",$hsp->query_string;
foreach ( $hsp->seq_inds('query','conserved') ){
  push @conserved,$str_array[$_ - 1];
}
</HOWTO>

$hsp->query_string will return
'IAVEEETKTTKKNKKQ-QQQANKNKNKNKKK--TTIAPEAAIDANIAAEVHTQVL'

In my example using the 'gap' class (instead of 'conserved'), @str_array
will contain 417, 431 and 432. The off-by-one indices do not exist in
that array.
Therefore, I still think the howto shows a special case where the hsp
query sequence starts at 1 (compared to 402 in my particular example). 


> 
> 
> > Query  402  IAVEEETKTTKKNKKQ-QQQANKNKNKNKKK--TTIAPEAAIDANIAAEVHTQVL
453
> >             +A+ E   TT K +KQ ++  NK  NK  KK  T+  P+AA+ + I AE  +Q L
> > Sbjct  139  VALREALSTTFKTQKQLKENLNKQPNKQWKKSQTSTNPKAALKSKIVAEFRSQAL
193
> >
> > Now, if you call seq_inds(query => 'gap') on that particular hsp
object,
> > you get these positions: 417, 431, 432. Obviously, there is no gap
in
> > the original query sequence at these positions.
> > How do you have to read these numbers ? Remark also that for
instance
> > 417 is the res just in front of the gap.
> 
> Its purpose is to let you know the position in query or subject
> coordinates where something interesting happened in the alignment. So
> seq_inds(query => 'gap') is telling you all the places that a gap
starts
> in the alignment in terms of the query coordinates. Hence 417 etc.

So, this means you have to interpret that as a gap is coming after 417 ?

> 
> 
> (Actually, does 432 make sense? Shouldn't it be 431 twice?)
Don't know, depends on how you have to 'read' this.
Thanks for looking into this.
Regards,
Marc





More information about the Bioperl-l mailing list