Skip to content
Snippets Groups Projects
Commit f6753cbf authored by Arne Stabenau's avatar Arne Stabenau
Browse files

handle features that only have seqnames (and no slices!)

parent 93cb056d
No related branches found
No related tags found
No related merge requests found
......@@ -792,9 +792,11 @@ sub seq {
Example : if( $this_feature->overlaps( $otherfeature ) ...
Description: Test if two features overlap. Strandedness is ignored. Both features have
to have slices. If they are not on the same coordinate system they need
as well database connection.
Returntype : boolean 0,1
Exceptions : on no slices, wrong argument type, more if not on same coord_system
as well database connection. Calculates how many bases overlap.
There is some backward compatibility for unsliced features built in here.
Its worth checking once a while wether this can go.
Returntype : int
Exceptions : wrong argument type, if not on same coord_system
Caller : general
=cut
......@@ -803,41 +805,63 @@ sub overlaps {
my $self = shift;
my $feature = shift;
if( ! $feature->isa( "Bio::EnsEMBL::Feature" )) {
throw( "Need a Bio::EnsEMBL::Feature as argument" );
my ( $sstart, $astart, $send, $aend, $sname, $aname );
my ( $scoordsystem, $acoordsystem, $overlap );
if( defined $self->{'slice'} ) {
$sname = $self->seq_region_name();
$sstart = $self->seq_region_start();
$send = $self->seq_region_end();
$scoordsystem = $self->slice->coord_system();
} else {
$sname = $self->seqname();
$sstart = $self->start();
$send = $self->end();
}
if( ! defined $self->slice() || ! defined $feature->slice() ) {
throw( "Features need to contain slices to support overlaps" );
if( defined $feature->{'slice'} ) {
$aname = $feature->seq_region_name();
$astart = $feature->seq_region_start();
$aend = $feature->seq_region_end();
$acoordsystem = $feature->slice()->coord_system();
} else {
$aname = $feature->seqname();
$astart = $feature->start();
$aend = $feature->end();
}
if( $self->slice()->coord_system()->equals
( $feature->slice()->coord_system())) {
if( $self->seq_region_name() eq $feature->seq_region_name() ) {
if( $feature->seq_region_start() <= $self->seq_region_end() &&
$self->seq_region_start() <= $feature->seq_region_end() ) {
return 1;
} else {
return 0;
}
} else {
return 0;
}
# simple case
if( $sname eq $aname ) {
$overlap = (($send>$aend) ? $aend : $send ) -
(($sstart>$astart) ? $sstart : $astart ) +1;
if( $overlap < 0 ) { $overlap = 0;}
return $overlap;
}
# not same coord system
# have to project
my $projection = $self->project( $feature->slice()->coord_system()->name(),
$feature->slice()->coord_system()->version());
for my $proj ( @$projection ) {
my $slice=$proj->[2];
if( $slice->seq_region_name() eq $feature->seq_region_name() &&
$slice->start() <= $feature->seq_region_end() &&
$feature->seq_region_start() <= $slice->end() ) {
return 1;
}
# full feature case
# if both have coord systems
if( defined $scoordsystem && defined $acoordsystem ) {
if( $scoordsystem->equals( $acoordsystem )) {
# they are really not on the same sequence
return 0;
} else {
my $projection = $self->project( $feature->slice()->coord_system()->name(),
$feature->slice()->coord_system()->version());
my $total_overlap = 0;
for my $proj ( @$projection ) {
my $slice=$proj->[2];
if( $slice->seq_region_name() eq $feature->seq_region_name()) {
$overlap = (($slice->end()>$aend) ? $aend : $slice->end() ) -
(($slice->start()>$astart) ? $slice->start() : $astart ) +1;
if( $overlap < 0 ) { $overlap = 0;}
$total_overlap += $overlap;
}
}
return $total_overlap;
}
} else {
return 0;
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment