#!/usr/bin/perl -w
# Run_orfs_through_model
#    VERSION: Version 1 (1 December 2003)
#    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 of chr(255)
#              Program will accept model of any order
#         Orfs file:
#              FastA formated DNA sequences, any number
#    OUTPUT FILES: Tab-delimited file, with following quantitites:
#         score: Log of probability of gene sequence occuring according to model, 
#                   vs according to background nucleotide frequencies
#                Score is increasingly positive for sequences increasingly close 
#                   to sequences of training set
#         orf_name: First field of FastA header
#         orf_length: Number of nucleotides in FastA sequence
#         orth_eval: E-value of closest cyanobacterial ortholog 
#                   or "none" if no ortholog
#         orf_descr: Remainder of FastA header

############## 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 $very_small_number = 10**-200;

#################### VARIABLES ######################

  my %model;                    # Markov Model for all text
  my $order;		        # Order of Markov Model used
  my $N;                        # Number of sequences in training set
  my $elements;                 # Number of windows considered in making the model

  my $number_of_orfs;           # Number of orfs in input file of orfs
  my $orf_number;               # Position of current gene in input file of orfs
  my $header;                   # Header line in current FastA-formatted orf
  my @orf_info;                 # (Extracted from header) Name, length, Eval of ortholog, annotation
  my $sequence;                 # Sequence of current gene

  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

###################### FILES ########################

  my $model_file_name = '6803PHXModel.dat';
  Get_Model($model_file_name);        # Restore Markov model

  my $orf_file_name = 'C:\1\genome\6803\6803orfs.txt';

  print "Reading orfs from $orf_file_name...";
  $number_of_orfs = Read_FastA_sequences($orf_file_name);
  print $LF, "$number_of_orfs orfs read", $LF, $LF;
  
  my $output_file_name = 'C:\temp\6803orfsVsPHX4.txt';
  open OUTPUT, ">$output_file_name" or die "Can't open $output_file_name: $!\n";

################### MAIN PROGRAM ####################
#  Extract statistics from model (used in later calculations)
#  Calculate background frequencies (log_p) based on nucleotide frequencies and 
#     frequencies of nucleotides based on model (log_q)
#  Then go through each orf in file, calculating and outputing score

   ($order, $N) = Stats_of_model(%model);
   $elements = Elements_of_model(%model);

   print "Model of order $order read from $model_file_name", $LF;
   print "based on $N sequences in training set and $elements elements", $LF, $LF;

   print "Calculating background frequencies...", $LF;
   calc_background_frequency();
   calc_motif_frequency();
   print $LF;
   print "Press any key to continue";
   <STDIN>;

   foreach $orf_number (1 .. $number_of_orfs) {   # Go through all orfs of Synechocystis
     ($header, $sequence) = Get_sequence_info ($orf_number);
     @orf_info = Parse_orf_annotation($header);
     $score = Calc_Markov_score($sequence);  # Calculate fit with Markov model
     Output_orf();     
   }

#################### SUBROUTINES ####################

##### GET_MODEL (file name)
#     Retrieves model stored in given file as hash
#     Uses retrieve method of module Storable
  sub Get_Model {
     my ($model_file_name) = @_;
     my $model_reference = retrieve $model_file_name;
     %model = %$model_reference;
  }

##### STATS_OF_MODEL (hash)
#     Determines order of retrieved model and number of sequences in training set
#     Strategy: 
#        Goes through each hash element consisting solely of begin marks: 
#              $hash{""}{""}{""}...
#           until it reaches nucleotides:
#              $hash{""}{""}{""}...{"A"}
#        With each layer, the order is incremented by one
#        Last layer recognized because it contains scalars, not arrays of hashes
#        Sets $N = number of elements of training set, determined as the sum
#           of the elements of the hash beginning {""}{""}{""}...
#     Returns the order and $N
  sub Stats_of_model {
     my (%hash) = @_;
     my $current_order = 0;
     my $N;
     while (defined $hash{""}) {
        %hash = %{$hash{""}};
        $current_order = $current_order + 1;
     }
     $N = $hash{$sum_symbol};
     return ($current_order, $N);
  }

##### ELEMENTS_OF_MODEL (model)
#     Recursive routine that strips off each layer of a hash
#     Last layer recognized because it contains scalars, not arrays of hashes
#     Adds up all sums in last layer (equal to every Markov chain considered by
#        the model)
#     Returns the sum of sums
  sub Elements_of_model {
     my (%hash) = @_;
     my @keys = sort keys %hash;   
     my $key;
     if (not defined $elements) {$elements = 0};
     foreach $key (@keys) {
        if (ref($hash{$key}) eq "HASH") {
           Elements_of_model( %{$hash{$key}} );
        } else {
           $elements = $elements + $hash{$sum_symbol};
           last;
        }
     }
     return $elements;
  }

##### PARSE_ORF_ANNOTATION (header line)
#     Extract various information from header line of FastA-formatted file
#     Returns orf name, orf length, e-value of ortholog (or "none"), description of orf
  sub Parse_orf_annotation {
     my ($annotation) = @_;
     $annotation =~ /(\w+)\s+\((\d+).+\)\s+((<-)(\d+)|(\d+)(->))\s+((ORTH:(\S+))|\s*)\s*(.*)$/;
          # example:  ssl1234 (555 nt) <-123456 ORTH:3e-30 Protein kinase
     my $orf_name = $1;
     my $orf_length = $2;
     my $orf_start = ($5 or $6);
     my $orf_direction = ($4 or $7);
     my $orth_eval = ($10 or "none");
     my $orf_descr = $11;
     return ($orf_name, $orf_length, $orth_eval, $orf_descr);
  }

