#610355590

Committed 22 Sep 2023 02:48PM UTC coverage: 83.241%. First build
Build # #610355590
Build Type
Pull #1040
travis-ci
Pull Request Pull Request #1040: Report consequences for each breakend
Run Details
134 of 134 new or added lines in 6 files covered. (100.0%)
22570 of 27114 relevant lines covered (83.24%)
1430.55 hits per line
Source File
Press 'n' to go to next uncovered line, 'b' for previous

87.48
/modules/Bio/EnsEMBL/Variation/Utils/VariationEffect.pm
=head1 LICENSE

Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
Copyright [2016-2025] EMBL-European Bioinformatics Institute

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file excepst in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=cut


=head1 CONTACT

 Please email comments or questions to the public Ensembl
 developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

 Questions may also be sent to the Ensembl help desk at
 <http://www.ensembl.org/Help/Contact>.

=cut

=head1 NAME

Bio::EnsEMBL::Variation::Utils::VariationEffect

=head1 DESCRIPTION

This module defines a set of predicate subroutines that check the effect of a
Bio::EnsEMBL::Variation::VariationFeature on some other Bio::EnsEMBL::Feature. 
All of these predicates take a VariationFeatureOverlapAllele as their first and
only argument and return a true or false value depending on whether the effect
being checked for holds or not. The link between these predicates and the 
specific effect is configured in the Bio::EnsEMBL::Variation::Utils::Config 
module and a list of OverlapConsequence objects that represent a link between,
for example, a Sequence Ontology consequence term, and the predicate that
checks for it is provided in the Bio::EnsEMBL::Variation::Utils::Constants
module. If you want to add a new consequence you should write a predicate in 
this module and then add an entry in the configuration file.

=cut

package Bio::EnsEMBL::Variation::Utils::VariationEffect;

use strict;
use warnings;

use base qw(Exporter);

our @EXPORT_OK = qw(overlap _intron_overlap _compare_seq_region_names within_feature within_cds MAX_DISTANCE_FROM_TRANSCRIPT within_intron stop_lost stop_retained start_lost frameshift $UPSTREAM_DISTANCE $DOWNSTREAM_DISTANCE $CHROMOSOME_SYNONYMS);

use constant MAX_DISTANCE_FROM_TRANSCRIPT => 5000;

our $UPSTREAM_DISTANCE = MAX_DISTANCE_FROM_TRANSCRIPT;
our $DOWNSTREAM_DISTANCE = MAX_DISTANCE_FROM_TRANSCRIPT;
our $CHROMOSOME_SYNONYMS = {};

#
# Interface with some of the module function XS reimplementation
# 
# If Bio::EnsEMBL::XS is installed, assign the function glob to
# the XS counterpart, otherwise assign to the original function
#
BEGIN {
  if (eval { require Bio::EnsEMBL::XS; 1 }) {
    *overlap = \&Bio::EnsEMBL::XS::Variation::Utils::VariationEffect::overlap;
  } else {
    *overlap = \&overlap_perl;
  }
}

# these two methods are perl implementations of the above Inline C
sub overlap_perl {
    my ( $f1_start, $f1_end, $f2_start, $f2_end ) = @_;
   
    return ( ($f1_end >= $f2_start) and ($f1_start <= $f2_end) );
}

sub _intron_overlap {
  my ($vf_start, $vf_end, $intron_start, $intron_end, $insertion) = @_;
  
  if(
           overlap($vf_start, $vf_end, $intron_start+2, $intron_start+7) ||
           overlap($vf_start, $vf_end, $intron_end-7,   $intron_end-2  ) ||
    overlap($vf_start, $vf_end, $intron_start-3, $intron_start-1) ||
           overlap($vf_start, $vf_end, $intron_end+1,   $intron_end+3  ) ||
           (
      $insertion && (
               $vf_start == $intron_start   ||
               $vf_end   == $intron_end     ||
               $vf_start == $intron_start+2 ||
               $vf_end   == $intron_end-2
      )
    )
  ) {
    return 1;
  }
  else {
    return 0;
  }
}

sub _compare_seq_region_names {
  my ($region1, $region2) = @_;
  if ($CHROMOSOME_SYNONYMS) {
    return 1 if
      (
        exists $CHROMOSOME_SYNONYMS->{$region1} &&
        grep /^$region2$/, keys %{$CHROMOSOME_SYNONYMS->{$region1}}
      ) || (
        exists $CHROMOSOME_SYNONYMS->{$region2} &&
        grep /^$region1$/, keys %{$CHROMOSOME_SYNONYMS->{$region2}}
      );
  }
  return lc $region1 eq lc $region2;
}

sub within_feature {
    my ($bvfoa, $feat, $bvfo, $bvf, $match_seq_region_names) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->feature;
    $match_seq_region_names ||= 0;

    my $cmp_chr = 1;
    if ($match_seq_region_names) {
      my $chr  = $bvf->{chr} || $bvf->{slice}->{seq_region_name};
      $cmp_chr = _compare_seq_region_names($chr, $feat->{slice}->{seq_region_name});
    }

    return $cmp_chr && overlap(
        $bvf->{start}, 
        $bvf->{end},
        $feat->{start}, 
        $feat->{end}
    );
}

sub partial_overlap_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->feature;
    
    return (
        within_feature(@_) and 
        (not complete_overlap_feature(@_)) and
        (($bvf->{end} > $feat->{end}) or ($bvf->{start} < $feat->{start}))
    );
}

sub complete_within_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->feature;
    
    return (
        ($bvf->{start} >= $feat->{start}) and 
        ($bvf->{end} <= $feat->{end})
    );
}

sub complete_overlap_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->feature;
    
    return ( 
        ($bvf->{start} <= $feat->{start}) and 
        ($bvf->{end} >= $feat->{end}) 
    );
}

sub _supporting_cnv_terms {
  #if variant is CNV, return class SO terms for its supporting variants
  my $bvf = shift;

  return if $bvf->class_SO_term(undef, 1) ne "copy_number_variation";
  return unless defined $bvf->structural_variation;

  my $support_vars  = $bvf->structural_variation->get_all_SupportingStructuralVariants;
  my @support_terms = map { $_->class_SO_term } @{$support_vars};
  return @support_terms;
}

