Bio-EnsEMBL
view release on metacpan or search on metacpan
lib/Bio/EnsEMBL/DBSQL/BaseFeatureAdaptor.pm view on Meta::CPAN
until it satisfies -LIMIT or hits -MAX_RANGE. Useful if you don't know how far away the features
might be, or if dealing with areas of high feature density. In the case of Variation Features, it is
conceivable that a 2000 bp window might contain very many features, resulting in a slow and bloated
response, thus the ability to explore outward in smaller sections can be useful.
Returntype : Listref of [$feature,$distance]
=cut
sub fetch_all_by_outward_search {
my $self = shift;
my ($ref_feature, $respect_strand, $opposite_strand, $downstream, $upstream, $start_search_range,
$limit,$not_overlapping,$five_prime,$three_prime, $max_range) =
rearrange([qw(FEATURE SAME_STRAND OPPOSITE_STRAND DOWNSTREAM UPSTREAM RANGE LIMIT NOT_OVERLAPPING FIVE_PRIME THREE_PRIME MAX_RANGE)], @_);
my $factor = 1;
$limit ||= 1;
$start_search_range ||= 1000;
my $current_search_range = $start_search_range;
$max_range ||= 10000;
my @results;
while (scalar @results < $limit && $current_search_range < $max_range) {
$current_search_range = $start_search_range * $factor;
@results = @{
$self->fetch_all_nearest_by_Feature(-RANGE => $current_search_range,
-FEATURE => $ref_feature,
-SAME_STRAND => $respect_strand,
-OPPOSITE_STRAND => $opposite_strand,
-DOWNSTREAM => $downstream,
-UPSTREAM => $upstream,
-LIMIT => $limit,
-NOT_OVERLAPPING => $not_overlapping,
-FIVE_PRIME => $five_prime,
-THREE_PRIME => $three_prime,
)
};
$factor++;
}
return \@results;
}
=head2 fetch_all_nearest_by_Feature
lib/Bio/EnsEMBL/DBSQL/BaseFeatureAdaptor.pm view on Meta::CPAN
a -DOWN/UPSTREAM value and/or -NOT_OVERLAPPING to fulfil its most common application.
Returntype : Listref containing an Arrayref of Bio::EnsEMBL::Feature objects and the distance
[ [$feature, $distance] ... ]
Caller : general
=cut
sub fetch_all_nearest_by_Feature{
my $self = shift;
my ($ref_feature, $respect_strand, $opposite_strand, $downstream, $upstream, $search_range,$limit,$not_overlapping,$five_prime,$three_prime) =
rearrange([qw(FEATURE SAME_STRAND OPPOSITE_STRAND DOWNSTREAM UPSTREAM RANGE LIMIT NOT_OVERLAPPING FIVE_PRIME THREE_PRIME)], @_);
if ( !defined($search_range)) {
$search_range ||= 1000;
}
$limit ||= 1;
unless (defined($ref_feature) && $ref_feature->isa('Bio::EnsEMBL::Feature')) {
throw ('fetch_all_nearest_by_Feature method requires a valid Ensembl Feature object to operate');
}
throw ('Do not specify both -upstream and -downstream. This is the same as no arguments') if ($upstream && $downstream);
lib/Bio/EnsEMBL/DBSQL/BaseFeatureAdaptor.pm view on Meta::CPAN
}
if (($downstream && $ref_feature->strand == 1)|| ($upstream && $ref_feature->strand == -1)) {
if ($five_prime) {
@candidates = grep { ($_->strand == 1) ? $_->seq_region_start > $ref_feature->end : $_->seq_region_end > $ref_feature->end } @candidates;
} elsif ($three_prime) {
@candidates = grep { ($_->strand == 1) ? $_->seq_region_end > $ref_feature->end : $_->seq_region_start > $ref_feature->end } @candidates;
}
}
# Then sort and prioritise the candidates
my $finalists; # = [[feature, distance, centre-weighted distance, length, dbID],..]
$finalists = $self->select_nearest($ref_feature,\@candidates,$limit,$not_overlapping,$five_prime,$three_prime);
$finalists = [ map { [ splice @$_,0,2 ]} @$finalists ]; # Remove the ugly bits from the sight of users.
return $finalists;
}
=head2 select_nearest
Arg [1] : Bio::Ensembl::Feature, a Feature to find the nearest neighbouring feature to.
Arg [2] : Listref of Features to be considered for nearness.
Arg [3] : Integer, limited number of Features to return. Equally near features are all returned in spite of this limit
Arg [4] : Boolean, Overlapping prohibition. Overlapped Features are forgotten
Arg [5] : Boolean, use the 5' ends of the nearby features for distance calculation
Arg [6] : Boolean, use the 3' ends of the nearby features for distance calculation
Example : $feature_list = $feature_adaptor->select_nearest($ref_feature,\@candidates,$limit,$not_overlapping)
Description: Take a list of possible features, and determine which is nearest. Nearness is a
tricky concept. Beware of using the distance between Features, as it may not be the number you think
it should be.
Returntype : listref of Features ordered by proximity
Caller : BaseFeatureAdaptor->fetch_all_nearest_by_Feature
=cut
sub select_nearest {
my $self = shift;
my $ref_feature = shift;
my $candidates = shift;
my $limit = shift;
my $not_overlapping = shift;
my $five_prime = shift;
my $three_prime = shift;
# Convert circular coordinates to linear ones for distance calculation
my $ref_start = ($ref_feature->start < $ref_feature->end) ? $ref_feature->start : $ref_feature->start - $ref_feature->length; # Not ->end, in case circular
my $ref_end = $ref_feature->end;
my $ref_midpoint = $self->_compute_midpoint($ref_feature);
my $position_matrix = [];
my $shortest_distance;
my $adjusted_distance; #Â for when features overlap and we need another metric to order by
foreach my $neighbour (@$candidates) {
# rearrange coordinates into forward-stranded orientation, plus flatten circular features.
my $neigh_start = ($neighbour->seq_region_start < $neighbour->seq_region_end)
? $neighbour->seq_region_start : $neighbour->seq_region_start - $neighbour->length; # Not ->end, in case it is circular
my $neigh_midpoint = $self->_compute_midpoint($neighbour);
my $neigh_end = $neighbour->seq_region_end;
# discard overlaps early if not required.
next if ( $not_overlapping
&& (
( $neigh_start >= $ref_start && $neigh_end <= $ref_end )
|| ( $neigh_end <= $ref_end && $neigh_end >= $ref_start )
)
);
my @args = ($ref_start,$ref_midpoint,$ref_end,$neigh_start,$neigh_midpoint,$neigh_end,$neighbour->strand);
if ($five_prime || $three_prime) {
my $five = $five_prime ? 1 : 0; # Choose 5' or 3' end
($shortest_distance,$adjusted_distance) = $self->_compute_prime_distance($five,@args);
next unless ($shortest_distance || $adjusted_distance);
}
else {
($shortest_distance,$adjusted_distance) = $self->_compute_nearest_end(@args);
}
push @$position_matrix,[ $neighbour, $shortest_distance, $adjusted_distance, $neighbour->length, $neighbour->display_id] unless ($not_overlapping && $shortest_distance == 0);
}
# Order by distance, then centre-to-centre distance, then smallest feature first, then an arbitrary ID.
# $position_matrix looks like this:
# [ [ $feature, closest measure of distance, size, dbID ] ]
# Not proud of this line. It re-emits the original array in sorted order to favour better measures of order.
my @ordered_matrix = map { [$_->[0],$_->[1],$_->[2],$_->[3],$_->[4]] }
sort { abs($a->[1]) <=> abs($b->[1]) || abs($a->[2]) <=> abs($b->[2]) || $a->[3] <=> $b->[3] || $a->[4] <=> $b->[4]} @$position_matrix;
# Knock off unwanted hits. Equal distances are included, so more can be returned than asked for.
@ordered_matrix = $self->_discard_excess_features_from_matrix(\@ordered_matrix,$limit);
lib/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm view on Meta::CPAN
Status : Stable
=cut
sub store {
my $self = shift;
my $slice = shift;
my $seqref = shift;
my $not_dna = shift;
#
# Get all of the sanity checks out of the way before storing anything
#
if(!ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice'))) {
throw('Slice argument is required');
}
my $cs = $slice->coord_system();
lib/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm view on Meta::CPAN
throw("Slice must have end==seq_region_length");
}
my $sr_len = $slice->length();
my $sr_name = $slice->seq_region_name();
if($sr_name eq '') {
throw("Slice must have valid seq region name.");
}
if($cs->is_sequence_level() && !$not_dna) {
if(!$seqref) {
throw("Must provide sequence for sequence level coord system.");
}
if(ref($seqref) ne 'SCALAR') {
throw("Sequence must be a scalar reference.");
}
my $seq_len = length($$seqref);
if($seq_len != $sr_len) {
throw("Sequence length ($seq_len) must match slice length ($sr_len).");
lib/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm view on Meta::CPAN
$sth->bind_param(3,$cs->dbID,SQL_INTEGER);
$sth->execute();
my $seq_region_id = $self->last_insert_id('seq_region_id', undef, 'seq_region');
if(!$seq_region_id) {
throw("Database seq_region insertion failed.");
}
if($cs->is_sequence_level() && !$not_dna) {
#store sequence if it was provided
my $seq_adaptor = $db->get_SequenceAdaptor();
$seq_adaptor->store($seq_region_id, $$seqref);
}
#synonyms
if(defined($slice->{'synonym'})){
foreach my $syn (@{$slice->{'synonym'}} ){
$syn->seq_region_id($seq_region_id); # set the seq_region_id
$syn->adaptor->store($syn);
lib/Bio/EnsEMBL/RepeatMaskedSlice.pm view on Meta::CPAN
Exceptions : none
Caller : RawComputes (PredictionTranscript creation code).
Status : Stable
=cut
sub new {
my $caller = shift;
my $class = ref($caller) || $caller;
my ($logic_names, $soft_mask, $not_default_masking_cases) = rearrange(['REPEAT_MASK',
'SOFT_MASK',
'NOT_DEFAULT_MASKING_CASES'], @_);
my $self = $class->SUPER::new(@_);
$logic_names ||= [''];
if(ref($logic_names) ne 'ARRAY') {
throw("Reference to list of logic names argument expected.");
}
$self->{'repeat_mask_logic_names'} = $logic_names;
$self->{'soft_mask'} = $soft_mask;
$self->{'not_default_masking_cases'} = $not_default_masking_cases;
$self->{'not_default_masking_cases'} ||= {};
return $self;
}
=head2 repeat_mask_logic_names
Arg [1] : reference to list of strings $logic_names (optional)
Example : $rm_slice->repeat_mask_logic_name(['repeat_masker']);
Description: Getter/Setter for the logic_names of the repeats that are used
lib/Bio/EnsEMBL/RepeatMaskedSlice.pm view on Meta::CPAN
Status : Stable
=cut
sub soft_mask {
my $self = shift;
$self->{'soft_mask'} = shift if(@_);
return $self->{'soft_mask'} || 0;
}
=head2 not_default_masking_cases
Arg [1] : hash reference $not_default_masking_cases (optional, default is {})
The values are 0 or 1 for hard and soft masking respectively
The keys of the hash should be of 2 forms
"repeat_class_" . $repeat_consensus->repeat_class,
e.g. "repeat_class_SINE/MIR"
"repeat_name_" . $repeat_consensus->name
e.g. "repeat_name_MIR"
depending on which base you want to apply the not default masking either
the repeat_class or repeat_name. Both can be specified in the same hash
at the same time, but in that case, repeat_name setting has priority over
repeat_class. For example, you may have hard masking as default, and
you may want soft masking of all repeat_class SINE/MIR,
but repeat_name AluSp (which are also from repeat_class SINE/MIR)
Example : $rm_slice->not_default_masking_cases({"repeat_class_SINE/MIR" => 1,
"repeat_name_AluSp" => 0});
Description: Getter/Setter which is used to escape some repeat class or name from the default
masking in place.
Returntype : hash reference
Exceptions : none
Caller : seq() and subseq() methods
Status : Stable
=cut
sub not_default_masking_cases {
my $self = shift;
$self->{'not_default_masking_cases'} = shift if (@_);
return $self->{'not_default_masking_cases'};
}
=head2 seq
Arg [1] : none
Example : print $rmslice->seq(), "\n";
Description: Retrieves the entire repeat masked sequence for this slice.
See also the B<Bio::EnsEMBL::Slice> implementation of this
method.
Returntype : string
lib/Bio/EnsEMBL/RepeatMaskedSlice.pm view on Meta::CPAN
=cut
sub seq {
my $self = shift;
#
# get all the features
#
my $repeats = $self->_get_repeat_features($self);
my $soft_mask = $self->soft_mask();
my $not_default_masking_cases = $self->not_default_masking_cases();
#
# get the dna
#
my $dna = $self->SUPER::seq(@_);
#
# mask the dna
#
$self->_mask_features(\$dna,$repeats,$soft_mask,$not_default_masking_cases);
return $dna;
}
=head2 subseq
Arg [1] : none
Example : print $rmslice->subseq(1, 1000);
Description: Retrieves a repeat masked sequence from a specified subregion
of this slice. See also the B<Bio::EnsEMBL::Slice>
implementation of this method.
lib/Bio/EnsEMBL/RepeatMaskedSlice.pm view on Meta::CPAN
$sub_end = $seq_region_slice->length ;
}
$subslice = $seq_region_slice->sub_Slice($sub_start, $sub_end);
}
else {
$subslice = $subsequence_slice;
}
my $repeats = $self->_get_repeat_features($subslice);
my $soft_mask = $self->soft_mask();
my $not_default_masking_cases = $self->not_default_masking_cases();
my $dna = $subsequence_slice->SUPER::seq();
$subsequence_slice->_mask_features(\$dna,$repeats,$soft_mask,$not_default_masking_cases);
return $dna;
}
=head2 _get_repeat_features
Args [1] : Bio::EnsEMBL::Slice to fetch features for
Description : Gets repeat features for the given slice
Returntype : ArrayRef[Bio::EnsEMBL::RepeatFeature] array of repeats
=cut
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
=cut
sub _get_VariationAdaptor {
my ($self, $object_type, $dbtype) = @_;
# very important to do this defaulting since we *have* to assume the variation
# database is called 'variation'. Using the current group will not work because
# that will be something like 'core' (most likely), 'ensembl' or 'vega'.
$dbtype ||= 'variation';
# Also we do not care about Registry->get_db() for variation DBs
my $do_not_check_db = 1;
return $self->_get_Adaptor($object_type, $dbtype, $do_not_check_db);
}
=head2 _get_CoreAdaptor
Arg [1] : String object_type to retrieve an adaptor for
Arg [2] : String dbtype to search for the given adaptor in. Defaults to core
Description : Searches for the specified adaptor in the Registry and returns it. Otherwise
it will return nothing if the adaptor was not found
ReturnType : Bio::EnsEMBL::DBSQL::BaseAdaptor derrived instance (specific to core-like dbs)
Exceptions : none
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
This method will warn when adaptors are missing but will never through an
exception. It is up to the calling code to decide how to handle the unavailablity
of an adaptor.
ReturnType : Bio::EnsEMBL::DBSQL::BaseAdaptor derrived instance. Otherwise it returns nothing
Exceptions : none
=cut
sub _get_Adaptor {
my ($self, $object_type, $dbtype, $do_not_check_db) = @_;
if(!$object_type) {
warning('Object type is a required parameter');
return;
}
my $adaptor = $self->adaptor();
if(!$adaptor) {
warning("Cannot get a ${object_type} adaptor without a SliceAdaptor attached to this instance of ".ref($self));
return;
}
my $local_db = $adaptor->db();
my $species = $local_db->species();
#First we query for the DBAdaptor using get_db(). This is a deprecated method
#call but "special" adaptors can be registered via this method. We must
#consult here 1st to find the possible special adaptor
if(!$do_not_check_db && $dbtype) {
my $db = $registry->get_db($local_db, $dbtype);
if($db) {
# If we got a return then use this DBAdaptor's species name, group and the given object type.
# Special adaptors can have different species names
$adaptor = $registry->get_adaptor($db->species(), $db->group(), $object_type);
}
else {
#Otherwise just use the current species, dbtype and object type
$adaptor = $registry->get_adaptor($species, $dbtype, $object_type);
}
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
if (my $adaptor = $self->_get_CoreAdaptor('KaryotypeBand')) {
return $adaptor->fetch_all_by_Slice($self);
}
return [];
}
=head2 get_repeatmasked_seq
Arg [1] : listref of strings $logic_names (optional)
Arg [2] : int $soft_masking_enable (optional)
Arg [3] : hash reference $not_default_masking_cases (optional, default is {})
The values are 0 or 1 for hard and soft masking respectively
The keys of the hash should be of 2 forms
"repeat_class_" . $repeat_consensus->repeat_class,
e.g. "repeat_class_SINE/MIR"
"repeat_name_" . $repeat_consensus->name
e.g. "repeat_name_MIR"
depending on which base you want to apply the not default
masking either the repeat_class or repeat_name. Both can be
specified in the same hash at the same time, but in that case,
repeat_name setting has priority over repeat_class. For example,
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
or soft-masked when arg[2] = 1 (sequence in lowercase)
Will only work with database connection to get repeat features.
Returntype : Bio::EnsEMBL::RepeatMaskedSlice
Exceptions : none
Caller : general
Status : Stable
=cut
sub get_repeatmasked_seq {
my ($self,$logic_names,$soft_mask,$not_default_masking_cases) = @_;
return Bio::EnsEMBL::RepeatMaskedSlice->new
(-START => $self->{'start'},
-END => $self->{'end'},
-STRAND => $self->{'strand'},
-ADAPTOR => $self->adaptor(),
-SEQ => $self->{'seq'},
-SEQ_REGION_NAME => $self->{'seq_region_name'},
-SEQ_REGION_LENGTH => $self->{'seq_region_length'},
-COORD_SYSTEM => $self->{'coord_system'},
-REPEAT_MASK => $logic_names,
-SOFT_MASK => $soft_mask,
-NOT_DEFAULT_MASKING_CASES => $not_default_masking_cases);
}
=head2 _mask_features
Arg [1] : reference to a string $dnaref
Arg [2] : array_ref $repeats
reference to a list Bio::EnsEMBL::RepeatFeature
give the list of coordinates to replace with N or with
lower case
Arg [3] : int $soft_masking_enable (optional)
Arg [4] : hash reference $not_default_masking_cases (optional, default is {})
The values are 0 or 1 for hard and soft masking respectively
The keys of the hash should be of 2 forms
"repeat_class_" . $repeat_consensus->repeat_class,
e.g. "repeat_class_SINE/MIR"
"repeat_name_" . $repeat_consensus->name
e.g. "repeat_name_MIR"
depending on which base you want to apply the not default masking either
the repeat_class or repeat_name. Both can be specified in the same hash
at the same time, but in that case, repeat_name setting has priority over
repeat_class. For example, you may have hard masking as default, and
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
(soft masking). The reference to a dna string which is passed
is changed in place.
Returntype : none
Exceptions : none
Caller : seq
Status : Stable
=cut
sub _mask_features {
my ($self,$dnaref,$repeats,$soft_mask,$not_default_masking_cases) = @_;
$soft_mask = 0 unless (defined $soft_mask);
$not_default_masking_cases = {} unless (defined $not_default_masking_cases);
# explicit CORE::length call, to avoid any confusion with the Slice
# length method
my $dnalen = CORE::length($$dnaref);
REP:foreach my $old_f (@{$repeats}) {
my $f = $old_f->transfer( $self );
my $start = $f->start;
my $end = $f->end;
my $length = ($end - $start) + 1;
lib/Bio/EnsEMBL/Slice.pm view on Meta::CPAN
# to add the following, and the other commented line few lines below.
my $rc_class;
my $rc_name;
if ($f->isa('Bio::EnsEMBL::RepeatFeature')) {
$rc_class = "repeat_class_" . $f->repeat_consensus->repeat_class;
$rc_name = "repeat_name_" . $f->repeat_consensus->name;
}
my $masking_type;
$masking_type = $not_default_masking_cases->{$rc_class} if (defined $not_default_masking_cases->{$rc_class});
$masking_type = $not_default_masking_cases->{$rc_name} if (defined $not_default_masking_cases->{$rc_name});
$masking_type = $soft_mask unless (defined $masking_type);
if ($masking_type) {
$padstr = lc substr ($$dnaref,$start,$length);
} else {
$padstr = 'N' x $length;
}
substr ($$dnaref,$start,$length) = $padstr;
}
lib/Bio/EnsEMBL/Utils/VegaCuration/Translation.pm view on Meta::CPAN
foreach my $trans (@$transcripts) {
my $tsi = $trans->stable_id;
my $tID = $trans->dbID;
my $tname = $trans->get_all_Attributes('name')->[0]->value;
if($has_log_object){
$hidden_remak_ttributes=$log_object->species_params->get_attributes->{$tsi}->{'hidden_remark'};
}else{
$hidden_remak_ttributes=$trans->get_all_Attributes('hidden_remark');
}
foreach my $rem (@$hidden_remak_ttributes) {
if ($rem->value =~ /not_for_Vega/) {
#$support->log_verbose("Skipping transcript $tname ($tsi) since 'not_for_Vega'\n",1);
$log_object->_save_log('log_verbose', '', $gsi, '', $tsi, '', "Skipping transcript $tname ($tsi) since 'not_for_Vega'");
next TRANS;
}
}
# go no further if there is a ribosomal framshift attribute
foreach my $attrib (@{$trans->get_all_Attributes('_rib_frameshift')}) {
if ($attrib->value) {
$log_object->_save_log('log', '', $gsi, '', $tsi, '', "Skipping $tsi ($tname) since it has a ribosomal frameshift attribute");
next TRANS;
}
lib/Bio/EnsEMBL/Utils/VegaCuration/Translation.pm view on Meta::CPAN
#$support->log_warning("INTERNAL STOPS EXTERNAL: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons[$mod_date]:\nSequence = $orig_seq\nStops at $positions)\n\n");
$log_object->_save_log('log_warning', '', $gsi, 'TRANSCRIPT', $tsi, 'VQCT_internal_stop', "INTERNAL STOPS EXTERNAL: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons[$mod_date]: Sequence = $orig_seq Stops at $positions...
}
}
#...yes - check remarks
else {
my $flag_remark = 0; # 1 if word seleno has been used
my $flag_remark2 = 0; # 1 if existing remark has correct numbering
my $alabel = 'Annotation_remark- selenocysteine ';
my $alabel2 = 'selenocysteine ';
my $annot_stops;
my $remarks;
my $att;
#get both hidden_remarks and remarks
foreach my $remark_type ('remark','hidden_remark') {
if($has_log_object){
$att=$log_object->species_params->get_attributes->{$trans->stable_id}->{$remark_type};
}else{
$att=$trans->get_all_Attributes($remark_type)
}
foreach my $attrib ( @$att) {
push @{$remarks->{$remark_type}}, $attrib->value;
}
}
#parse remarks to check syntax for location of edits
while (my ($attrib,$remarks)= each %$remarks) {
foreach my $text (@{$remarks}) {
if ( ($attrib eq 'remark') && ($text=~/^$alabel(.*)/) ){
#$support->log_warning("seleno remark for $tsi stored as Annotation_remark not hidden remark) [$mod_date]\n");
$log_object->_save_log('log_warning', '', $gsi, '', $tsi, 'VQCT_wrong_selC_coord', "seleno remark for $tsi stored as Annotation_remark not hidden remark) [$mod_date]");
$annot_stops=$1;
}
elsif ($text =~ /^$alabel2(.*)/) {
my $maybe = $1;
if($maybe =~ /^\s*\d+(\s+\d+)*\s*$/) {
$annot_stops=$maybe;
} else {
$log_object->_save_log('log', '', $gene->stable_id, '', $tsi, '', "Maybe annotated stop in incorrect format, maybe just a remark that happens to begin '$alabel2'".
" -- might need to investigate: '$alabel2$maybe' [$mod_date]");
}
}
}
}
#check the location of the annotated edits matches actual stops in the sequence
my @annotated_stops;
if ($annot_stops){
my $i = 0;
foreach my $offset (split(/\s+/, $annot_stops)) {
#OK if it matches a known stop
if (
defined($found_stops[$i]) && defined($found_stops[$i]->[1]) && ($found_stops[$i]->[1] == $offset)) {
push @annotated_stops, $offset;
}
# catch old annotations where number was in DNA not peptide coordinates
elsif (defined($found_stops[$i]) && defined($found_stops[$i]->[1]) && (($found_stops[$i]->[1] * 3) == $offset)) {
$log_object->_save_log('log_warning', '', $gene->stable_id, 'DNA', $tsi, 'VQCT_wrong_selC_coord', "DNA: Annotated stop for transcript tsi ($tname) is in DNA not peptide coordinates) [$mod_date]");
}
# catch old annotations where number off by one
( run in 2.341 seconds using v1.01-cache-2.11-cpan-0a987023a57 )