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

[Sendu Bala]

Marc Logghe wrote:
> Hi all,
> 
> I am a little bit confused about the above mentioned seq_inds() method.
> At first, I had the impression that the method returns an array of
> positions in the hsp (hit or query) sequence.

Yes...


> 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...


> 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.


(Actually, does 432 make sense? Shouldn't it be 431 twice?)