[Bioperl-l] SearchIO.pm: recalculate position between hit and qry

[Jason Stajich]

Coordinate mapping is done in Bio::Coordinate objects.  New code which is
only in CVS and not in 1.4 release will take an alignment and allow you to
ask what base in the alignment corresponds to position in the query or
hit.

This *should* work but is untested.

# convert the HSP to a Bio::SimpleAlign
my $hsp_aln = $hsp->get_aln;
# Build a mapper object
my $mapper = Bio::Align::Utils->from_seq_to_alignmentpos($hsp_aln);
print "SNP at base 501 is base ", $mapper->map(501), " in the alignment\n";
If you want to connect that to the Hit sequence
my $base_mapped = $mapper->map(501);

print "501 in the query is base ",
$hsp_aln->location_from_column($base_mapped)->start,"\n";

(you get back a location in the event the remapped base lies on a gap).

Someday this sort of stuff needs to go into an Align HOWTO....
-jason

On Mon, 7 Jun 2004, Jan Aerts wrote:

> Hi all,
>
> Every once in a while I BLAST a sequence against dbSNP to annotate SNPs
> on it. To find the exact position of the SNP on the query, I find myself
> writing code to check if the SNP position (which is known on the 'hit'
> sequence) is within the alignment and recalculating the position on the
> query, taking into account
> *start and end positions of both query and hit
> *strand of hit
> *any gaps on query or hit.
>
> I was not able to find code in SearchIO.pm to do this for me. Would it
> be a good idea to include it?
>
> Possible use:
>   print $hsp->qry_pos('501');
> where 501 is the position on the hit (e.g. the position of a SNP within
> its flanking sequence). And the other way around:
>   print $hsp->hit_pos('2965');
> where 2965 is the position on the query.
>
> Any suggestions?
> jan.
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l
>

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