sub deletion {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;
    
    # sequence variant will have alleles
    if($bvf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
        my ($ref_allele, $alt_allele) = _get_alleles(@_);
        return (
            (defined($ref_allele) && ($alt_allele eq '' || length($alt_allele) < length($ref_allele)) and $ref_allele) or
            $bvf->allele_string =~ /deletion/i
        );
    }
    
    # structural variant depends on class
    if($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {
      return (
          copy_number_loss(@_) or
          ($bvf->class_SO_term(undef, 1) eq 'deletion') or
          ($bvf->class_SO_term(undef, 1) =~ /deletion/i) or
          grep(/deletion/i, _supporting_cnv_terms($bvf))
      );
    }
    
    else { return 0; }
}

sub insertion {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;
    
    # sequence variant will have alleles
    if($bvf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
        my ($ref_allele, $alt_allele) = _get_alleles(@_);
        return (
            (defined($ref_allele) && ($ref_allele eq '' || length($alt_allele) > length($ref_allele)) and $alt_allele) or
            $bvf->allele_string =~ /insertion/i
        );
    }
    
    # structural variant depends on class
    if($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {
      my $class_SO_term = $bvf->class_SO_term(undef, 1);

      return (
          copy_number_gain(@_) or
          ($bvf->class_SO_term(undef, 1) eq 'insertion') or
          ($bvf->class_SO_term(undef, 1) =~ /insertion/i) or
          grep(/insertion/i, _supporting_cnv_terms($bvf))
      );
    }
    
    else { return 0; }
}

sub copy_number_gain {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;

    return (
        duplication(@_) or
        tandem_repeat(@_) or
        $bvf->class_SO_term(undef, 1) =~ /gain/i or
        grep(/gain/i, _supporting_cnv_terms($bvf))
    );
}

sub copy_number_loss {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;
    
    return (
      $bvf->class_SO_term(undef, 1) =~ /loss/i or
      grep(/loss/i, _supporting_cnv_terms($bvf))
    );
}

sub duplication {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;
    
    return (
        (
            ($bvf->class_SO_term(undef, 1) eq 'duplication') or
            ($bvf->class_SO_term(undef, 1) =~ /duplication/i) or
            grep(/duplication/i, _supporting_cnv_terms($bvf))
        ) and
        (not tandem_repeat(@_))
    );
}

sub tandem_repeat {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf ||= $bvfoa->base_variation_feature;
    
    # for sequence variants, check sequence vs ref
    if($bvf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
        my ($ref_allele, $alt_allele) = _get_alleles(@_);
        
        return 0 unless $ref_allele and $alt_allele;
        return 0 unless
            length($alt_allele) > length($ref_allele) and
            length($alt_allele) % length($ref_allele) == 0;
        
        my $copies = length($alt_allele) / length($ref_allele);
        
        return $alt_allele eq $ref_allele x $copies;
    }
    
    # structural variant depends on class
    if($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {
        return (
            ($bvf->class_SO_term(undef, 1) eq 'tandem_duplication') or
            ($bvf->class_SO_term(undef, 1) eq 'tandem_repeat') or
            ($bvf->class_SO_term(undef, 1) =~ /tandem/i) or
            grep(/tandem/i, _supporting_cnv_terms($bvf))
        );
    }
}

sub chromosome_breakpoint {
  my ($bvfoa, $feat, $bvfo, $bvf) = @_;
  $bvf ||= $bvfoa->base_variation_feature;

  if ($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {
    return (
      ($bvf->class_SO_term(undef, 1) eq 'chromosome_breakpoint') or
      ($bvf->class_SO_term(undef, 1) =~ /chromosome_breakpoint/i)
    );
  }
}

sub feature_ablation {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;
    
    return (complete_overlap_feature($bvfoa, $feat, $bvfo, $bvf) and deletion(@_));
}

sub feature_amplification {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;
    
    return (complete_overlap_feature($bvfoa, $feat, $bvfo, $bvf) and copy_number_gain(@_));
}

sub feature_elongation {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;
    
    return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele');
    
    return (
        within_cdna(@_) and
        complete_within_feature($bvfoa, $feat, $bvfo, $bvf) and
        (copy_number_gain(@_) or insertion(@_))
    );
}

sub feature_truncation {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele');
    
    if(chromosome_breakpoint(@_)) {
        return 1 if within_feature($bvfoa, $feat, $bvfo, $bvfoa->breakend, 1);
    }

    # require transcripts (but not other feature types) to be within cDNA
    return 0 if $feat->isa('Bio::EnsEMBL::Transcript') and not within_cdna(@_);

    return (
        (
            partial_overlap_feature($bvfoa, $feat, $bvfo, $bvf) or
            complete_within_feature($bvfoa, $feat, $bvfo, $bvf)
        ) and
        (
            copy_number_loss(@_) or deletion(@_)
        )
    );
}

sub protein_altering_variant{

    ## check in protein
    my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);

    return 0 unless defined $ref_pep && defined $alt_pep;

    ## don't assign if child term appropriate

    return 0 if  length($alt_pep) eq length($ref_pep);       # synonymous_variant(@_);  missense_variant(@_);
    return 0 if  $ref_pep =~/^\*/  || $alt_pep =~/^\*/;      # stop lost/ gained/ retained
    return 0 if  $alt_pep =~/^\Q$ref_pep\E|\Q$ref_pep\E$/;   # inframe_insertion(@_);

    return 0 if inframe_deletion(@_);  
    return 0 if start_lost(@_);
    return 0 if frameshift(@_);

    return 1;
}

#sub transcript_fusion {
#    #my ($bvfoa, $feat, $bvfo, $bvf) = @_;
#    #my $bvf   = $bvfoa->base_variation_feature;
#    
#    return 0;
#    
#    #my $transcripts = $bvf->_get_overlapping_Transcripts();
#}

sub _before_start {
    my ($bvf, $feat, $dist) = @_;
    
    return ( ($bvf->{end} >= ($feat->{start} - $dist)) and 
        ($bvf->{end} < $feat->{start}) );
}

sub _after_end {
    my ($bvf, $feat, $dist) = @_;
    return ( ($bvf->{start} <= ($feat->{end} + $dist)) 
            and ($bvf->{start} > $feat->{end}) );
}

sub _upstream {
    my ($bvf, $feat, $dist) = @_;
    $dist += get_max_shift_length($bvf);
    return $feat->strand == 1 ? 
        _before_start($bvf, $feat, $dist) : 
        _after_end($bvf, $feat, $dist);
}

sub _downstream {
    my ($bvf, $feat, $dist) = @_;
    $dist += get_max_shift_length($bvf);
    return $feat->strand == 1 ? 
        _after_end($bvf, $feat, $dist) : 
        _before_start($bvf, $feat, $dist);
}

sub get_max_shift_length {
  my ($bvf) = shift;
  my $max_shift_length = 0;
  foreach my $hash (@{$bvf->{tva_shift_hashes}}){
    $max_shift_length = $hash->{shift_length} if $hash->{shift_length} > $max_shift_length;
  }
  return $max_shift_length;
}

#package Bio::EnsEMBL::Variation::TranscriptVariationAllele;

sub upstream {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return _upstream($bvf, $feat, $UPSTREAM_DISTANCE);
}

sub downstream {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvf  ||= $bvfoa->base_variation_feature;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return _downstream($bvf, $feat, $DOWNSTREAM_DISTANCE);
}

sub affects_transcript {
    my ($bvf, $tran) = @_;
    
    return 0 unless $tran->isa('Bio::EnsEMBL::Transcript');
    
    return overlap(
        $bvf->{start}, 
        $bvf->{end},
        $tran->{start} - 5000, 
        $tran->{end} + 5000
    );
}

sub within_transcript {
    return within_feature(@_);
}

sub within_nmd_transcript {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return ( within_transcript(@_) and ($feat->biotype eq 'nonsense_mediated_decay') );
}

sub within_coding_gene {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return ( within_transcript(@_) and $feat->translation );
}

sub coding_transcript_variant {
    return ( (not coding_unknown(@_)) and complete_overlap_feature(@_) and within_coding_gene(@_) );
}

sub within_non_coding_gene {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;

    return ( within_transcript(@_) and (not $feat->translation) and (not within_mature_miRNA(@_)) and (not non_coding_exon_variant(@_)) );
}

sub non_coding_exon_variant {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    
    return 0 if complete_overlap_feature(@_) or $feat->translation or within_mature_miRNA(@_);
    
    # get overlapped exons
    # this may include some non-overlapping ones in the case of transcripts with frameshift introns
    # so we double check with overlap()
    my $exons = $bvfo->_overlapped_exons;
    
    if(scalar grep {overlap($bvf->{start}, $bvf->{end}, $_->{start}, $_->{end})} @$exons) {
        return 1;
    }
    else {
        return 0;
    }
}

sub within_miRNA {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    
    # don't call this for now
    
    return 0;
    $feat ||= $bvfoa->base_variation_feature_overlap->feature;
    
    return ( ($feat->biotype eq 'miRNA') and within_transcript(@_) );
}

sub within_mature_miRNA {
  my ($bvfoa, $feat, $bvfo, $bvf) = @_;
  $bvfo ||= $bvfoa->base_variation_feature_overlap;
  $bvf  ||= $bvfo->base_variation_feature;
  $feat ||= $bvfo->feature;

  return 0 unless ( ($feat->biotype eq 'miRNA') and within_transcript(@_) );

  foreach my $attribute(@{ $feat->get_all_Attributes('miRNA') }) {

    if (defined $attribute && $attribute->value =~ /(\d+)-(\d+)/) {
      for my $coord ($bvfo->_mapper->cdna2genomic($1, $2)) {
        if ($coord->isa('Bio::EnsEMBL::Mapper::Coordinate')) {
          if (overlap(
              $bvf->seq_region_start(), 
              $bvf->seq_region_end(), 
              $coord->start, 
              $coord->end) ) {
            return 1;
          }
        }
      }
    }
  }

  return 0;
}

sub donor_splice_site {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;

    my $ie = $bvfoa->_intron_effects($feat, $bvfo, $bvf);
    
    return $feat->strand == 1 ? 
        $ie->{start_splice_site} :
        $ie->{end_splice_site};
}

sub acceptor_splice_site {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;

    my $ie = $bvfoa->_intron_effects($feat, $bvfo, $bvf);
    
    return $feat->strand == 1 ? 
        $ie->{end_splice_site} :
        $ie->{start_splice_site};
}

sub essential_splice_site {
    return ( acceptor_splice_site(@_) or donor_splice_site(@_) );
}

sub splice_donor_5th_base_variant {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    my $ie = $bvfoa->_intron_effects($feat, $bvfo, $bvf);

    return $feat->strand == 1 ? $ie->{fifth_base_splice_site} : $ie->{fifth_base_splice_site_reverse};
}

sub splice_donor_region_variant {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    my $ie = $bvfoa->_intron_effects($feat, $bvfo, $bvf);
    
    return 0 if splice_donor_5th_base_variant(@_);
    return $feat->strand == 1 ? $ie->{donor_region_splice_site} : $ie->{donor_region_splice_site_reverse};
}

sub splice_polypyrimidine_tract_variant  {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    my $ie = $bvfoa->_intron_effects($feat, $bvfo, $bvf);

    return $feat->strand == 1 ? $ie->{polypyrimidine_splice_site} : $ie->{polypyrimidine_splice_site_reverse};
}

sub splice_region {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    
    return 0 if donor_splice_site(@_);
    return 0 if acceptor_splice_site(@_);
    return 0 if essential_splice_site(@_);
    return 0 if splice_donor_region_variant(@_);
    return 0 if splice_donor_5th_base_variant(@_);
    
    return $bvfoa->_intron_effects($feat, $bvfo, $bvf)->{splice_region};
}

sub within_intron {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;

    return $bvfoa->_intron_effects($feat, $bvfo, $bvf)->{intronic};
}

sub within_cds {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    $bvf  ||= $bvfo->base_variation_feature;
    
    my $cds_coords = $bvfo->cds_coords;
    
    if (@$cds_coords > 0) {
        for my $coord (@$cds_coords) {
            if ($coord->isa('Bio::EnsEMBL::Mapper::Coordinate')) {
                if ($coord->end > 0 && $coord->start <= length($bvfo->_translateable_seq)) { 
                    return 1;
                }
            }
        }
    }

    # we also need to check if the vf is in a frameshift intron within the CDS

    if (defined $feat->translation &&
        $bvfoa->_intron_effects($feat, $bvfo, $bvf)->{within_frameshift_intron}) {
 
        return overlap(
            $bvf->{start}, 
            $bvf->{end}, 
            $feat->coding_region_start,
            $feat->coding_region_end,
        );
    }
        
    return 0;
}

sub within_cdna {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    
    my $cdna_coords = $bvfo->cdna_coords;
    
    if (@$cdna_coords > 0) {
        for my $coord (@$cdna_coords) {
            if ($coord->isa('Bio::EnsEMBL::Mapper::Coordinate')) {
                if ($coord->end > 0 && $coord->start <= $feat->length) {
                    return 1;
                }
            }
        }
    }
    
    # we also need to check if the vf is in a frameshift intron within the cDNA

    if ($bvfoa->_intron_effects->{within_frameshift_intron}) {
        return within_transcript(@_); 
    }
    
    return 0;
}

sub _before_coding {
    my ($bvf, $tran) = @_;
    return 0 unless defined $tran->translation;
    
    my $bvf_s  = $bvf->{start};
    my $bvf_e  = $bvf->{end};
    my $t_s    = $tran->{start};
    my $cds_s  = $tran->coding_region_start;
    
    # we need to special case insertions just before the CDS start
    if ($bvf_s == $bvf_e+1 && $bvf_s == $cds_s) {
        return 1;
    }
   
    return overlap($bvf_s, $bvf_e, $t_s, $cds_s-1);    
}

sub _after_coding {
    my ($bvf, $tran) = @_;
    return 0 unless defined $tran->translation;
    
    my $bvf_s  = $bvf->{start};
    my $bvf_e  = $bvf->{end};
    my $t_e    = $tran->{end};
    my $cds_e  = $tran->coding_region_end;
    
    # we need to special case insertions just after the CDS end
    if ($bvf_s == $bvf_e+1 && $bvf_e == $cds_e) {
        return 1;
    }
    
    return overlap($bvf_s, $bvf_e, $cds_e+1, $t_e);
}

sub within_5_prime_utr {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    $bvf  ||= $bvfo->base_variation_feature;

    my $five_prime_of_coding = 
        $feat->strand == 1 ? 
            _before_coding($bvf, $feat) : 
            _after_coding($bvf, $feat);
    
    return ( $five_prime_of_coding and within_cdna(@_) );
}

sub within_3_prime_utr {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $feat ||= $bvfo->feature;
    $bvf  ||= $bvfo->base_variation_feature;
    
    my $three_prime_of_coding = 
        $feat->strand == 1 ? 
            _after_coding($bvf, $feat) : 
            _before_coding($bvf, $feat);
    
    return ( $three_prime_of_coding and within_cdna(@_) );
}

sub complex_indel {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $bvf  ||= $bvfo->base_variation_feature;
   
    # pass the no_db flag to var_class to ensure we don't rely on the database for it 
    # as it may not have been set at this stage in the pipeline
    my $class = $bvf->var_class(1);

    return 0 unless $class =~ /insertion|deletion|indel/;

    return @{ $bvfo->cds_coords } > 1;
}

sub _get_peptide_alleles {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{_get_peptide_alleles})) {
        my @alleles = ();

        #return () if frameshift(@_);

        if(
            defined $bvfoa &&
            (my $alt_pep = $bvfoa->peptide) &&        
            (my $ref_pep = _get_ref_pep(@_))
        ) {
            $ref_pep = '' if $ref_pep eq '-';
            $alt_pep = '' if $alt_pep eq '-';
        
            @alleles = ($ref_pep, $alt_pep);
        }

        $cache->{_get_peptide_alleles} = \@alleles;
    }

    return @{$cache->{_get_peptide_alleles}};
}

sub _get_ref_pep {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    my $ref_tva = $bvfo->get_reference_TranscriptVariationAllele;

    ## The shift hash has to be added to the reference tva before calculating the relevant peptide string.
    ## This is because we don't precalculate the shift_hash for reference tvas for speed and for
    ## variants with multiple alternate alleles with potentially different shift lengths.
    $ref_tva->{shift_hash} = $bvfoa->{shift_hash} if (defined($bvfoa->{shift_hash}));
    return $ref_tva->peptide;
}

sub _get_codon_alleles {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    
    return () if frameshift(@_);

    my $alt_codon = $bvfoa->codon;
    
    return () unless defined $alt_codon;

    $bvfo ||= $bvfoa->base_variation_feature_overlap;    
    my $ref_codon = $bvfo->get_reference_TranscriptVariationAllele->codon;
    
    return () unless defined $ref_codon;
    
    $ref_codon = '' if $ref_codon eq '-';
    $alt_codon = '' if $alt_codon eq '-';
    
    return ($ref_codon, $alt_codon);
}

sub _get_alleles {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    
    my $ref_tva = $bvfo->get_reference_VariationFeatureOverlapAllele;
    
    return () unless defined ($ref_tva);
    
    my $ref_allele = $ref_tva->variation_feature_seq;
    my $alt_allele = $bvfoa->variation_feature_seq;
    
    return () unless defined($ref_allele) and defined($alt_allele);
    
    $ref_allele = '' if $ref_allele eq '-';
    $alt_allele = '' if $alt_allele eq '-';
    
    return ($ref_allele, $alt_allele);
}

sub start_lost {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{start_lost})) {

        # default
        $cache->{start_lost} = 0;

        return 0 unless _overlaps_start_codon(@_);

        $bvfo ||= $bvfoa->base_variation_feature_overlap;
        $feat ||= $bvfo->feature;
        $bvf  ||= $bvfo->base_variation_feature;
        
        # sequence variant
        if($bvfo->isa('Bio::EnsEMBL::Variation::TranscriptVariation')) {
            return $cache->{start_lost} = 1 if _ins_del_start_altered(@_) && !(inframe_insertion(@_) || inframe_deletion(@_));
            return $cache->{start_lost} = 1 if _inv_start_altered(@_);
            my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
        
            return 0 unless $ref_pep;
            return 0 unless $alt_pep && $alt_pep ne 'X';
            
            # allow for introducing additional bases that retain start codon e.g. atg -> aCGAtg
            $cache->{start_lost} = (
                ($bvfo->translation_start == 1) and
                ($alt_pep !~ /\Q$ref_pep\E$/) and 
                ($alt_pep !~ /^\Q$ref_pep\E/)
            );
        }
        
        # structural variant
        elsif($bvfo->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariation')) {        
            my ($tr_crs, $tr_cre) = ($feat->coding_region_start, $feat->coding_region_end);
            return 0 unless defined($tr_crs) && defined($tr_cre);
            
            if($feat->strand == 1) {
                $cache->{start_lost} = overlap($tr_crs, $tr_crs + 2, $bvf->{start}, $bvf->{end});
            }
            else {
                $cache->{start_lost} = overlap($tr_cre-2, $tr_cre, $bvf->{start}, $bvf->{end});
            }
        }
        
        else {
            return 0;
        }
    }

    return $cache->{start_lost};
}

