#!/usr/bin/perl -w
# Vote
#    VERSION: Version 1 (11 November 2003)
#    PURPOSE: Uses set of best predictive genes determined from training set
#             Goes through independent set of Golub et al
#             Predicts ALL vs AML distinction
#             Prediction based on weighted voting scheme of Golub et al
#    INPUT FILES:  Training set of Golub et al
#             FORMAT: (tab-delimited)
#                 Header line
#                 Name_of_gene  Accession_number  expression_value [P|A|M]
#                                                 |------ 38 times ------|
#                    1 - 27 = ALL data for patients
#                   28 - 38 = AML data for patients
#                   (each line is data for one gene)
#    OUTPUT FILES: Description of format of output file
#             FORMAT: (tab-delimited)
#                 Header line
#                 Line_number  Name_of_gene  expression_value (38 times)  correlation_value
#    INCOMPLETE! SUBROUTINE CALCULATE_VOTE NEEDS TO BE COMPLETED!

############## LIBRARIES AND PRAGMAS ################
  use strict;

#################### CONSTANTS ######################

  my $LF = "\n";
  my $tab = "\t";

  my $number_of_ALL_patients = 27;
  my $number_of_AML_patients = 11;
  my $number_of_patients;

# Columns in training set devoted to expression data for each patient
  my @items_of_ALL_patients = Multiply(2, (1 .. 27) );
  my @items_of_AML_patients = Multiply(2, (28 .. 38) );
  my @items_of_patients;

  my $PS_threshold = 0.3;

#################### VARIABLES ######################

  my @indept_set;           # Holds all data for all patients and genes
  my @patient_names;        # First line of training set. Will be modified for output
  my $line;                   # Current line of training set under consideration
  my $name;                   # Name of current gene in training set
  my $accession;              # Name of current accession
  my $correlation;
  my @accessions;
  my @ALL_values;             # Expression values of current gene, ALL patients only
  my @AML_values;             # Expression values of current gene, AML patients only
  my %mean_ALL;               # For each predictor gene, mean of ALL expression 
                              #    keyed by accession: $mean_ALL{$accession}
  my %mean_AML;               # For current predictor gene, mean of AML expression 
                              #    keyed by accession: $mean_AML{$accession}
  my %avg_mean;               # For each predictor gene, (mean_ALL + mean_AML)/2
                              #    keyed by accession: $avg_mean{$accession}
  my %correlations;           # For each predictor gene, correlation value
                              #    keyed by accession: $mean_ALL{$accession}
  my $mean_ALL;               # Mean of putative ALL values for current gene
  my $mean_AML;               # Mean of putative AML values for current gene
  my $number_of_genes = 0;    # Count of predictor genes 
  my $patient;                # Current patient of independent set
  my @items;                  # Current list of expression values
  my $decision;               # ALL or AML or ?
  my $prediction_strength;    # (Vmin - Vloss)/(Vwin + Vloss) per Golub et al

###################### FILES ########################

  my $independent_file = 'data_set_ALL_AML_independent.txt';
  open INDEPT_SET, "<$independent_file" or die "Can't open $independent_file: $!\n";

  my $best_genes_file = 'data_set_ALL_AML_best.txt';

################### MAIN PROGRAM ####################
#    STRATEGY
#       1. Read set of predictor genes and calculate needed quantities
#       2. Read independent set
#          3. Extract expression data from one patient
#          4. Calculate vote, based on comparison with predictor gene expression pattern
#          5. Repeat steps 3-4 for each patient of independent set

  print $LF, "Reading set of best genes...";
  Calculate_means_from_best_genes ( );

  print $LF, "Reading independent set...";
  Read_independent_set ( );
  print " ", Size(@indept_set), " genes read.";
  
  print $LF, "Calculating correlations...";

  foreach $patient (@items_of_patients) {
     @items = Get_data_from_patient ($patient);
     ($decision, $prediction_strength) = Calculate_vote(@items);
     print join($tab, $patient_names[$patient], $decision, $prediction_strength), $LF;
  }     

#################### SUBROUTINES ####################

##### CALCULATE_MEANS_FROM_BEST_GENES
#     Opens and reads list of predictor genes
#     Calculates and saves %mean_ALL, %mean_AML, %avg_mean
#     Saves correlations for each gene in %correlations
#     All keyed to accession number
  sub Calculate_means_from_best_genes {
     my $line;                        # Temporary storage for one line of data
     my $gene;
     my @items;                       # Temporary storage for data from one line
     my ($mean_ALL, $mean_AML);
     
     open BEST_GENES, "<$best_genes_file" or die "Can't open $best_genes_file: $!\n";

     $line = <BEST_GENES>;          # Read second line of file
     while (defined $line) {
        $number_of_genes = $number_of_genes + 1;
        chomp $line;                  # Remove linefeed from end of line
        @items = split($tab, $line);  # Extract data from line
        Get_data_from_line(@items);
        $mean_ALL = Calculate_mean(@ALL_values);
        $mean_AML = Calculate_mean(@AML_values);
        $mean_ALL{$accession} = $mean_ALL;
        $mean_AML{$accession} = $mean_AML;
        $avg_mean{$accession} = ($mean_ALL + $mean_AML) / 2 ;
        $correlations{$accession} = $correlation;
        $line = <BEST_GENES>;       # Read next line
     }
  }

