#!/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.


if ($#ARGV !=3) {
  die "This script creates two files to be used by SVMLight, built from the file containing phylogenetic profiles as a list of log(E-value) and a list of genes considered positive. There is one file for positive genes, one for negative. Each gene is transformed into a 0/1 phylogenetic profiles, and bits are permuted in order to use the phylogenetic kernel function in SVMLight.\n
Usage: prep.pl <phyloprofiles.txt> <positiveExamples.txt> <outpositive> <outnegative>\n";
}


my $orgnb = 24;   # number of organisms

my @permutation = (24,1,14,18,16,14,2,11,12,20,10,9,17,5,13,15,6,7,23,4,21,8,22,3);  # permutation of organism to have the correct order for the phylogenetic kernel


$phyloFile = $ARGV[0];         # File containing the original phylogenetic profiles
$posExampleFile = $ARGV[1];    # File containing a list of "positive genes" (a MIPS functional category)
$outpos = $ARGV[2];            # File to output positive genes
$outneg = $ARGV[3];            # File to output negative genes



# Read the list of "positive" genes in the hash 'pos'
open (INFILE , $posExampleFile)
  or die "Can't open '$posExampleFile': $!";
my $line;
my %pos;

while ($line = <INFILE>) {
  chomp $line;
  $line =~ /sce:(.*)$/;
  $pos{$1} = 1;
}
close (INFILE);


# Open the file containing the phylogenetic profiles as log(E-value), and the output files
open PHYLOF , $phyloFile
  or die "Can't open '$phyloFile': $!";
open FPOS , ">".$outpos
  or die "Can't open '$outpos': $!";
open FNEG , ">".$outneg
  or die "Can't open '$outneg': $!";

# Skip the first two lines
$line = <PHYLOF>;
$line = <PHYLOF>;

# Read each line (each gene) one by one
while ($line = <PHYLOF>) {
  chomp $line;

  # split the line
  my @l = split ' ',$line;

  # get the gene name (first column)
  $line = shift @l;

  if ($pos{$line}) {
    # If the gene is positive, print its 0/1 phylogenetic profile in the positive output file
    print FPOS  "+1 ";
    my $i;

    for ($i=1 ; $i<= $orgnb ; $i++){
      if ($l[$permutation[$i-1]-1] > 0) {
	# if -log(E-value)>0, print '+1'
	print FPOS "$i:1 ";
      }
      else {
	print FPOS "$i:0 ";
      }
    }
    print FPOS "\n";
  } else {
    # If the gene is negative, do the same on the negative output file
    print FNEG  "-1 ";
    my $i;

    for ($i=1 ; $i<= $orgnb ; $i++){
      if ($l[$permutation[$i-1]-1] > 0) {
	print FNEG "$i:1 ";
      }
      else {
	print FNEG "$i:0 ";
      }
    }
    print FNEG "\n";
  }
}


close PHYLOF;
close FPOS;
close FNEG;

