#!/usr/bin/perl -w
# Permute_training_set
#    VERSION: Version 1 (10 November 2003)
#    PURPOSE: Permutes patient labels in data set of Golub et al
#             Calculates measure of correlation for each gene, each permutation
#             For each number of genes, determines n% best measure of correlation
#    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

######## LIBRARIES, SUBROUTINES, AND PRAGMAS ########
  use strict;
  sub Permute (\@\@);
 
#################### CONSTANTS ######################

  my $LF = "\n";
  my $tab = "\t";

  my $number_of_patients = 38;
  my $number_of_ALL_patients = 27;
  my $number_of_AML_patients = 11;

  ### Columns in training set devoted to expression data for each patient
  my @columns_with_expression_data = Multiply(2, (1 .. 38) );
  my @items_of_ALL_patients = (0 .. 26);
  my @items_of_AML_patients = (27 .. 37);

  my $first_informative_gene = 59;  # Previous genes are controls

  my $significance = .5;            # CHANGE for specific level of significance
  my $number_of_permutations = 20;  # Golub et al used 400

#################### VARIABLES ######################

  my @training_set;           # Holds all data for all patients and genes
                              #    2D array: $training_set[gene][patient]
                              #    Contains only expression data
  my $header;                 # 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 @values;                 # Expression values of current gene
  my @ALL_values;             # Expression values of current gene, ALL patients only
  my @AML_values;             # Expression values of current gene, AML patients only
  my @presence;               # [P|A|M] designation of current gene
  my @ALL_presence;           # [P|A|M] designation of current gene, ALL patients only
  my @AML_presence;           # [P|A|M] designation of current gene, AML patients only
  my @current_correlations;   # Measures of correlation for current permutation
                              #    Eventually gets sorted lowest--> highest
  my @correlations;           # Measures of correlation for all permutations
                              #    2D array: $correlations[ranking][permutation]
  my $permutation;            # Current permutation

  my $ALL_mean;
  my $AML_mean; 
  my $ALL_stdev;
  my $AML_stdev;
  my $correlation;

###################### FILES ########################

  my $training_file = 'data_set_ALL_AML_train.txt';
  my $output_file = 'ALL_AML_50pct.txt';

################### MAIN PROGRAM ####################
# STRATEGY:
#    1. Read training set
#       2. Permute patient labels
#       3. Calculate measure of correlation for each gene
#       4. Sort and store correlations for this permutation
#       5. Repeat steps 2-4 for as many permutations as specified
#    6. For each rank (i.e. best correlation, 2nd best, 3rd best...)
#           Sort the correlations, lowest --> highest
#    7. Output the correlation for each rank that is better than n% of the total,
#           where n is given significance

  print "Reading training set...";
  Read_training_set ( );
  print Size(@training_set), " genes processed.", $LF;

  print $LF, "Calculating correlations...";

  foreach $permutation (1 .. $number_of_permutations) {
     print $permutation, " ";
     Permute (@items_of_ALL_patients, @items_of_AML_patients);
     @current_correlations = ();
     foreach $line ($first_informative_gene - 1 .. Size(@training_set) - 1) {
        Get_data_from_line ($line);
        $ALL_mean = Calculate_mean(@ALL_values);
        $AML_mean = Calculate_mean(@AML_values);
        $ALL_stdev = Calculate_stdev($ALL_mean, @ALL_values);
        $AML_stdev = Calculate_stdev($AML_mean, @AML_values);
        $correlation = Calculate_correlation($ALL_mean, $ALL_stdev, $AML_mean, $AML_stdev); 
        Save_correlation_info($correlation);
     }

     @current_correlations = sort { $a <=> $b } @current_correlations; 
     Save_rankings ( );
  }

  Sort_rankings ( );
  Output_curve ( );

#################### SUBROUTINES ####################

##### READ_TRAINING_SET 
#     Reads training set file
#     Produces two dimensional array consisting of expression data in all lines:
#        @training_set[line][patient]
  sub Read_training_set {
     my $line;                        # Temporary storage for one line of data
     my @items;                       # Temporary storage for data from one line
     my @all_items;
     my $item;

     open TRAINING_SET, "<$training_file" or die "Can't open $training_file: $!\n";
     $header = <TRAINING_SET>;        # Read header line
     chomp $header;                   # Remove linefeed from end of line

     $line = <TRAINING_SET>;          # Read second line of file
     while (defined $line) {          # Do only while data left in file
        chomp $line;                  # Remove linefeed from end of line
        @all_items = split($tab, $line);  # Extract data from line
        @items = ( );
        foreach $item (@columns_with_expression_data) {
           push @items, $all_items[$item];
        }
        push @training_set, [@items]; # Store data in array
        $line = <TRAINING_SET>;       # Read next line
     }
  }

##### GET_DATA_FROM_LINE (line number)
#     Processes a line in the training set indicated by the given line number
#     Extracts information on gene expression for ALL patients and AML patients
#         and all patients (the two subsets combined)
#     Extracts: values     The gene expression values
#                          @ALL_values, @AML_values
#     Strategy: The specific array containing info for the gene is extracted
#                  from the training set, using the given line number
#               Since the training set is an array of references to each gene's
#                  array, it is necessary to do Perl magic 
#                         @{$training_set[$line_number]}
#                  to translate the reference into an actual array
  sub Get_data_from_line {
     my ($line_number) = @_;
     my @items = @{$training_set[$line_number]};
     my $item;

     (@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];
     }
 
  }

##### 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
     }
  }

##### CALCULATE_STDEV (mean, values)
#     Returns standard deviation of set of values
#     Could calculate mean itself but instead asks for it
  sub Calculate_stdev {
     my ($mean, @values) = @_;
     my $value;
     my $sum = 0;
     my $n = 0;

     foreach $value (@values) {
        $sum = $sum + ($value - $mean)**2;
        $n = $n + 1;
     }

     if ($n <2) {
        print "Need at least 2 values in sub Calculate_stdev";
        die;
     } else {     
        return sqrt($sum/($n-1));
     }
  }

##### CALCULATE_CORRELATION (first mean, first std dev, second mean, second std dev)
#     Measure of correlation per Golub et al (1999)
  sub Calculate_correlation {
     my ($mean1, $stdev1, $mean2, $stdev2) = @_;
     return (($mean1 - $mean2) / ($stdev1 + $stdev2));
  }

##### SAVE_CORRELATION_INFO (correlation)
#     Saves the correlation value in @current_correlations
  sub Save_correlation_info {
     my ($correlation) = @_;
     push @current_correlations, $correlation;
  }

##### SAVE_RANKINGS
#     Saves correlations into @correlations
#     Since @correlations is arranged with ranking as the first subscript
#        it has to be saved item by item
  sub Save_rankings {
     foreach my $ranking (0 .. Size(@current_correlations) - 1) {
        push @{$correlations[$ranking]}, $current_correlations[$ranking];
     }
  }

##### SORT_RANKINGS
#     Sorts correlations for each permutation by score (low --> high)
#     The set of correlations for a given ranking is stored in 
#        $correlations[$ranking]
  sub Sort_rankings {
     foreach my $ranking (0 .. Size(@current_correlations) - 1) {
        @{$correlations[$ranking]} = reverse sort @{$correlations[$ranking]}
     }
  }

##### OUTPUT_CURVE
#     Output one rank at a time, the measure of correlation then the ranking
#     This order makes it easier for Excel to make curve with X-axis as correlation
#     Output file is tab-delimited
  sub Output_curve {
     open OUTPUT, ">$output_file" or die "Can't open $output_file: $!\n";
     my $index_of_significance = ($number_of_permutations * $significance) - 1;
     foreach my $ranking (0 .. Size(@current_correlations) - 1) {
        print OUTPUT $correlations[$ranking][$index_of_significance], $tab, $ranking+1, $LF;
     }
  }

#### 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(@_)}

##### PERMUTE (array1, array2)
#     (Actually, prototype forces arguments into list of reference)
#     Combines two arrays and redistributes elements at random
#     Changes values of argument arrays, but not their sizes
  sub Permute (\@\@) {
     my ($array1_ref, $array2_ref) = @_;
     my @array1 = @$array1_ref;
     my @array2 = @$array2_ref;
     my @whats_left = (@array1, @array2);
     my @new_array1 = ();

     my $array1_size = @array1;
     my $item;

     foreach $item (1 .. $array1_size) {
        my $index = int( rand( Size(@whats_left) ) );
          ### random number between 0 and size of @whats_left
        push @new_array1, splice(@whats_left, $index, 1);
          ### Extracts the $index'th element of @whats_left
          ###   pushing it onto @new_array1
          ###   and knitting together the hole left behind in @whats_left
     }
     @$array1_ref = @new_array1;
     @$array2_ref = @whats_left;
       ### give constructed arrays back to arrays passed to subroutine
  }