sub _inv_start_altered {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;   
    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};
    unless(exists($cache->{inv_start_altered})) {
        $cache->{inv_start_altered} = 0;

        return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele');
        return 0 unless $bvfoa->seq_is_unambiguous_dna();
        return 0 unless _overlaps_start_codon(@_);

        $bvfo ||= $bvfoa->base_variation_feature_overlap;

        # get cDNA coords
        my ($cdna_start, $cdna_end) = ($bvfo->cdna_start, $bvfo->cdna_end);
        return 0 unless $cdna_start && $cdna_end;

        # make and edit UTR + translateable seq
        my $translateable = $bvfo->_translateable_seq();
        my $utr = $bvfo->_five_prime_utr();
        return 0 unless $utr;
        my $utr_and_translateable = ($utr ? $utr->seq : '').$translateable;
        my $shifting_offset = defined($bvfoa->{shift_hash}) ? $bvfoa->{shift_hash}->{shift_length} : 0;
        $cdna_start += $shifting_offset;
        $cdna_end += $shifting_offset;

        return 0 if($cdna_end > length($utr_and_translateable));

        my $vf_feature_seq = $bvfoa->feature_seq;
        $vf_feature_seq = '' if $vf_feature_seq eq '-';
        my $atg_start = length($utr->seq);

        substr($utr_and_translateable, $cdna_start - 1, ($cdna_end - $cdna_start) + 1) = $vf_feature_seq;
        my $new_sc = substr($utr_and_translateable, $atg_start, 3);

        return $cache->{inv_start_altered} = 1 if $new_sc ne 'ATG';
    }

    return $cache->{inv_start_altered};
}

