[Bioperl-l] Alignment->slice() issue?
Kevin Brown
Kevin.M.Brown at asu.edu
Wed Jun 17 15:41:18 UTC 2009
Warning: This is very ugly code and makes a few assumptions, such as the
alignment objects are stored in order of their start position. I made
this assumption as that is how I put them into the object to begin with.
=head1 C<slice>
Function to slice up an alignment sequence based on start and end
parameters
and returns a new alignment object.
slice($alignment, $start, $end)
=cut
sub slice
{
my ($alignment, $start, $end, $new_align) = @_;
$$new_align = new Bio::SimpleAlign;
print $$alignment->no_sequences() . "\n";
$$new_align->add_seq(
new Bio::LocatableSeq(
-seq =>
substr(
$$alignment->get_seq_by_pos(1)->seq(),
$start - 1, $end
- $start + 1
),
-id =>
$$alignment->get_seq_by_pos(1)->display_id(),
-start =>
max($$alignment->get_seq_by_pos(1)->start - $start + 1, 1),
-end => min(
$$alignment->get_seq_by_pos(1)->end - $start + 1,
$end - $start
+ 1
),
-alphabet => 'dna',
-strand =>
$$alignment->get_seq_by_pos(1)->strand()
)
);
# implement a binary search to determine a decent offset into
the alignment
my $probe;
if ($$alignment->no_sequences() <= 2) {
$probe = $$alignment->no_sequences();
}
else {
my ($L, $R) = (1, $$alignment->no_sequences());
while (($R - $L) > 1)
{
$probe = floor(($R + $L) / 2);
# gotta watch this. Had the check backwards and so was
never going
# in the right direction for the search. If I reverse
these two
# variables, then I have to either reverse the
conditions or change
# the > to a <.
if ($$alignment->get_seq_by_pos($probe)->start() >
$start)
{
$R = $probe;
}
else
{
$L = $probe;
}
}
}
# now go through the results that are after that point
for (my $i = $probe ; $i <= $$alignment->no_sequences() ; $i++)
{
my $seq = $$alignment->get_seq_by_pos($i);
last if ($seq->start() > $end);
# Only concern ourselves with primers that land inside
the desired region
# other primers will show up in the image maps for each
gene.
if ($seq->start() >= $start && $seq->end() <= $end)
{
# values for the substr pullout of a given
sequence
my $offset = max($start - $seq->start(), 0);
my $length =
min($end, $seq->end()) - max($start,
$seq->start()) + 1;
$$new_align->add_seq(
new Bio::LocatableSeq(
-seq => $seq->seq(),
-id =>
$seq->display_id(),
-start =>
max($seq->start - $start + 1, 1),
-end => min($seq->end -
$start + 1, $end - $start + 1),
-alphabet => 'dna',
-strand =>
$seq->strand()
)
);
}
}
return 1;
}
> -----Original Message-----
> From: bioperl-l-bounces at lists.open-bio.org
> [mailto:bioperl-l-bounces at lists.open-bio.org] On Behalf Of
> Malcolm Cook
> Sent: Tuesday, June 16, 2009 1:07 AM
> To: bioperl-l at lists.open-bio.org
> Subject: [Bioperl-l] Alignment->slice() issue?
>
> Kevin,
>
> I'm getting struck by this old issue you once coded around.
>
> http://bioperl.org/pipermail/bioperl-l/2007-January/024665.html
>
> Any chance you could share your implementation with fellow
> traveller...
>
> ??
>
> Thanks,
>
> Malcolm Cook
> Stowers insitute for Medical research
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>
More information about the Bioperl-l
mailing list