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

[Kevin Brown]

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
>