#!/usr/bin/perl
#
# Written by Jean-Philippe.Vert@mines.org , november 4, 2002
#
# Reference:
# Jean-Philippe Vert, "A tree kernel to analyze phylogenetic profiles", Bioinformatics, vol. 18, p. S276-S284, 2002.
#
# This script performs the experiment reported in the paper cited above. For all functional categories of the MIPS database, it computes the performance of two SVM (one with a linear kernel and one with a tree kernel) to predict the function of a gene from its phylogenetic profile.

##########################################################
# Parameters: change them accroding to your configuration
##########################################################

# File containing the phylogenetic profiles (as series of -log(E-value))
my $phyloFile = "../data/profiles.txt";

# Directory containing the MIPS classes (one file in text format for each functional category)
my $mipsPath = "../data/mips/txt";

# Directory to save the results
my $resultDir = "../result";

# Directory to write temporary files
my $tmpDir = "../tmp";

# Perl script to prepare files for SVMLight
my $perlPrep = "./prep.pl";

# Perl script to randomly split a file into two files
my $perlSplit = "./randsplit.pl";

# SVM programs
my $svmLearn = "../svmlight/svm_learn";
my $svmClassify = "../svmlight/svm_classify";

# Proportion of example used as training example during cross-validation
my $ratio = 0.66666;

# Number of iteration for each experiment (to average the ROC curves on a number of random train/test splits)
my $iter = 10;

# Maximum number of MIPS categories to process (useful to test the script on only a few classes)
my $nmax = 1000;



##########################################
# Program starts here
##########################################

# Variable initialization
my %resultLinear = ();   # hash table to store the ROC50 index for each class with a linear SVM
my %resultTree = ();     # hash table to store the ROC50 index for each class with a tree kernel SVM

my $proce = 0;            # number of categories processed

my $posTotal = $tmpDir."/posTotal.dat";
my $negTotal = $tmpDir."/negTotal.dat";
my $posTrain = $tmpDir."/posTrain.dat";
my $negTrain = $tmpDir."/negTrain.dat";
my $posTest = $tmpDir."/posTest.dat";
my $negTest = $tmpDir."/negTest.dat";
my $trainTotal = $tmpDir."/train.dat";
my $testTotal = $tmpDir."/test.dat";
my $svmModel = $tmpDir."/svm.model";
my $svmResult = $tmpDir."/svm.result";


# Store the list of MIPS category files
my @files = glob("$mipsPath/*.txt");

# Count the total number of genes
my $genesNb=0;
open (INFILE , $phyloFile)
  or die "Can't open $phyloFile: $!";
while (<INFILE>) {
  $genesNb++;
}
close INFILE;
print "$genesNb genes read in the file $phyloFile.\n";


