#!/usr/bin/perl -w
# Use_model
#    VERSION: Version 1 (16 June 2004)
#    PURPOSE:  Applies Markov model to sets of gene sequences
#         Outputs scores of fitness to model
#    INPUT FILES: 
#         Model file:
#              Hash made by MakeMarkov, stored by store method of Storable
#              Sums are stored in hash with key ending with chr(255)
#              Program will accept model of any order
#         Sequence file:
#              FastA formated DNA sequence
#    OUTPUT FILES: Tab-delimited file, with following quantitites:
#         mean-coordinate: coordinate in middle of moving window
#         Markov-score:  SOMETHING NEEDS TO GO HERE
#                Score is increasingly positive for sequences increasingly 
#                   close to sequences of training set

############## LIBRARIES AND PRAGMAS ################

  use strict;
  use Storable;
  use FastA_module;

#################### CONSTANTS ######################

  my $true = 1;	  	        # Perl doesn't have logical variables so
  my $false = ! $true;          #    we simulate them
  my $LF = "\n";                # Linefeed
  my $tab = "\t";               # Tab
  my $sum_symbol = chr(255);    # Key representing sum of elements of hash
  my $fill_char = chr(1);       # Used to pad keys to full length
  my $window_size = 200;        # Number of nucleotides used to calculate score
  my $window_increment = 25;    # Number of nucleotides by which to shift window

#################### VARIABLES ######################

  my %model;                    # Markov Model for all text
  my %nuc_frequencies;          # Frequency in sequence of each nucleotide
  my %model_bg;
  my $order;		        # Order of Markov Model used
  my $N;                        # Number of sequences in training set
  my $number_of_keys;           # Number of letter combinations in model
  my $number_of_windows;  	# Number of windows considered in making the model

  my $header;                   # Header line in current FastA-formatted orf
  my $sequence;                 # Given sequence 

  my %bg_frequency;             # nucleotide frequencies of overall sequences
  my %log_q;                    # $log_q{$nucleotide}[$position] is the log of the
                                #    frequency of the nucleotide at that position
  my %log_p;                    # $log_p{$nucleotide} is the background frequency
                                #    of the nucleotide
  my $score;                    # Score representing fit of gene with model 
                                #    based on training set
  my @scores;                   # Set of scores derived from window moving 
                                #    through sequence

###################### FILES ########################

  my $model_file_name = 'model-Tery.dat';
  %model = Get_Model($model_file_name);        # Restore Markov model

  my $input_file_name = 'te3022-region.nt';
  ($header, $sequence) = Read_FastA_sequence ($input_file_name);

  my $output_file_name = 'te3022-scores.txt';
  open OUTPUT, ">$output_file_name" or die "Can't open $output_file_name: $!\n";

################### MAIN PROGRAM ####################
####  Calculate nucleotide frequencies (may or may not be used)
####  Calculate background frequencies (log_p) based on ...
####  Calculate transition frequencies (log_q) based on model 
####  Go through given sequence, bit by bit, calculating score for each bit,
####     and save the scores

   %nuc_frequencies = Calc_nuc_frequencies ($sequence);
   %log_p = Calc_background_frequency();
   %log_q = Calc_motif_frequency(%model);
   @scores = Calculate_set_of_Markov_scores ($sequence);

   Save_scores (@scores);

#################### SUBROUTINES ####################

##### GET_MODEL (file name)
#     Retrieves model stored in given file as hash
#     Uses retrieve method of module Storable
#     Calls for calculation of statistics 
#        (only one, the order of the model, is used later)

  sub Get_Model {
     my ($file_name) = @_;
     my $model_reference = retrieve $file_name;
     my %model = %$model_reference;
     print $LF, "MODEL from $file_name:", $LF;
     $order = Stats_of_model(%model);
     return %model;
  }

##### STATS_OF_MODEL (hash)
#     Determines order of retrieved model and number of sequences in training set
#     Strategy: 

  sub Stats_of_model {
     my (%model) = @_;
     my $biggest_length = 0;
     my $current_order;
     my $number_of_instances = 0;
     my $number_of_keys = scalar (keys %model);
     my $number_of_sequences = 0;

     foreach my $key (keys %model) {
        if (substr ($key, -1, 1) eq $sum_symbol) 
           {$number_of_instances = $number_of_instances + $model{$key}}
        elsif (length ($key) == 1) 
           {$number_of_sequences = $number_of_sequences + $model{$key}}
        elsif (length ($key) > $biggest_length) 
           {$biggest_length = length ($key)}
     }
     $current_order = $biggest_length - 1;

     print "   order $current_order", $LF;
     print "   $number_of_sequences sequences in training set", $LF;
     print "   $number_of_keys keys and $number_of_instances instances";
     print $LF, $LF;

     return $current_order;
  }

                ########## SCORING SUBROUTINES ##########