##### READ_INDEPENDENT_SET 
#     Reads training set file
#     Returns two dimensional array consisting of all relevant items in all lines
#        @indept_set[line][item]
#     Saves expression values only for those genes in predictor set
#     Saves accession numbers also
  sub Read_independent_set {
     my $line;                        # Temporary storage for one line of data
     my @items;                       # Temporary storage for data from one line

     open INDEPT_SET, "<$independent_file" or die "Can't open $independent_file: $!\n";
     my $header = <INDEPT_SET>;                # Read header line
     chomp $header;
     @patient_names = split($tab, $header);    # Save patient names from header
     shift @patient_names;                     # Ignore 1st column of header
     shift @patient_names;                     # Ignore 2nd column of header

     $line = <INDEPT_SET>;                     # Read second line of file
     while (defined $line) {                   # Do only while data left in file
        chomp $line;                           # Remove linefeed from end of line
        @items = split($tab, $line);           # Extract data from line
        shift @items;
        $accession = shift @items;
        if (defined $avg_mean{$accession}) {   # if gene part of predictor set...
           push @accessions, $accession;
           push @indept_set, [@items]; # Store data in array
        }
        $line = <INDEPT_SET>;       # Read next line
     }
     $number_of_patients = Size(@items) / 2;
     @items_of_patients = Multiply(2, (0 .. $number_of_patients-1) );
     close INDEPT_SET;
  }

##### GET_DATA_FROM_LINE (list of expression values)
#     Processes a given set of expression values
#     Extracts information on gene expression for ALL patients and AML patients
#         and all patients (the two subsets combined)
#     Extracts: name       The name of the gene
#               accession  The accession number of the gene
#               values     The gene expression values
#                          @ALL_values, @AML_values
  sub Get_data_from_line {
     my @items = @_;
     my $item;

     $name = $items[0];
     $accession = $items[1];
     $correlation = $items[-1];
     (@ALL_values, @AML_values) = ( );

     ### Extract info from ALL patients
     foreach $item (@items_of_ALL_patients) {
        push @ALL_values, $items[$item];
     }

     ### Extract info from AML patients
     foreach $item (@items_of_AML_patients) {
        push @AML_values, $items[$item];
     }
 
  }

##### GET_DATA_FROM_PATIENT (patient)
#     Gets expression data for given patient from indpendent set
  sub Get_data_from_patient {
     my ($patient) = @_;
     my @items = ();
     my $gene;
   
     foreach $gene (0 .. $number_of_genes - 1) {
        push @items, $indept_set[$gene][$patient];
     }
     return @items;
  }

##### CALCULATE_VOTE (@items)
#     Calculates vote per Golub et al from expression values of given gene for patient
  sub Calculate_vote {
     my @items = @_;        # Expression values for current patient

     ### GUTS OF SUBROUTINE GOES HERE!

  }

##### CALCULATE_MEAN (values)
#
  sub Calculate_mean {
     my @values = @_;
     my $value;
     my $sum = 0;
     my $n = 0;

     foreach $value (@values) {
        $sum = $sum + $value;
        $n = $n + 1;
     }

     if ($n == 0) {
        print "No values given in sub Calculate_mean";
        die;
     } else {     
        return $sum/$n
     }
  }

#### MULTIPLY (scalar, array)
#    Multiplies each element of the array by the scalar value
#    Returns array holding product of multiplication
  sub Multiply {
    my ($scalar,@array) = @_;
    my @product;
    foreach my $i (0 .. @array - 1) {
       $product[$i] = $scalar * $array[$i];
    }
    return @product;
  } 

##### SIZE (array)
#     Returns size of array 
#     Synonym for scalar (array)
  sub Size {return scalar(@_)}

##### MIN (list)
#     Returns minimum value in a given list
  sub Min {
     my @list = @_;
     my $min;

     foreach my $x (@list) {
        if (not defined $min or $x < $min) {$min = $x}
     }
     return $min;
  }

##### MAX (list)
#     Returns maximum value in a given list
  sub Max {
     my @list = @_;
     my $max;

     foreach my $x (@list) {
        if (not defined $max or $x > $max) {$max = $x}
     }
     return $max;
  }

