my ($self, $object, $to_assembly) = @_;
throw("Need an assembly version to project to.") unless ($to_assembly);
throw("Need an object to project.") unless ($object and ref($object));
my ($slice, $object_type);
if ($object->isa('Bio::EnsEMBL::Feature')) {
$object_type = 'feature';
} elsif ($object->isa('Bio::EnsEMBL::Slice')) {
$object_type = 'slice';
} else {
throw("Need a Feature or Slice to project.");
}
# if the feature or slice is sourced from another db, we have to "transfer"
# it to the db that contains the assembly mapping. the transfer is very
# shallow but that should do for our purposes
if ($self->external_source) {
my $slice_adaptor = $self->adaptor->get_SliceAdaptor;
if ($object_type eq 'feature') {
# createa a new slice from the target db
my $f_slice = $object->slice;
my $target_slice = $slice_adaptor->fetch_by_name($f_slice->name);
# now change the feature so that it appears it's from the target db
$object->slice($target_slice);
} else {
# createa a new slice from the target db
$object = $slice_adaptor->fetch_by_name($object->name);
}
}
if ($object_type eq 'feature') {
$slice = $object->feature_Slice;
} else {
$slice = $object;
}
# warn if trying to project to assembly version the object already is on
if ($slice->coord_system->version eq $to_assembly) {
warning("Assembly version to project to ($to_assembly) is the same as your object's assembly (".$slice->coord_system->version.").");
}
# now project the slice
my $cs_name = $slice->coord_system_name;
my @segments = @{ $slice->project($cs_name, $to_assembly) };
# we need to reverse the projection segment list if the orignial
if ($slice->strand == -1) {
@segments = reverse(@segments);
}
# apply rules to projection results
#
# discard the projected feature/slice if
# 1. it doesn't project at all (no segments returned)
# 2. [unless MERGE_FRAGMENTS is set] the projection is fragmented (more
# than one segment)
# 3. [if CHECK_LENGTH is set] the projection doesn't have the same length
# as the original feature/slice
# 4. all segments are on same chromosome and strand
# keep track of the status of applied rules
my @status = ();
# test (1)
return undef unless (@segments);
#warn "DEBUG: passed test 1\n";
# test (2)
return undef if (!($self->merge_fragments) and scalar(@segments) > 1);
push @status, 'fragmented' if (scalar(@segments) > 1);
#warn "DEBUG: passed test 2\n";
# test (3)
my $first_slice = $segments[0]->to_Slice;
my $last_slice = $segments[-1]->to_Slice;
my $length_mismatch = (($last_slice->end - $first_slice->start + 1) !=
$object->length);
return undef if ($self->check_length and $length_mismatch);
push @status, 'length_mismatch' if ($length_mismatch);
#warn "DEBUG: passed test 3\n";
# test (4)
my %sr_names = ();
my %strands = ();
foreach my $seg (@segments) {
my $sl = $seg->to_Slice;
$sr_names{$sl->seq_region_name}++;
$strands{$sl->strand}++;
}
return undef if (scalar(keys %sr_names) > 1 or scalar(keys %strands) > 1);
#warn "DEBUG: passed test 4\n";
# remember rule status
$self->last_status(join('|', @status));
# everything looks fine, so adjust the coords of your feature/slice
my $new_slice = $first_slice;
$new_slice->{'end'} = $last_slice->end;
if ($object_type eq 'slice') {
return $new_slice;
} else {
$object->start($new_slice->start);
$object->end($new_slice->end);
$object->strand($new_slice->strand);
$object->slice($new_slice->seq_region_Slice);
# undef dbID and adaptor so you can store the feature in the target db
$object->dbID(undef);
$object->adaptor(undef);
return $object;
}
}