##### OUTPUT_ORF
#     Output score and orf info to screen and to file
  sub Output_orf {
     printf "%6.1f", $score;
     print $tab, join($tab, @orf_info), $LF;
     printf OUTPUT "%6.1f", $score;
     print OUTPUT $tab, join($tab, @orf_info), $LF;
  }

                ########## SCORING SUBROUTINES ##########

##### CALC_BACKGROUND_FREQUENCY
#     Counts nucleotides in all sequences
#     Calculates nucleotide frequencies
#     Calculates log of frequencies, used in later calculations
  sub calc_background_frequency {
     my %nuc_sum;            # Temporary hash to count nucleotides
     my $count;              # Temporary count of specific nucleotide in specific sequence
     my $header;             # Returned by Get_sequence_info (ignored)
     my $sequence;           # Sequence of specific orf
     my $nucleotide;         # Current nucleotide 
     my $orf;                # Current number of orf

     foreach $orf (1 .. $number_of_orfs) {
        ($header, $sequence) = Get_sequence_info($orf);
        foreach $nucleotide ("A", "C", "G", "T") {
           $count = ($sequence =~ s/$nucleotide/$nucleotide/g);
           $nuc_sum{$nucleotide} += $count;
           $nuc_sum{$sum_symbol} += $count;
        }
     }

     print "BACKGROUND FREQUENCIES:", $LF;
     foreach $nucleotide ("A", "C", "G", "T") {
        $bg_frequency{$nucleotide} = $nuc_sum{$nucleotide}/$nuc_sum{$sum_symbol};   
        $log_p{$nucleotide} = log($bg_frequency{$nucleotide});
        print "   $nucleotide: ";
        printf "%5.3f$LF", $bg_frequency{$nucleotide};
     }
  } 

##### CALC_MOTIF_FREQUENCY
#     Currently hardwired for 3rd-order model
#     Calculates adjusted frequency (q) of nucleotide at last position of Markov chain,
#        expressed as log(q):
#
#           log_q = log_e ( model{a}{b}{c}{d} / sum )
#
#     where:
#        model{a}{b}{c}{d} is the number of times the nucleotide abcd appeared in the training set
#        sum is the number of times the nucleotides abc appeared in the training set
#   
  sub calc_motif_frequency {
     my ($nuc1, $nuc2, $nuc3, $nuc4);
     my $sum;
     foreach $nuc1 ("", "A", "C", "G", "T") {
        foreach $nuc2 ("A", "C", "G", "T") {
           foreach $nuc3 ("A", "C", "G", "T") {
              if (defined($model{$nuc1}{$nuc2}{$nuc3}{$sum_symbol})) { 
                 $sum = $model{$nuc1}{$nuc2}{$nuc3}{$sum_symbol};
                 foreach $nuc4 ("A", "C", "G", "T") {
                    if (defined($model{$nuc1}{$nuc2}{$nuc3}{$nuc4})) {   
                      $log_q{$nuc1}{$nuc2}{$nuc3}{$nuc4} =
                           log( $model{$nuc1}{$nuc2}{$nuc3}{$nuc4} / $sum );
                    }
                 }
              }
           }
        }
     }
     ##### Beginning sequences "...T", "..TT", etc are made neutral so as not to
     ##### penalize rare start codons more than once
     foreach $nuc3 ("", "A", "C", "G", "T") {
        foreach $nuc4 ("A", "C", "G", "T") {
           $log_q{""}{""}{$nuc3}{$nuc4} = 0;
        }
     }
  }
  
##### CALC_MARKOV_SCORE (sequence)
#     Hard-wired for 3rd order model
#     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
#     Looks at frequency of nucleotide just beyond the window
#     Begins window at boundary, so if sequence begins ATGCCG, window begins
#        at " AT"
  sub Calc_Markov_score {
     my ($seq) = @_;
     my ($nuc1, $nuc2, $nuc3) = (0, 1, 2);          # To avoid confusion by Perl coordinates
     my $nuc4;                                      # order+1th nucleotide
     my @nucs = ( "", split("", substr($seq,0,$order-1)) );  # Positions window at beginning
     my $score = 0;
     foreach my $position ($order-1..length($seq)-1) {       # Moves order+1th nucleotide through sequence
        $nuc4 = substr($seq, $position, 1);
        if (defined($log_q{$nucs[$nuc1]}{$nucs[$nuc2]}{$nucs[$nuc3]}{$nuc4})) {
           $score += $log_q{$nucs[$nuc1]}{$nucs[$nuc2]}{$nucs[$nuc3]}{$nuc4} - $log_p{$nuc4};
        } else {
           print "WARNING: |", join("|", $nucs[$nuc1], $nucs[$nuc2], $nucs[$nuc3], $nuc4);
           print "| not found in model!", $LF;
           print "Press any key to continue";
           <STDIN>;
           $score += log( $very_small_number ) - $log_p{$nuc4};
        }
        push @nucs, $nuc4;                           # Add order+1th nucleotide to window
        shift @nucs;			             # Move window by one
     }
     return $score;
  }