sub start_retained_variant {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $bvf  ||= $bvfo->base_variation_feature;
    # structural variants don't have an allele string 
    return 0 if ($bvf->allele_string && ($bvf->allele_string eq 'COSMIC_MUTATION' || $bvf->allele_string eq 'HGMD_MUTATION'));

    my $pre = $bvfoa->_pre_consequence_predicates;

    return 0 unless _overlaps_start_codon(@_);
    if ($pre->{snp}) {
        return !_snp_start_altered(@_);
    } else { 
        return !_ins_del_start_altered(@_);
    }
}

sub _overlaps_start_codon {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{overlaps_start_codon})) {
        $cache->{overlaps_start_codon} = 0;

        $bvfo ||= $bvfoa->base_variation_feature_overlap;
        $feat ||= $bvfo->feature;
        return 0 if grep {$_->code eq 'cds_start_NF'} @{$feat->get_all_Attributes()};

        my ($cdna_start, $cdna_end) = ($bvfo->cdna_start_unshifted, $bvfo->cdna_end_unshifted);
        my $shifting_offset = defined($bvfoa->{shift_hash}) ? $bvfoa->{shift_hash}->{shift_length} : 0;
        $cdna_start += $shifting_offset;
        $cdna_end += $shifting_offset;
        return 0 unless $cdna_start && $cdna_end;
        
        $cache->{overlaps_start_codon} = overlap(
            $cdna_start, $cdna_end,
            $feat->cdna_coding_start, $feat->cdna_coding_start + 2
        );
    }

    return $cache->{overlaps_start_codon};
}

