[Bioperl-l] Alignment->slice() issue?

[bill at genenformics.com]

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
>