[Bioperl-l] subtract for Bio::RangeI.pm

[Stephen Montgomery]

Hi -

It is a subtract function for the Bio::RangeI class.  (To be added if
interested)

All the best!
Stephen Montgomery


//ADD TO BIO::RANGEI



=head2 subtract

  Title   : subtract
  Usage   : my @subtracted = $r1->subtract($r2)
  Function: Subtract range r2 from range r1
  Args    : arg #1 = a range to subtract from this one (mandatory)
            arg #2 = strand option ('strong', 'weak', 'ignore')
(optional)
  Returns : undef if they do not overlap or r2 contains this RangeI,
            or an arrayref of Range objects (this is an array since some
            instances where the subtract range is enclosed within this
range
            will result in the creation of two new disjoint ranges)

=cut

sub subtract() {
   my ($self, $range, $so) = @_;
    $self->throw("missing arg: you need to pass in another feature")
      unless $range;
    return unless $self->_testStrand($range, $so);
    
    if ($self eq "Bio::RangeI") {
	$self = "Bio::Range";
	$self->warn("calling static methods of an interface is
deprecated; use $self instead");
    }
    $range->throw("Input a Bio::RangeI object") unless
$range->isa('Bio::RangeI');
    
    if (!$self->overlaps($range)) {
        return undef;
    }
    
    ##Subtracts everything
    if ($range->contains($self)) {
        return undef;   
    }
    
    my ($start, $end, $strand) = $self->intersection($range, $so);
    ##Subtract intersection from $self range
    
    my @outranges = ();
    if ($self->start < $start) {
        push(@outranges, 
		 $self->new('-start'=>$self->start,
			    '-end'=>$start - 1,
			    '-strand'=>$self->strand,
			   ));   
    }
    if ($self->end > $end) {
        push(@outranges, 
		 $self->new('-start'=>$end + 1,
			    '-end'=>$self->end,
			    '-strand'=>$self->strand,
			   ));   
    }
    return \@outranges;
}



//UNIT TEST

#!/usr/bin/perl
use strict;
use Bio::SeqFeature::Generic;
use Data::Dumper;
use Test;

plan tests => 13;

my $feature1 =  new Bio::SeqFeature::Generic ( -start => 1, -end =>
1000, -strand => 1);
my $feature2 =  new Bio::SeqFeature::Generic ( -start => 100, -end =>
900, -strand => -1);

my $subtracted = $feature1->subtract($feature2);
ok(defined($subtracted));
ok(scalar(@$subtracted) == 2);
foreach my $range (@$subtracted) {
    ok($range->start == 1 || $range->start == 901);
    ok($range->end == 99 || $range->end == 1000);
}

my $subtracted = $feature2->subtract($feature1);
ok(!defined($subtracted));
my $subtracted = $feature1->subtract($feature2, 'weak');
ok(!defined($subtracted));
my $subtracted = $feature1->subtract($feature2, 'strong');
ok(!defined($subtracted));

my $feature3 =  new Bio::SeqFeature::Generic ( -start => 500, -end =>
1500, -strand => 1);
my $subtracted = $feature1->subtract($feature3);
ok(defined($subtracted));
ok(scalar(@$subtracted) == 1);
my $subtracted_i = @$subtracted[0];
ok($subtracted_i->start == 1);
ok($subtracted_i->end == 499);