sub _snp_start_altered {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{snp_start_altered})) {
        $cache->{snp_start_altered} = 1;

        return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele');
        return 0 unless $bvfoa->seq_is_unambiguous_dna();

        $bvfo ||= $bvfoa->base_variation_feature_overlap;

        # get cDNA coords
        my ($cdna_start, $cdna_end) = ($bvfo->cdna_start, $bvfo->cdna_end);
        return 0 unless $cdna_start && $cdna_end;

        # make and edit UTR + translateable seq
        my $translateable = $bvfo->_translateable_seq();
        my $utr = $bvfo->_five_prime_utr();
        my $utr_and_translateable = ($utr ? $utr->seq : '').$translateable;

        my $vf_feature_seq = $bvfoa->feature_seq;
        $vf_feature_seq = '' if $vf_feature_seq eq '-';

        substr($utr_and_translateable, $cdna_start - 1, ($cdna_end - $cdna_start) + 1) = $vf_feature_seq;
        
        # get the aa now in place of prevoius start codon and see it is ATG
        # it can happen for non-ATG variant (e.g - gene WT1, transcript ENST00000332351)
        my $aa_replacing_start_codon = substr($utr_and_translateable, 0 - length($translateable), 3);
        $cache->{snp_start_altered} = 0 if $aa_replacing_start_codon eq "ATG";
    }

    return $cache->{snp_start_altered};
}

