#!/usr/bin/perl -w
# Class_predictor-complete
#    VERSION: Version 3 (9 November 2003)
#    PURPOSE: Goes through data set of Golub et al
#             Finds genes that best predict ALL/AML distinction
#             Saves to file correlation value, name, and expression data of each gene
#             Prints to screen correlation value and name of each gene
#    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 AND PRAGMAS ################
  use strict;

#################### 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 @items_of_ALL_patients = Multiply(2, (1 .. 27) );
  my @items_of_AML_patients = Multiply(2, (28 .. 38) );

  my $First_informative_gene = 59;  # Previous genes are controls

#################### VARIABLES ######################

  my @training_set;           # Holds all data for all patients and genes
  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

###################### FILES ########################

  my $training_file = 'E:\IntroBioinf2002\Microarray\data_set_ALL_AML_train.txt';
# my $output_file = 'PUT SOMETHING HERE';

################### MAIN PROGRAM ####################

  print "Reading training set...";
  Read_training_set ( );

  print $LF, "Calculating correlations...";

  ### SOMETHING USEFUL GOES HERE
  ### You'll no doubt make good use of the subroutine Get_data_from_line
  ###    which extracts for you expression values and other things (take a look)
  ### You'll also want to use the subroutine Save_correlation_info
  ###    after every calculation of a correlation measure

  print $LF, "Sorting correlations...";

  ### SOMETHING USEFUL GOES HERE
  ### Note that there is already in this program a subroutine called
  ###    Sort_training_set_by_correlation, which could well be pertinent.

  Write_best_genes ( );

#################### SUBROUTINES ####################

##### READ_TRAINING_SET 
#     Reads training set file
#     Produces two dimensional array consisting of all items in all lines:
#        @training_set[line][item]
  sub Read_training_set {
     my $line;                        # Temporary storage for one line of data
     my @items;                       # Temporary storage for data from one line

     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
        @items = split($tab, $line);  # Extract data from line
        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: name       The name of the gene
#               accession  The accession number of the gene
#               values     The gene expression values
#                          @ALL_values, @AML_values, @values
#               presence   The P|A|M code
#                          @ALL_presence, @AML_presence, @presence
#     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;

     $name = $items[0];
     $accession = $items[1];
     (@values, @ALL_values, @AML_values) = ( );
     (@presence, @ALL_presence, @AML_presence) = ( );

     ### Extract info from ALL patients
     foreach $item (@items_of_ALL_patients) {
        push @ALL_values, $items[$item];
        push @ALL_presence, $items[$item+1];
     }

     ### Extract info from AML patients
     foreach $item (@items_of_AML_patients) {
        push @AML_values, $items[$item];
        push @AML_presence, $items[$item+1];
     }
 
     ###  Combine patient information into arrays with all info   
     push @values, @ALL_values, @AML_values;
     push @presence, @ALL_presence, @AML_presence;
  }


##### SAVE_CORRELATION_INFO (line, correlation)
#     Saves the correlation value at beginning of line of @training_set 
#         corresponding to a gene 
  sub Save_correlation_info {
     my ($gene_number, $correlation) = @_;
     unshift @{$training_set[$gene_number]}, $correlation;
  }

##### SORT_TRAINING_SET_BY_CORRELATION ("Ascending" | "Descending")
#     Sorts @training_set by the correlation (position 0 of the array)
#     Sorts ascending (low to high) or descending (high to low) depending
#         on given argument to subroutine
#     Actually, only the first letter of the argument matters
  sub Sort_training_set_by_correlation {
     my ($direction) = @_;
     if (uc(substr($direction, 0, 1)) eq "A") {
        @training_set = sort { $$a[0] <=> $$b[0] } @training_set;
     } else {
        @training_set = sort { $$b[0] <=> $$a[0] } @training_set;
     }
  }

##### WRITE_BEST_GENES
#     This subroutine doesn't do anything useful at the moment
#     When it does, it will include a call to Print_gene_info ( )
#     You'll also need to define $output_file in the FILE section
#         and remove the # in the open statement below
  sub Write_best_genes {
#    open OUTPUT, ">$output_file" or die "Can't open $output_file: $!\n";
#    print OUTPUT "correlation", $tab, $header, $LF;
     print $LF, "correlation", $tab, $header, $LF;
     print $LF, $LF, "This routine should print the list of best genes to a file " .
               "and perhaps to the screen", $LF;
  }

##### PRINT_GENE_INFO (line number)
#     Prints one gene's worth of information, including the correlation measure
#     Take out the pause and the # in front of the last line when ready to go
  sub Print_gene_info {
     my ($line_number) = @_;
     my ($correlation, $name, @items) = @{$training_set[$line_number]};
     print join($tab, $correlation, $name), $LF;
     <STDIN>               # TESTING PURPOSES ONLY!!!
#    print OUTPUT join($tab, $line_number, $name, @items, $correlation), $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(@_)}

