#!/usr/bin/perl
#
# Written by Jean-Philippe Vert
# July 3rd, 2001
#
# This script is used in the paper:
# J.-P. Vert, "Support vector machine prediction of signal peptide cleavage site using a new class of kernels for strings", proceedings of the Pacific Symposium on Biocomputing 2002.
#
# Its goal is to test a method to predict the signal peptide cleavage site of proteins using a support vector machine based on a new kernel. More details can be found in the original paper.
#
# Contact information: Jean-Philippe.Vert@mines.org
# http://web.kuicr.kyoto-u.ac.jp



######################################################################
# Change the following parameters to reflect your files organization #
######################################################################

# The file which contains the original signal protein data (obtained by concatenating files from http://www.cbs.dtu.dk/ftp/signalp/)
my $protfilename = "SIG.red";

# A path to write temporary files
my $tmpfilepath = "data/";

# The 'mySVM' executable for training the SVM (called 'mysvm' in the mySVM(v2.0) distribution)
my $mysvm = "bin/Linux-2.2.18-0vl4.2smp/mysvm";

# The 'mySVM' executable for prediction (called 'predict' in the mySVM(v2.0) distribution):
my $predict = "bin/Linux-2.2.18-0vl4.2smp/predict";

# The parameter file for mySVM
my $paramfile = "param.data";

############################################################################
# You should not have to change anything below this line to run the script #
############################################################################


# Reading the command line parameters #

if (!($#ARGV == 3)) {
    die "Syntax : ./kernel.pl n1 n2 iteration proportion\nn1 is the number of amino-acids before the cleavage site\nn2 is the number of amino-acid after the cleavage site\nhence the window has a size n1+n2\niteration is the number of iteration (different training/test sets)\nproportion (between 0 and 1) is the proportion of elements in the training set\n";
}


my $n1 = $ARGV[0];            # number of amino-acids before the cleavage site
my $n2 = $ARGV[1];            # number of amino-acids after the cleavage site
my $totaliter = $ARGV[2];     # number of iterations for the main loop
my $proportion = $ARGV[3];    # proportion of elements in the training set vs test set
my $windowsize = $n1 + $n2;   # window size


# Result files containing positive / negative examples, and training / test sets for mySVM
my $outfilepos="positive.exa";
my $outfileneg="negative.exa";
my $mysvmfiletrain="svmtrain.data";
my $mysvmfiletest="svmtest.data";


# Global variables

my %aminoacids = ("A" => 0,
		  "R" => 1,
		  "N" => 2,
		  "D" => 3,
		  "C" => 4,
		  "Q" => 5,
		  "E" => 6,
		  "G" => 7,
		  "H" => 8,
		  "I" => 9,
		  "L" => 10,
		  "K" => 11,
		  "M" => 12,
		  "F" => 13,
		  "P" => 14,
		  "S" => 15,
		  "T" => 16,
		  "W" => 17,
		  "Y" => 18,
		  "V" => 19);
my $aminonb = 20;

my @result = ();      # contains the results for the weight matrix method
my @result2 = ();     # contains the results for the SVM method








####################################
# Start of the main iteration loop #
####################################

my $ite;
for ($ite = 0 ; $ite < $totaliter ; $ite++) {

  # Each iteration of the main loop consists of the following:
  # 1 - construction of a training set and a test set (by random splitting)
  # 2 - estimation of the weight matrix on the training set
  # 3 - test of the weight matrix on the test set
  # 4 - Creation of training and test set files for mySVM
  # 5 - SVM training on the training set
  # 6 - SVM test on the test set
  
  print "Iteration $ite\n";
  
  # Variable initialization to count amino-acids at different positions
  my @count1;   # count the amino-acids in the positive training set
  my @count2;   # count the amino-acids in the negative training set
  my @weight;   # weight matrix
  my $i;
  for ($i=0 ; $i<$windowsize ; $i++) {
    @count1 = (@count1 , [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]);
  }
  for ($i=0 ; $i<$windowsize ; $i++) {
    @count2 = (@count2 , [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]);
  }
  for ($i=0 ; $i<$windowsize ; $i++) {
    @weight = (@weight , [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]);
  }
  
  
  
  ###########################################################################
  # 1 - Construction of a training set and a test set (by random splitting) #
  # Each file contains a list aminoacid sequences, each of size $windowsize #
  ###########################################################################
  
  # Open files to write the different sets:
  # POSFILE1 : positive training set
  # POSFILE2 : positive test set
  # NEGFILE1 : negative training set
  # NEGFILE2 : negative test set
  open (POSFILE1 , ">".$tmpfilepath.$outfilepos."1") || die "Cannot open $tmpfilepath$outfilepos 1 : $!";
  open (POSFILE2 , ">".$tmpfilepath.$outfilepos."2") || die "Cannot open $tmpfilepath$outfilepos 2 : $!";
  open (NEGFILE1 , ">".$tmpfilepath.$outfileneg."1") || die "Cannot open $tmpfilepath$outfileneg 1 : $!";
  open (NEGFILE2 , ">".$tmpfilepath.$outfileneg."2") || die "Cannot open $tmpfilepath$outfileneg 2 : $!";
  
  # Open the protein file
  open (INFILE , $protfilename) || die "Cannot open $protfilename\n";
    
  my $line;
  while ($line = <INFILE>) {

    # The protein are read one by one
    # '$line1' contains the description of the protein

    my $line2;
    my $line3;
    my $length;
    
    # Read the amino-acid sequence (i.e., next line in the file)
    ($line2 = <INFILE>) || die "Problem with the file $protfilename: I can't read the amino-acid sequence of $line1\n";	
    
    # Read the corresponding signal indicator (i.e., next line in the file)
    ($line3 = <INFILE>) || die "Problem with the file  $protfilename : I can't read the signal indicator sequence of $line1\n";
    
    # Length of the sequence
    $length = length $line2;
    
    # For each sequence, all possible windows of size '$windowsize' are extracted
    for ($i = 0 ; $i <= ($length - $windowsize) ; $i++) {
	
	my $r = rand;    # '$r' determines whether the window is in the training or in the test set

	my $toprint = substr($line2 , $i , $windowsize);  # '$toprint' corresponds to one window
	
	if (substr($line3 , $i + $n1 , 1) eq 'C') {
	  # In that case the window is a cleavage site
	  if ($r < $proportion) {
	    print POSFILE1 "$toprint\n";
	  }
	  else {
	    print POSFILE2 "$toprint\n";
	  }

	}
	else {
	  # in that case it is a negative example
	  if ($r < $proportion) {
	    print NEGFILE1 "$toprint\n";
	  }
	  else {
	    print NEGFILE2 "$toprint\n";
	  }
	}
	
	
	if ($r < $proportion) {
	  # If it is a training example then update the corresponding counters
	  my $j;
	  for ($j=0 ; $j < $windowsize ; $j++) {
	    my $c = substr($toprint , $j , 1);
	    if (substr($line3 , $i + $n1 , 1) eq 'C') {
	      $count1[$j][$aminoacids{$c}]++;   # count1 is the counter for the positive examples
	    }
	    else {
	      $count2[$j][$aminoacids{$c}]++;   # count2 is the counter for the negative examples
	    }
	  }
	}
      }
  }

  close (INFILE);
  close (POSFILE1);
  close (POSFILE2);
  close (NEGFILE1);
  close (NEGFILE2);

  print "Training and test sets created\n";


  ###########################################################
  # 2 - estimation of the weight matrix on the training set #
  ###########################################################

  # The weight of an amino-acid at a particuliar position is obtained by dividing its frequency in the positive training set by its "natural" frequency (in the positive AND negative training set), and then taking the logarithm of this quotient.

  my $sum1 = 0;
  my $sum2 = 0;
  my $j;

  # Compute the total number of examples in the positive training set ('$sum1') and in the total training set ('$sum2')
  for ($j = 0 ; $j < $aminonb ; $j++) {
    $sum1 += $count1[0][$j];
    $sum2 += $count2[0][$j];
  }

  print "$sum1 positive examples in the training set\n";
  print "$sum2 negative examples in the training set\n";

  # Compute the weight matrix
  for ($i=0 ; $i<$windowsize ; $i++) {
    for ($j = 0 ; $j < $aminonb ; $j++) {
      # '@weight' is the weight matrix
      $weight[$i][$j] = log( $count1[$i][$j] + 1 ) - log( $sum1 + $aminonb ) - log(  $count1[$i][$j] + $count2[$i][$j] + 1 ) + log( $sum1 + $sum2 + $aminonb );

      # '@count1' is set to contain the logarithm of the estimated probability on the positive training set:
      $count1[$i][$j] = log( $count1[$i][$j] + 1 ) - log( $sum1 + $aminonb );
    }
  }
  print "Weight matrix computed\n";

  #    for ($i=0 ; $i<$windowsize ; $i++) {
  #	print "@{$weight[$i]}\n";
  #    }


  #################################################
  # 3 - test of the weight matrix on the test set #
  #################################################

  my %scores =();   # a hash variable to contain the score of each protein in the test set

  # First on the positive test set
  open (POSFILE2 , $tmpfilepath.$outfilepos."2") || die "Cannot open $tmpfilepath$outfilepos 2 : $!";
  my $number = 0;
  while ($line = <POSFILE2>) {
    # For each example in the test set:
    ++$number;
    my $score = 0;
    for ($i = 0 ; $i < $windowsize ; $i++) {      #scan over the window
      my $c = substr($line , $i , 1);
      $score += $weight[$i][$aminoacids{$c}];     # add the weight of the observed amino-acid
    }
    # store the score in a hash variable
    my $name = "p".$number;
    $scores{$name} = $score;
  }
  close (POSFILE2);

  # Then on the negative test set
  open (NEGFILE2 , $tmpfilepath.$outfileneg."2") || die "Cannot open $tmpfilepath$outfilepos 2 : $!";
  $number = 0;
  while ($line = <NEGFILE2>) {
    # For each example in the test set:
    ++$number;
    my $score = 0;
    for ($i = 0 ; $i < $windowsize ; $i++) {      #scan over the window
      my $c = substr($line , $i , 1);
      $score += $weight[$i][$aminoacids{$c}];     # add the weight of the observed amino-acid
    }
    my $name = "n".$number;
    $scores{$name} = $score;
    }
  close (NEGFILE2);

  # sort the test elements by decreasing score
  my @sorted_keys = sort {
    $scores{$b} <=> $scores{$a} || $a cmp $b
  } keys (%scores);

  my $key;
  $i=0;
  $j=0;

  # Compute the true positive vector
  my @localresul = ();

  foreach $key (@sorted_keys) {
    #	print "$key :\t$scores{$key}\n";
    if (substr($key,0,1) eq "p") {
      $i++;
    }
    else {
      $localresult[$j++] = $i;
    }
  }

  # Normalize it and store it in the global variable $result
  my $k;
  for ($k=0 ; $k < 100 ; $k++) {
    $result[$k] += $localresult[int ($k * $j / 100)] * 100 / $i;
  }
  #    print "@{$result[$ite]}\n";
  print "Weighted matrix method done\n";


  ##########################################################
  #  4 - Creation of training and test set files for mySVM #
  ##########################################################

  # Creation of the training file
  open (SVMFILE , ">".$tmpfilepath.$mysvmfiletrain) || die "Cannot open $tmpfilepath$mysvmfiletrain: $!";

  # The first line contains information about the format of the data ('@examples' and 'format yx')
  print SVMFILE "\@examples\nformat yx\n";

  # Process the training file of positive examples
  open (INFILE , $tmpfilepath.$outfilepos."1") || die "Cannot open $tmpfilepath$outfilepos 1 : $!";
  while ($line = <INFILE>) {
    # this is a positive example so print '+1'
    print SVMFILE "+1 ";
    for ($j=0 ; $j < $windowsize ; $j++) {
      my $c = substr($line , $j , 1);
      # Print the amino-acid index followed by its log-probability
      printf SVMFILE "%d %.4f ",$aminoacids{$c}+1,$count1[$j][$aminoacids{$c}];
    }
    print SVMFILE "\n";
  }
  close (INFILE);

  # Process the training negative examples
  open (INFILE , $tmpfilepath.$outfileneg."1") || die "Cannot open $tmpfilepath$outfilepos 1 : $!";
  while ($line = <INFILE>) {
    print SVMFILE "-1 ";
    for ($j=0 ; $j < $windowsize ; $j++) {
      my $c = substr($line , $j , 1);
      printf SVMFILE "%d %.4f ",$aminoacids{$c}+1,$count1[$j][$aminoacids{$c}];
    }
    print SVMFILE "\n";
  }
  close (INFILE);
  close (SVMFILE);    # training file done

  # Creation of the test file
  open (SVMFILE , ">".$tmpfilepath.$mysvmfiletest) || die "Cannot open $tmpfilepath$mysvmfiletest: $!";
  print SVMFILE "\@examples\nformat yx\n";

  # Test positive examples
  open (INFILE , $tmpfilepath.$outfilepos."2") || die "Cannot open $tmpfilepath$outfilepos 2 : $!";
  while ($line = <INFILE>) {
    print SVMFILE "+1 ";
    for ($j=0 ; $j < $windowsize ; $j++) {
      my $c = substr($line , $j , 1);
      printf SVMFILE "%d %.4f ",$aminoacids{$c}+1,$count1[$j][$aminoacids{$c}];
    }
    print SVMFILE "\n";
  }
  close (INFILE);

  # Training negative examples
  open (INFILE , $tmpfilepath.$outfileneg."2") || die "Cannot open $tmpfilepath$outfilepos 1 : $!";
  while ($line = <INFILE>) {
    print SVMFILE "-1 ";
    for ($j=0 ; $j < $windowsize ; $j++) {
      my $c = substr($line , $j , 1);
      printf SVMFILE "%d %.4f ",$aminoacids{$c}+1,$count1[$j][$aminoacids{$c}];
    }
    print SVMFILE "\n";
  }
  close (INFILE);
  close (SVMFILE);    # test file done

  print "Files created for mySVM\n";
  
  ########################################
  # 5 - SVM training on the training set #
  ########################################
  
  print "SVM learning...\n";

  system $mysvm , $paramfile , $tmpfilepath.$mysvmfiletrain ;




  ################################
  # 6 - SVM test on the test set #
  ################################
    
  print "SVM testing...\n";


  system $predict , $paramfile , $tmpfilepath.$mysvmfiletrain.".svm" , $tmpfilepath.$mysvmfiletest ; 




  # Fill the vector of results for the SVM method
  %scores =();
  my @names;

  open (INFILE , $tmpfilepath.$mysvmfiletest) || die "Cannot open file  $tmpfilepath$svmtestfile\n";
  my $pos=0;
  my $neg=0;
    
  # Create a list of names for each example
  while ($line = <INFILE>) {
    if (substr($line,0,1) eq "+") {
      # positive example;
      $names[($pos++)+($neg)] = "p".$pos;
    }
    if (substr($line,0,1) eq "-") {
      # negative example;
      $names[($pos)+($neg++)] = "n".$neg;
    }
  }
  close (INFILE);
    
  # Read the corresponding scores and store them in a hash variable
    
  open (INFILE , $tmpfilepath.$mysvmfiletest.".pred") || die "Cannot open file  $tmpfilepath$svmtestfile .pred\n";  
  $k=0;
  while ($line = <INFILE>) {
    if (!(substr($line,0,1) eq "#")) {
      chomp $line;
      $scores{$names[$k++]} = $line; 
    }
  }

  if (!( $pos + $neg == $k)) {
    die "Error : i read $pos positive and $neg negative examples in ".$tmpfilepath.$mysvmfiletest."\nBut $k scores in ".$tmpfilepath.$mysvmfiletest.".pred\n";
  }

  # sort the test elements by decreasing score
  my @sorted_keys = sort {
    $scores{$b} <=> $scores{$a} || $a cmp $b
  } keys (%scores);

  my $key;
  $pos = 0;
  $neg = 0;
    
  # Compute the true positive vector
  my @localresul = ();
    
  foreach $key (@sorted_keys) {
    #	print "$key :\t$scores{$key}\n";
    if (substr($key,0,1) eq "p") {
      $pos++;
    }
    else {
      $localresult[$neg++] = $pos;
    }
  }

  for ($k=0 ; $k < 100 ; $k++) {
    $result2[$k] += $localresult[int ($k * $neg / 100)] * 100 / $pos;
  }

    
    
  ##########################
  # Print out the results
  ########################
  print "Average performance after $ite iterations:\n";
  print "Average result for the weight matrix method 1:\n";
  my $k;
  for ($k=0 ; $k < 100 ; $k++) {
    printf "%.1f ", $result[$k] / ($ite + 1);
  }
  print "\nAverage result for the SVM method :\n";
  for ($k=0 ; $k < 100 ; $k++) {
    printf "%.1f ", $result2[$k] / ($ite + 1);
  }
  print "\n";
    


  
  ##################################
  # End of the main iteration loop #
  ##################################
}