sub _ins_del_start_altered {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{ins_del_start_altered})) {
        $cache->{ins_del_start_altered} = 0;

        return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele');
        return 0 unless $bvfoa->seq_is_unambiguous_dna();
        return 0 unless _overlaps_start_codon(@_);

        my $pre = $bvfoa->_pre_consequence_predicates;
        return 0 unless $pre->{increase_length} || $pre->{decrease_length};

        $bvfo ||= $bvfoa->base_variation_feature_overlap;

        # get cDNA coords
        my ($cdna_start, $cdna_end) = ($bvfo->cdna_start, $bvfo->cdna_end);
        return 0 unless $cdna_start && $cdna_end;
        
        # make and edit UTR + translateable seq
        my $translateable = $bvfo->_translateable_seq();
        my $utr = $bvfo->_five_prime_utr();
        my $utr_and_translateable = ($utr ? $utr->seq : '').$translateable;

        my $vf_feature_seq = $bvfoa->feature_seq;
        $vf_feature_seq = '' if $vf_feature_seq eq '-';

        substr($utr_and_translateable, $cdna_start - 1, ($cdna_end - $cdna_start) + 1) = $vf_feature_seq;

        # check if still retain start
        if ($utr) {
            my $atg_start = length($utr->seq);
            my $new_sc = substr($utr_and_translateable, $atg_start, 3);
            my $new_utr = substr($utr_and_translateable, 0, length($utr->seq));
            return $cache->{ins_del_start_altered} if ($new_utr eq $utr->seq && $new_sc eq 'ATG');
        }

        # sequence shorter, we know it has been altered
        return $cache->{ins_del_start_altered} = 1 if length($utr_and_translateable) < length($translateable);

        $cache->{ins_del_start_altered} = $translateable ne substr($utr_and_translateable, 0 - length($translateable));
    }

    return $cache->{ins_del_start_altered};
}

sub synonymous_variant {
    my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
    
    return 0 unless $ref_pep;

    return ( ($alt_pep eq $ref_pep) and (not stop_retained(@_) and not start_retained_variant(@_) and ($alt_pep !~ /X/) and ($ref_pep !~ /X/)) );
}

sub missense_variant {
    my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
    
    return 0 unless defined $ref_pep;
    
    return 0 if start_lost(@_);
    return 0 if stop_lost(@_);
    return 0 if stop_gained(@_);
    return 0 if partial_codon(@_);
    return 0 if stop_retained(@_); # added because a variant can not be stop_retained and misense 
    
    # also need to check that the translated sequence is not the same i think . 
    return ( $ref_pep ne $alt_pep ) && ( length($ref_pep) == length($alt_pep)  );
}

sub inframe_insertion {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $bvf  ||= $bvfo->base_variation_feature;

    # sequence variant
    if($bvf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
        my ($ref_codon, $alt_codon) = _get_codon_alleles(@_);

        return 0 if start_lost(@_);
        return 0 unless defined $ref_codon;
        return 0 unless ( length($alt_codon) > length ($ref_codon) );

        my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);

        return 0 unless defined($ref_pep) && defined($alt_pep);

        # not an inframe insertion if the inserted AA is before the start codon
        # we can use start_retained to check this
        return 0 if start_retained_variant(@_) && $alt_pep =~ /\Q$ref_pep\E$/;

        # if we have a stop codon in the alt peptide
        # trim off everything after it
        # this allows us to detect inframe insertions that retain a stop
        $alt_pep =~ s/\*.+/\*/;

        return 1 if ($alt_pep =~ /^\Q$ref_pep\E/) || ($alt_pep =~ /\Q$ref_pep\E$/);

    }
    
    # structural variant
    elsif($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {

        # TO BE DONE, NO WAY OF KNOWING WHAT SEQUENCE IS INSERTED YET
        return 0;
        
        # must be an insertion
        return 0 unless insertion(@_);
        
        my $cds_coords = $bvfo->cds_coords;
        
        if(scalar @$cds_coords == 1) {
            
            # wholly within exon
            if($cds_coords->[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) {
                1;
            }
        }
        
        # variant partially overlaps
        else {
            return 0;
        }
    }
    
    else {
        return 0;
    }
}

sub inframe_deletion {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $bvf  ||= $bvfo->base_variation_feature;
    
    # sequence variant
    if($bvf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
        return 0 if partial_codon(@_);

        my ($ref_codon, $alt_codon) = _get_codon_alleles(@_);
        my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
        return 0 unless defined $ref_codon;
        return 0 unless length($alt_codon) < length ($ref_codon);
        
        # simple string match
        return 1 if ($ref_codon =~ /^\Q$alt_codon\E/) || ($ref_codon =~ /\Q$alt_codon\E$/);

        # check for internal match
        ($ref_codon, $alt_codon) = @{Bio::EnsEMBL::Variation::Utils::Sequence::trim_sequences($ref_codon, $alt_codon)};
        
        # if nothing remains of $alt_codon,
        # then it fully matched a part in the middle of $ref_codon
        return length($alt_codon) == 0 && length($ref_codon) % 3 == 0;
    }
    
    # structural variant
    elsif($bvf->isa('Bio::EnsEMBL::Variation::StructuralVariationFeature')) {
        
        # must be a deletion
        return 0 unless deletion(@_);
        
        my $cds_coords = $bvfo->cds_coords;
        my $exons      = $bvfo->_exons;
        
        # in exon
        return (
           scalar @$cds_coords == 1 and
           $cds_coords->[0]->isa('Bio::EnsEMBL::Mapper::Coordinate') and
           scalar grep {complete_within_feature($bvfoa, $_, $bvfo, $bvf)} @$exons and
           $bvf->length() % 3 == 0
        );
    }
    
    else {
        return 0;
    }
}

sub stop_gained {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{stop_gained})) {
        $cache->{stop_gained} = 0;
        
        ## check for inframe insertion before stop 
        return 0 if stop_retained(@_);

        my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
        
        return 0 unless defined $ref_pep;

        $cache->{stop_gained} = ( ($alt_pep =~ /\*/) and ($ref_pep !~ /\*/) );
    }

    return $cache->{stop_gained};
}

sub stop_lost {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};
    
    unless(exists($cache->{stop_lost})) {
        $cache->{stop_lost} = 0;

        $bvfo ||= $bvfoa->base_variation_feature_overlap;
        $bvf  ||= $bvfo->base_variation_feature;
        $feat ||= $bvfo->feature;
        
        # sequence variant
        if($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele')) {
            
            # special case frameshift
    #        if(frameshift(@_)) {
    #          my $ref_pep = _get_ref_pep(@_);
    #          return $ref_pep && $ref_pep =~ /\*/;
    #        }
            
            my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);
            if(defined($ref_pep) && defined($alt_pep)) {
                $cache->{stop_lost} = ( ($alt_pep !~ /\*/) and ($ref_pep =~ /\*/) );
            }
            else {
                $cache->{stop_lost} = _ins_del_stop_altered(@_);
            }
        }
        
        # structural variant
        elsif($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele')) {
            return 0 unless deletion(@_);
            
            my ($tr_crs, $tr_cre) = ($feat->coding_region_start, $feat->coding_region_end);
            return 0 unless defined($tr_crs) && defined($tr_cre);
            
            if($feat->strand == 1) {
                $cache->{stop_lost} = overlap($tr_cre-2, $tr_cre, $bvf->{start}, $bvf->{end});
            }
            else {
                $cache->{stop_lost} = overlap($tr_crs, $tr_crs + 2, $bvf->{start}, $bvf->{end});
            }
        }
        
        else {
            return 0;
        }
    }

    return $cache->{stop_lost};
}

sub stop_retained {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    return 0 if partial_codon(@_);
    return 0 if stop_lost(@_);

    unless(exists($cache->{stop_retained})) {
        $bvfo ||= $bvfoa->base_variation_feature_overlap;
        $bvf  ||= $bvfo->base_variation_feature;
        # structural variants don't have an allele string
        return 0 if ($bvf->allele_string && ($bvf->allele_string eq 'COSMIC_MUTATION' || $bvf->allele_string eq 'HGMD_MUTATION'));

        $cache->{stop_retained} = 0;

        $bvfo ||= $bvfoa->base_variation_feature_overlap;

        my $pre = $bvfoa->_pre_consequence_predicates;

        my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);

        if(defined($alt_pep) && $alt_pep ne '') {
         
          ## handle inframe insertion of a stop just before the stop (no ref peptide)
          $cache->{stop_retained} = ref_eq_alt_sequence(@_);
        }
        else {
            $cache->{stop_retained} = ($pre->{increase_length} || $pre->{decrease_length}) && _overlaps_stop_codon(@_) && !_ins_del_stop_altered(@_);
        }

    }
    
    return $cache->{stop_retained};
}

sub ref_eq_alt_sequence {
   my ($bvfoa, $feat, $bvfo, $bvf) = @_; 
   
   $bvfo ||= $bvfoa->base_variation_feature_overlap;
   my $ref_seq = $bvfo->_peptide;

   my $mut_seq = $ref_seq;
   my $tl_start = $bvfo->translation_start;
   my $tl_end = $bvfo->translation_end;
   
   my ($ref_pep, $alt_pep) = _get_peptide_alleles(@_);

   return 0 if $ref_pep eq "X" && $alt_pep eq "X"; # this is to account for incomplete coding terminal;

   # this is to account for synonymous variant if $ref_pep eq $alt_pep 
   # as there is no resulting change to the amino acid sequence, it is not stop_retained
   #return 0 if $ref_pep ne "*" && $alt_pep ne "*" && $ref_pep eq $alt_pep;

   # this is a logic from the former logic 
   return 1 if ($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele') && defined($ref_seq) && $tl_start > length($ref_seq) && $alt_pep =~ /^\*/);

   substr($mut_seq, $tl_start-1, $tl_end - $tl_start + 1) = $alt_pep; # creating a mutated sequence from the ref sequence. 

   my $mut_substring = substr($mut_seq, 0, length($ref_seq)); # getting a substring up to the length of the ref sequence for comparison from index 0 to the length of the ref seq;
   
   my $final_stop = substr($mut_seq, length($ref_seq)) if length($ref_seq) < length($mut_seq); # getting the length of the $mut_seq from the length of the ref_seq to the end 
   
   my $final_stop_length = length($final_stop) if defined($final_stop) ne '';
   
   # 1 is if the ref_pep and the first letter of the alt_pep is the same and the alt_pep has * in it 
   # 2 is the ref_seq eq $mut_substring and the final stop length is less than 3
   # 3 is * in ref_pep and the same index position exists for both the ref and alt pep
   return 1 if ( ($ref_pep eq substr($alt_pep, 0, 1) && $alt_pep =~ /\*/) ||
       ($ref_seq eq $mut_substring && defined($final_stop_length) && $final_stop_length < 3) || ( $ref_pep =~ /\*/ && (index($ref_pep, "*") + 1 == index($alt_pep, "*") + 1) ));
   return 0;
}

sub _overlaps_stop_codon {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{overlaps_stop_codon})) {
        $cache->{overlaps_stop_codon} = 0;

        $bvfo ||= $bvfoa->base_variation_feature_overlap;
        $feat ||= $bvfo->feature;
        return 0 if grep {$_->code eq 'cds_end_NF'} @{$feat->get_all_Attributes()};

        my ($cdna_start, $cdna_end) = ($bvfo->cdna_start, $bvfo->cdna_end);
        return 0 unless $cdna_start && $cdna_end;
        
        $cache->{overlaps_stop_codon} = overlap(
            $cdna_start, $cdna_end,
            $feat->cdna_coding_end - 2, $feat->cdna_coding_end
        );
    }

    return $cache->{overlaps_stop_codon};
}

sub _ins_del_stop_altered {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{ins_del_stop_altered})) {
        $cache->{ins_del_stop_altered} = 0;

        return 0 if $bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele');
        return 0 unless $bvfoa->seq_is_unambiguous_dna();
        return 0 unless _overlaps_stop_codon(@_);

        my $pre = $bvfoa->_pre_consequence_predicates;
        return 0 unless $pre->{increase_length} || $pre->{decrease_length};

        $bvfo ||= $bvfoa->base_variation_feature_overlap;

        # get cDNA coords and CDS start
        my ($cdna_start, $cdna_end, $cds_start) = ($bvfo->cdna_start, $bvfo->cdna_end, $bvfo->cds_start);
        return 0 unless $cdna_start && $cdna_end && $cds_start;

        # make and edit UTR + translateable seq
        my $translateable = $bvfo->_translateable_seq();
        my $utr = $bvfo->_three_prime_utr();

        my $utr_and_translateable = $translateable.($utr ? $utr->seq : '');

        my $vf_feature_seq = $bvfoa->feature_seq;
        $vf_feature_seq = '' if $vf_feature_seq eq '-';

        # use CDS start to anchor the edit
        # and cDNA coords to get the length (could use VF, but have already retrieved cDNA coords)
        substr($utr_and_translateable, $cds_start - 1, ($cdna_end - $cdna_start) + 1) = $vf_feature_seq;

        # new sequence shorter, we know it has been altered
        return $cache->{ins_del_stop_altered} = 1 if length($utr_and_translateable) < length($translateable);

        # now we need the codon from the new seq at the equivalent end pos from translateable
        # and to translate it to check if it is still a stop
        my $pep = Bio::Seq->new(
            -seq        => substr($utr_and_translateable, length($translateable) - 3, 3),
            -moltype    => 'dna',
            -alphabet   => 'dna',
        )->translate(
            undef, undef, undef, $bvfo->_codon_table
        )->seq;
        $cache->{ins_del_stop_altered} = !($pep && $pep eq '*');
    }

    return $cache->{ins_del_stop_altered};
}

sub frameshift {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;

    # sequence variant
    if($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele')) {

        return 0 if partial_codon(@_);
        return 0 if stop_retained(@_);
    
        return 0 unless defined $bvfo->cds_start && defined $bvfo->cds_end;
        
        my $var_len = $bvfo->cds_end - $bvfo->cds_start + 1;
    
        my $allele_len = $bvfoa->seq_length;
    
        # if the allele length is undefined then we can't call a frameshift
    
        return 0 unless defined $allele_len;
    
        return abs( $allele_len - $var_len ) % 3;
    }
    
    # structural variant
    elsif($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele')) {
        my $exons = $bvfo->_exons;
        
        return (
            (
                deletion(@_) or
                copy_number_loss(@_)
            ) and
            scalar grep {complete_within_feature($bvfoa, $_, $bvfo, $bvf)} @$exons and
            $bvf->length % 3 != 0
        );
        
        # TODO INSERTIONS
    }
    
    else {
        return 0;
    }
}

sub partial_codon {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;

    # use cache for this method as it gets called a lot
    my $cache = $bvfoa->{_predicate_cache} ||= {};

    unless(exists($cache->{_partial_codon})) {

        # default
        $cache->{_partial_codon} = 0;

        $bvfo ||= $bvfoa->base_variation_feature_overlap;
   
        return 0 unless defined $bvfo->translation_start;

        my $cds_length = length $bvfo->_translateable_seq;


        my $codon_cds_start = ($bvfo->translation_start * 3) - 2;

        my $last_codon_length = $cds_length - ($codon_cds_start - 1);
        
        $cache->{_partial_codon} = ( $last_codon_length < 3 and $last_codon_length > 0 );
    }

    return $cache->{_partial_codon};
}

sub coding_unknown {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    $bvfo ||= $bvfoa->base_variation_feature_overlap;
    $bvf  ||= $bvfo->base_variation_feature;

    return 0 if complete_overlap_feature(@_);

    # sequence variant
    if($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptVariationAllele')) {
        return (
            within_cds(@_) and (
                (not $bvfoa->peptide) or
                (not _get_peptide_alleles(@_)) or
                ($bvfoa->peptide =~ /X/) or
                ((_get_peptide_alleles(@_))[0] =~ /X/)
            ) and (
                not (
                    frameshift(@_) or inframe_deletion(@_) or protein_altering_variant(@_) or
                    start_retained_variant(@_) or start_lost(@_) or
                    stop_retained(@_) or stop_lost(@_)
                )
            )
        );
    }
    
    # structural variant
    elsif($bvfoa->isa('Bio::EnsEMBL::Variation::TranscriptStructuralVariationAllele')) {
        return (within_cds(@_) and not(inframe_insertion(@_) or inframe_deletion(@_) or frameshift(@_)));
    }
    
    else {
        return 0;
    }
}

#package Bio::EnsEMBL::Variation::RegulatoryFeatureVariationAllele;

sub within_regulatory_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    return within_feature(@_);
}

#package Bio::EnsEMBL::Variation::ExternalFeatureVariationAllele;

sub within_external_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    return (within_feature(@_) and (not within_miRNA_target_site(@_)));
}

#sub within_miRNA_target_site {
#    my $efva = shift;
#
#    my $fset = $efva->variation_feature_overlap->feature->feature_set;
#
#    return ($fset && $fset->name eq 'miRanda miRNA targets');
#}

#package Bio::EnsEMBL::Variation::MotifFeatureVariationAllele;

#sub within_motif_feature {
#    my $mfva = shift;
#    return (
#        within_feature($mfva) and
#        !increased_binding_affinity($mfva) and
#        !decreased_binding_affinity($mfva) 
#    );
#}

sub within_motif_feature {
    my ($bvfoa, $feat, $bvfo, $bvf) = @_;
    return within_feature(@_);
}

#sub increased_binding_affinity {
#    my $mfva = shift;
#    my $change = $mfva->binding_affinity_change;
#    return (within_feature($mfva) and (defined $change) and ($change > 0));
#}
#
#sub decreased_binding_affinity {
#    my $mfva = shift;
#    my $change = $mfva->binding_affinity_change;
#    return (within_feature($mfva) and (defined $change) and ($change < 0));
#}

sub contains_entire_feature {
    my $vfo = shift;

    my $bvf  = $vfo->base_variation_feature;
    my $feat = $vfo->feature;

    return ( ($bvf->{start} <= $feat->{start}) && ($bvf->{end} >= $feat->{end}) ); 
}

1;
Ensembl / ensembl-variation / #610355590

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous
