File Coverage

blib/lib/Bio/RNA/SpliceSites/Scoring/MaxEntScan.pm

Criterion	Covered	Total	%
statement	13	15	86.6
branch			n/a
condition			n/a
subroutine	5	5	100.0
pod			n/a
total	18	20	90.0

package Bio::RNA::SpliceSites::Scoring::MaxEntScan;

use 5.008;
use strict;
use warnings;
use Carp;

#Data submodules; these names are from the original maxEntScan distribution.
use Bio::RNA::SpliceSites::Scoring::SpliceModels::me2x3acc;
use Bio::RNA::SpliceSites::Scoring::SpliceModels::me2x5;
use Bio::RNA::SpliceSites::Scoring::SpliceModels::splice5sequences;

require Exporter;

our @ISA = qw/ Exporter /;

my $functions = [ qw/ score5 score3 / ];
our %EXPORT_TAGS = ( 'all' => $functions , );
our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
our $VERSION = '0.04';

sub log2 {
  my $number = shift;
  log($number)/log(2);
}

sub is_scalar_reference {
  my $reference_to_validate = shift;
  ref( $reference_to_validate ) eq "SCALAR" ? 1 : 0;
}

sub is_kmer {
  my ( $reference_to_validate , $valid_length ) = @_;
  length $$reference_to_validate == $valid_length ? 1 : 0;
}

sub is_genetic_alphabet {
  my $reference_to_validate = shift;
  $$reference_to_validate =~ /^[ACTGactg]+$/ ? 1 : 0;
}

sub split_sequence {
  #Split the provided splice site sequence into the splice donor/acceptor dinucleotide and the concatenated remainder.
  my ( $reference_to_sequence , $splice_site_type ) = @_;
  my @sequence_array = split // , uc( $$reference_to_sequence );
  my ( $dinucleotide , $outer_portion );
  if ( $splice_site_type == 5 ) {
    $dinucleotide  =   join "" , @sequence_array[3..4]; #Positions 4 and 5 are the splice donor dinucleotide.
    $outer_portion = ( join "" , @sequence_array[0..2] )  . ( join "" , @sequence_array[5..8] );
  }
  elsif ( $splice_site_type == 3 ) {
    $dinucleotide  =   join "" , @sequence_array[18..19]; #Positions 19 and 20 are the splice acceptor dinucleotide.
    $outer_portion = ( join "" , @sequence_array[0..17] ) . ( join "" , @sequence_array[20..22] ); #Join nucleotides 1-18 and 21-23.
  }
  else {
    croak "Invalid type of splice site to split: must be either 5 or 3.\n";
  }
  return ( $dinucleotide , $outer_portion );
}

sub get_splice_5_sequence_matrix_value {
  my $outer_portion = shift;
  exists $Bio::RNA::SpliceSites::Scoring::SpliceModels::splice5sequences::table->{ $outer_portion } ?
    return $Bio::RNA::SpliceSites::Scoring::SpliceModels::splice5sequences::table->{ $outer_portion } : carp "Unable to find sequence matrix for key $outer_portion.\n";
}

sub get_splice_5_score_matrix_value {
  my $sequence_matrix_value = shift; #The term "matrix" was used in the original maxEntScan programming and is retained here for clarity.
  exists $Bio::RNA::SpliceSites::Scoring::SpliceModels::me2x5::table->[ $sequence_matrix_value ] ?
    return $Bio::RNA::SpliceSites::Scoring::SpliceModels::me2x5::table->[ $sequence_matrix_value ] : carp "Index out of score matrix range: $sequence_matrix_value\n";
}

sub score_consensus {
  my ( $dinucleotide , $type ) = @_;
  my %bgd =  ( 'A' => 0.27   , 'C' => 0.23   , 'G' => 0.23   , 'T' => 0.27 ); #Shared between each type of splice site.
  my ( %cons1 , %cons2 ); #Populate conditional to splice site type: donor or acceptor.

  if ( $type == 5 ) {
    return 15.7507349436393 if $dinucleotide eq 'GT'; #Short circuit for the perfect GT splice donor dinucleotide.
    %cons1 = ( 'A' => 0.004  , 'C' => 0.0032 , 'G' => 0.9896 , 'T' => 0.0032 );
    %cons2 = ( 'A' => 0.0034 , 'C' => 0.0039 , 'G' => 0.0042 , 'T' => 0.9884 );
  }
  elsif ( $type == 3 ) {
    #return ????? if $dinucleotide eq 'AG'; #Short circuit for perfect AG splice acceptor dinucleotide.
    %cons1 = ( 'A' => 0.9903 , 'C' => 0.0032 , 'G' => 0.0034 , 'T' => 0.003 );
    %cons2 = ( 'A' => 0.0027 , 'C' => 0.0037 , 'G' => 0.9905 , 'T' => 0.003 );
  }
  else {
    croak "Invalid type of splice site to score consensus for: must be either 5 or 3.\n";
  }

  my ( $first_nucleotide , $second_nucleotide ) = split // , $dinucleotide;
  return ( $cons1{$first_nucleotide} * $cons2{$second_nucleotide} ) / ( $bgd{$first_nucleotide} * $bgd{$second_nucleotide} );
}

sub score5 {
  my $sequence_reference = shift;
  #Validate argument:
  unless ( is_scalar_reference( $sequence_reference ) ) {
    carp "Not a scalar reference.\n";
    return 'invalid_invocation';
  }
  unless ( is_kmer( $sequence_reference , 9 ) ) {
    carp "Invalid 5'ss length: must be 9 nucleotides long.\n";
    return 'invalid_length';
  }
  unless ( is_genetic_alphabet( $sequence_reference ) ) {
    carp "Invalid alphabet: must be only [ACTGactg] with no 'n' nucleotides.\n";
    return 'invalid_alphabet';
  }

  my ( $dinucleotide , $outer_portion ) = split_sequence( $sequence_reference , 5 );

  #Compute the log2 of the product of the score_5_consensus() subroutine return for the entire sequence (left_side_of_product)
  #  and the me2x5 (score matrix) value for the sequence matrix value of the "outer portion" of the splice site, which is the heptamer with the donor dinucleotide removed
  #  from the splice site nonamer, which is assigned to the scalar variable $outer_portion

  my $score = log2( score_consensus( $dinucleotide , 5 ) * get_splice_5_score_matrix_value( get_splice_5_sequence_matrix_value( $outer_portion ) ) );
  return sprintf( "%.2f" , $score ); #2 decimals.
}

sub hash_seq {
  #Returns a hash key for a sequence as a 4-radix integer.
  #E.g. given sequence 'CAGAAGT', returns 4619 as a scalar.
  my $sequence = shift;
  $sequence=~ y/ACGT/0123/;
  my @sequence_array = split // , $sequence;
  my $sum = 0;
  my $end = length( $sequence ) - 1;

  my @four_radix = qw/ 1 4 16 64 256 1024 4096 16384 /;
  $sum += $sequence_array[$_] * $four_radix[ $end - $_ ] for 0 .. $end;
  return $sum;
}

sub get_max_ent_score {
  my ( $sequence , $table_ref ) = @_; #Table ref is a reference to an array of references to hashes prepared in the me2x3acc submodule.
  my @partial_score = ( $table_ref->[0]{ hash_seq( substr $sequence , 0  , 7 ) } ,
			$table_ref->[1]{ hash_seq( substr $sequence , 7  , 7 ) } ,
			$table_ref->[2]{ hash_seq( substr $sequence , 14 , 7 ) } ,
			$table_ref->[3]{ hash_seq( substr $sequence , 4  , 7 ) } ,
			$table_ref->[4]{ hash_seq( substr $sequence , 11 , 7 ) } ,
			$table_ref->[5]{ hash_seq( substr $sequence , 4  , 3 ) } ,
			$table_ref->[6]{ hash_seq( substr $sequence , 7  , 4 ) } ,
			$table_ref->[7]{ hash_seq( substr $sequence , 11 , 3 ) } ,
			$table_ref->[8]{ hash_seq( substr $sequence , 14 , 4 ) } );
  my $final_score = $partial_score[0] * $partial_score[1] * $partial_score[2] * $partial_score[3] * $partial_score[4] /
    ( $partial_score[5] * $partial_score[6] * $partial_score[7] * $partial_score[8] );
  return $final_score;
}

sub score3 {
  my $sequence_reference = shift; #Should be 23nt long.
  #Validate argument:
  unless ( is_scalar_reference( $sequence_reference ) ) {
    carp "Not a scalar reference.\n";
    return 'invalid_invocation';
  }
  unless ( is_kmer( $sequence_reference , 23 ) ) {
    carp "Invalid 3'ss length: must be 23nt long.\n";
    return 'invalid_length';
  }
  unless ( is_genetic_alphabet( $sequence_reference ) ) {
    carp "Invalid alphabet: must be only [ACTGactg] with no 'n' nucleotides.\n";
    return 'invalid_alphabet';
  }

  my ( $dinucleotide , $outer_portion ) = split_sequence( $sequence_reference , 3 );
  my $score = log2( score_consensus( $dinucleotide , 3 ) * get_max_ent_score( $outer_portion , $Bio::RNA::SpliceSites::Scoring::SpliceModels::me2x3acc::table ) );
  return sprintf( "%.2f" , $score ); #2 decimals.
}

1;

__END__