Skip to content
Snippets Groups Projects
Commit 1f714d13 authored by Ian Longden's avatar Ian Longden
Browse files

_merge_ajoining_coords added to make sure that having more than one coord is...

_merge_ajoining_coords added to make sure that having more than one coord is wrong. Andy Yates wrote the code i just checked it and added it.
parent c25b660f
No related branches found
No related tags found
No related merge requests found
...@@ -1231,29 +1231,39 @@ sub peptide { ...@@ -1231,29 +1231,39 @@ sub peptide {
unless($tr && ref($tr) && $tr->isa('Bio::EnsEMBL::Transcript')) { unless($tr && ref($tr) && $tr->isa('Bio::EnsEMBL::Transcript')) {
throw("transcript arg must be Bio::EnsEMBL:::Transcript not [$tr]"); throw("transcript arg must be Bio::EnsEMBL:::Transcript not [$tr]");
} }
#convert exons coordinates to peptide coordinates #convert exons coordinates to peptide coordinates
my $tmp_exon = $self->transfer($tr->slice); my $tmp_exon = $self->transfer($tr->slice);
if (!$tmp_exon) { if (!$tmp_exon) {
throw("Couldn't transfer exon to transcript's slice"); throw("Couldn't transfer exon to transcript's slice");
} }
my @coords = my @coords =
$tr->genomic2pep($tmp_exon->start, $tmp_exon->end, $tmp_exon->strand); $tr->genomic2pep($tmp_exon->start, $tmp_exon->end, $tmp_exon->strand);
#filter out gaps #filter out gaps
@coords = grep {$_->isa('Bio::EnsEMBL::Mapper::Coordinate')} @coords; @coords = grep {$_->isa('Bio::EnsEMBL::Mapper::Coordinate')} @coords;
#if this is UTR then the peptide will be empty string #if this is UTR then the peptide will be empty string
my $pep_str = ''; my $pep_str = '';
if(scalar(@coords) > 1) {
throw("Error. Exon maps to multiple locations in peptide." . if(scalar(@coords) > 1) {
" Is this exon [$self] a member of this transcript [$tr]?"); my $coord = $self->_merge_ajoining_coords(\@coords);
} elsif(scalar(@coords) == 1) { if($coord) {
@coords = ($coord);
}
else {
my ($e_id, $tr_id) = ($self->stable_id(), $tr->stable_id());
throw("Error. Exon maps to multiple locations in peptide and those".
" locations are not continuous." .
" Is this exon [$e_id] a member of this transcript [$tr_id]?");
}
}
elsif(scalar(@coords) == 1) {
my $c = $coords[0]; my $c = $coords[0];
my $pep = $tr->translate; my $pep = $tr->translate;
#bioperl doesn't give back residues for incomplete codons #bioperl doesn't give back residues for incomplete codons
#make sure we don't subseq too far... #make sure we don't subseq too far...
my ($start, $end); my ($start, $end);
...@@ -1261,7 +1271,7 @@ sub peptide { ...@@ -1261,7 +1271,7 @@ sub peptide {
$start = ($c->start < $end) ? $c->start : $end; $start = ($c->start < $end) ? $c->start : $end;
$pep_str = $tr->translate->subseq($start, $end); $pep_str = $tr->translate->subseq($start, $end);
} }
return return
Bio::Seq->new( -seq => $pep_str, Bio::Seq->new( -seq => $pep_str,
-moltype => 'protein', -moltype => 'protein',
...@@ -1269,6 +1279,46 @@ sub peptide { ...@@ -1269,6 +1279,46 @@ sub peptide {
-id => $self->display_id ); -id => $self->display_id );
} }
=head2 _merge_ajoining_coords
Arg [1] : ArrayRef of Bio::EnsEMBL::Mapper::Coordinate objects
Example :
Description : Merges coords which are ajoining or overlapping
Returntype : Bio::EnsEMBL::Mapper::Coordinate or undef if it cannot happen
Exceptions : Exception if the cooords cannot be condensed into one location
Caller : internal
Status : Development
=cut
sub _merge_ajoining_coords {
my ($self, $coords) = @_;
my $okay = 1;
my $coord = shift @{$coords};
my $start = $coord->start();
my $last_end = $coord->end();
foreach my $other_coord (@{$coords}) {
if( ($last_end + 1) >= $other_coord->start() ) {
$last_end = $other_coord->end();
}
else {
$okay = 0;
last;
}
}
if(!$okay) {
return;
}
my $new_coord = Bio::EnsEMBL::Mapper::Coordinate->new(
$coord->id(), $start, $last_end, $coord->strand(), $coord->rank());
return $new_coord;
}
=head2 seq =head2 seq
......
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