# process MIPS categories one by one:
foreach (@files) {

  # process no more than $nmax categories:
  if ($proce++ < $nmax) {

    my $filename=$_;        # Name of the file containing the list of genes of the current category

    # Create the files of positive and negative examples
#    print "now prep : $perlPrep  $phyloFile $filename $posTotal $negTotal\n";
    system "$perlPrep  $phyloFile $filename $posTotal $negTotal";

    # Count the number of positive examples
    my $i=0;
    my $line;
    open INFILE , $posTotal;
    while ($line = <INFILE>) {
      $i++;
    }
    close INFILE;

    # Process the category if it contains at least 10 positive examples
    if ($i < 10) {
      print "Skip $filename ($i positive examples)\n";
    } else {
      print "Processing $filename ($i positive examples)\n";


      # Compute the factor to balance the number of positive and negative examples in SVMLight
      my $balance = ($genesNb-$i)/$i;

      # Local variable to store the performance on the current category
      my $roc50Linear = 0 ,$roc50Tree = 0 ;            # roc50 index
      my @roc = ();                                    # roc curve
      my @tp50Linear = () , @tp50Tree = ();            # number of true positive as a function of false positive (up to 50 false positives)


      # Perform $iter iteration of the classification experiment
      my $it;
      for ($it = 0 ; $it < $iter ; $it++) {
	
	# Randomly split positive and negative examples into training and test sets
	system "$perlSplit $posTotal $posTrain $posTest $ratio";
	system "$perlSplit $negTotal $negTrain $negTest $ratio";
	system "cat $posTrain $negTrain > $trainTotal";
	system "cat $posTest $negTest > $testTotal";
	
	# Test the linear SVM
	print "$svmLearn -j $balance -v 0 -m 300 $trainTotal $svmModel\n";
	system "$svmLearn -j $balance -v 0 -m 300 $trainTotal $svmModel";
	print "$svmClassify -v 0 $testTotal $svmModel $svmResult\n";
	system "$svmClassify -v 0 $testTotal $svmModel $svmResult";

	# Compute the performance of the linear SVM
	@roc = computeTp($testTotal , $svmResult);
	
#	print "@roc \n";

	for ($i=0 ; $i<=50 ; $i++) {
	  $tp50Linear[$i] += $roc[$i];
	}
#	$roc50Linear += computeRoc50(@tp50Linear);

#	print $roc50Linear."\n";

	

	# test the phylo SVM
	print "$svmLearn -j $balance -v 0 -t 4  -m 300 $trainTotal $svmModel";
	system "$svmLearn -j $balance  -t 4  -m 300 $trainTotal $svmModel";
	print "$svmClassify -v 0 $testTotal $svmModel $svmResult";
	system "$svmClassify  $testTotal $svmModel $svmResult";
	
	# Compute the performance of the tree kernel SVM
	@roc = computeTp($testTotal , $svmResult);
	for ($i=0 ; $i<=50 ; $i++) {
	  $tp50Tree[$i] += $roc[$i];
	}
#	$roc50Tree += computeRoc50(@tp50Tree);

      } # end of the iteration loop

      # Compute the ROC 50 index
      $roc50Linear = computeRoc50(@tp50Linear)/$iter;
      $roc50Tree = computeRoc50(@tp50Tree)/$iter;


      # Print out the results
      $filename =~ /\/([^\/]+)\.txt/;

#      print "File = $1 instead of $filename\n";
      my $resultFilename = $1;
      open (OUTFILE , ">$resultDir/$resultFilename")
	or die "Can't open $resultDir/$1: $!";

      print OUTFILE "Number of iterations : $iter\nRatio of examples in the training set : $ratio\n\nLinear kernel:\nROC50 = $roc50Linear\nROC50 curve : \n";
      for ($i=0 ; $i <= 50 ; $i++) {
	printf (OUTFILE  "%.1f ",$tp50Linear[$i] * 100 / $iter);
      }
      print OUTFILE "\n\nTree kernel:\nROC50 = $roc50Tree\nROC50 curve : \n";
      for ($i=0 ; $i <= 50 ; $i++) {
	printf (OUTFILE  "%.1f ",$tp50Tree[$i] * 100 / $iter);
      }
      close OUTFILE;

      # Store the ROC50 index in a hash table
      $resultLinear{$filename} = $roc50Linear;
      $resultTree{$filename} = $roc50Tree;

    }
  }
}

# Sort the functional categories by roc50 index obtained by the linear SVM
my @s = sort { $resultLinear{$b} <=> $resultLinear{$a} } keys %resultLinear;

open OUTFILE , ">$resultDir/roc50";
foreach $line (@s) {
  printf OUTFILE  "%.3f\t%.3f\t$line\n" , $resultLinear{$line} , $resultTree{$line};
}
close OUTFILE;




sub computeTp {
  # compute the vector of recall (true positives / total positive) as a function of false positives (for up to 50 false positives)
  # two arguments: test file used by SVM Light (with annotation) and predicted file

  # array to store the TP vector
  my @result=();

  # Extract from the files a hash with the class and score of each example
  open (TRUEFILE,$_[0])
    or die "Can't open $_[0]: $!";

  open (PREDFILE,$_[1])
    or die "Can't open $_[0]: $!";

  my %score = ();
  my $line, $p=0, $n=0;

  while ($line = <TRUEFILE>) {
    my $sc = <PREDFILE>;
    chomp $sc;
    if ($line =~ /\+/) {    #positive example
      $p++;
      $score{"p".$p} = $sc;
    } else {                # negative example
      $n++;
      $score{"n".$n} = $sc;
    }
  }
  close TRUEFILE;
  close PREDFILE;

  # sort examples by decreasing score
  my @sortedKeys =  sort { $score{$b} <=> $score{$a} } keys (%score);

  # Compute the true positive vector
  $p = 0;
  $n = 0;

  foreach $key (@sortedKeys) {
#    print $key."\t".$score{$key}."\n";
    if (substr($key,0,1) eq "p") {
      $p++;
    } else {
      if ($n <= 50) {
	$result[$n++] = $p;
      }
    }
  }

#  print "total : $p positive examples\n"; 

  # divide the number of true positives by the total number of positive to get the vector of recall
  for ($n=0 ; $n<=50 ; $n++) {
    $result[$n] /= $p;
  }

  # return the vector of recall as a function of false positives
  return @result;
}


sub computeRoc50 {
  # compute the roc50 index from the vector of recall as a function of FN

  my $r = 0;

  # Sum up the number of TP at 1,..,50 FP:
  my $i;
  for ($i=0 ; $i<50 ; $i++) {
    $r += $_[$i];
  }

  # Normalize to get the ROC50 index (normalized area under the ROC curve, between 0 and 1, 1 for a perfect classifier)
  $r /= 50;

  return $r;
}