##### CALC_NUC_FREQUENCIES
#     Counts nucleotides in given sequence
#     Calculates and returns nucleotide frequencies

  sub Calc_nuc_frequencies {
     my ($seq) = @_;         # Given sequence
     my %nuc_sum;            # Temporary hash to count nucleotides
     my %frequency;          # Frequency of each nucleotide
     my $count;              # Temporary count of specific nucleotide 
                             #     in specific sequence
     my $nucleotide;         # Current nucleotide 

     foreach $nucleotide ("A", "C", "G", "T") {
        $count = ($sequence =~ s/$nucleotide/$nucleotide/g);
                             # Perl magic counts nucleotide in sequence
        $nuc_sum{$nucleotide} += $count;
        $nuc_sum{$sum_symbol} += $count;
     }

     print "BACKGROUND FREQUENCIES:", $LF;
     foreach $nucleotide ("A", "C", "G", "T") {
        $frequency{$nucleotide} = $nuc_sum{$nucleotide}/$nuc_sum{$sum_symbol};   
        print "   $nucleotide: ";
        printf "%5.3f$LF", $frequency{$nucleotide};
     }
     return %frequency;
  } 

##### CALC_BACKGROUND_FREQUENCY
#     Counts nucleotides in all sequences
#     Calculates nucleotide frequencies
#     Calculates log of frequencies, used in later calculations

  sub Calc_background_frequency {
     my %log_p;

     ### SOMETHING NEEDS TO GO HERE ###

     return %log_p;
  } 

##### CALC_MOTIF_FREQUENCY
#     Calculates adjusted frequency (q) of nucleotide at last position 
#        of Markov chain, expressed as log(q)
    
  sub Calc_motif_frequency {
     my (%model) = @_;
     my %log_adj_freq;
     my $given_letters;    # First letters in letter combination
     my $predicted_letter; # Last letter in letter combination
     my $letters_key;      # Full letter combination
     my $instances;        # Number of times letters occur
     my $sum;              # Number of times given letters appear, followed by
                           #   any letter
     my $sum_key;          # Key accessing sum 

     foreach my $key (keys %model) {
        $given_letters = substr ($key, 0, length($key)-1);
        $predicted_letter = substr ($key, -1, 1);
        $sum_key = $given_letters . $sum_symbol;
        $sum = $model {$sum_key};
        $letters_key = $given_letters . $predicted_letter;
        $instances = $model{$letters_key};

        if (not $predicted_letter eq $sum_symbol) {
           $log_adj_freq{$key} =
                log( ### SOMETHING NEEDS TO GOES HERE ### )
        }
     }
     return %log_adj_freq;
  }
  
##### CALC_MARKOV_SCORE (sequence)
#     Calculates score by comparing probability of arriving at given sequence
#        by model frequencies vs background nucleotide frequencies
#     Runs a window through the sequence, where the size of the window is the
#        order of the model

  sub Calc_Markov_score {
     my ($seq) = @_;        # Sequence of which to calculate score
     my $letters;           # Letters (length order + 1) under consideration
     my $given_letters;     # First letters
     my $predicted_letter;  # Last letter
     my $score = 0;         # Will be sum of scores for each letter group
 
    foreach my $position (0 .. length($seq) - ($order+1)) {    
                   # Moves order+1th nucleotide through sequence
        $given_letters = substr ($seq, $position, $order);
        $predicted_letter = substr ($seq, $position + $order, 1);
        $letters = substr ($seq, $position, $order+1);
        if (defined($log_q{$letters})) 
          {$score += $log_q{$letters} - $log_p{$predicted_letter}}
        else 
          {print "WARNING: $letters not found in model!", $LF;
           $score += log(1/$number_of_windows) - $log_p{$predicted_letter}}
     }
     return $score;
  }

##### CALCULATE_SET_OF_MARKOV_SCORES (sequence)
#     Runs a window the size of $window_size through sequence
#     Calculates and stores Markov score for each window
#     Returns set of scores and associated coordinate

  sub Calculate_set_of_Markov_scores {
    my ($sequence) = @_;
    my $window;          # Window of size $window_size moving through sequence
    my $mean_coordinate; # Middle of window
    my @scores;          # Set of scores, one for each window

    for (my $coordinate = 0; 
         $coordinate <= length ($sequence) - $window_size; 
         $coordinate = $coordinate + $window_increment) {
       $window = substr( $sequence, $coordinate, $window_size);
       $mean_coordinate = $coordinate + ($window_size/2);
       push @scores, [$mean_coordinate, Calc_Markov_score ($window)];
    }
    return @scores;
  }

##### SAVE_SCORES (scores)
#     Outputs list of scores to a tab-delimited file

  sub Save_scores {
    my (@scores) = @_;

    foreach my $score (@scores) {
       print OUTPUT join ($tab, @$score), $LF;
    }
    
    close OUTPUT;
    print $LF, "OUTPUT:", $LF;
    print "   Scores saved: ", Size(@scores), $LF;
    print "   Window size: $window_size nucleotides", $LF;
    print "   Distance between coordinates: $window_increment", $LF;

  }

##### SIZE (array)
#     Returns size of array (synonymous with scalar)
  sub Size {return scalar(@_)}  