[Bioperl-l] Alignment->slice() issue?
bill at genenformics.com
bill at genenformics.com
Wed Jun 17 17:03:16 UTC 2009
Hopefully this is helpful.
http://www.ncbi.nlm.nih.gov/IEB/ToolBox/CPP_DOC/lxr/source/src/objects/seqalign/Dense_seg.cpp#L648
Bill at genenformics
> 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
>>
>
> _______________________________________________
> 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