source: trunk/se/lu/thep/wenni/bin/base_plugin_script/wenni.pl @ 115

#!/usr/bin/perl -w

# $Id: wenni.pl 115 2006-06-22 08:43:07Z jari $

# Copyright (C) 2005, 2006 Jari Häkkinen
#
# This file is part of WeNNI,
# http://lev.thep.lu.se/trac/baseplugins/wiki/WeNNI
#
# WeNNI is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# WeNNI is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
# License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
# 02111-1307, USA.

use strict;
use FindBin;
use Getopt::Long;

sub basefile_options($$);
sub BaseFileConverter($);
sub cleanup($);
sub cmdline_options($);
sub nni($);
sub NNIFileConverter($$);
sub result($$);

# read command line options
my %cmdline_option;
cmdline_options(\%cmdline_option);

# write STDIN to a file
my $basefile="BASEfile.data";
open(BASEFILE,">$basefile");
while(<>) {
    print BASEFILE;
}
close(BASEFILE);

# determine beta and number of neighbours. Presedence order is (low to
# high priority): default value, basefile value, command line value.
my %option=( beta => 0.6, datatype => "derived",
             neighbours => 4, nodelete => 0 );
basefile_options(\%option,$basefile);
foreach my $key (keys %cmdline_option) {
    # override any options if set at command line
    if (defined($cmdline_option{$key})) {
        $option{$key}=$cmdline_option{$key};
    }
}

BaseFileConverter($option{'datatype'});
NNIFileConverter($option{'beta'},$option{'datatype'});
nni($option{'neighbours'});
result($basefile,$option{'datatype'});
if (!$option{'nodelete'}) {
    cleanup($basefile);
}

sub basefile_options($$) {
    my $option=$_[0];
    open(BASEFILE,"<$_[1]");
    my $parameter_section_found=0;
    while(<BASEFILE>) {
        chomp;
        if (/section\tWeNNIParams/) {
            $parameter_section_found=1;
            last;
        }
    }
    if ($parameter_section_found) {
        while(<BASEFILE>) {
            chomp;
            if (/%/) {
                last;
            }
            my @param=split "\t";
            $$option{$param[0]}=$param[1];
        }
    }
    close(BASEFILE);
}

sub BaseFileConverter($) {
    my $datatype=$_[0];
    my $command="$FindBin::Bin/BaseFileConverter";
    $command.=" BASEfile.data";
    $command.=" wenni_";
    if ($datatype eq "derived") {
      $command.=" -assayFields intensity1";
      $command.=" -assayFields intensity2";
    }
    else {
      $command.=" -assayFields FCh1Mean";
      $command.=" -assayFields BCh1Mean";
      $command.=" -assayFields FCh2Mean";
      $command.=" -assayFields BCh2Mean";
    }
    $command.=" -assayFields BCh1SD";
    $command.=" -assayFields BCh2SD";
    return system("$command");
}

sub cleanup($) {
    my $basefile=$_[0];
    unlink($basefile);
    unlink("wenni_BCh1Mean.data");
    unlink("wenni_BCh1SD.data");
    unlink("wenni_BCh2Mean.data");
    unlink("wenni_BCh2SD.data");
    unlink("wenni_FCh1Mean.data");
    unlink("wenni_FCh2Mean.data");
    unlink("wenni_imputed.data");
    unlink("wenni_intensity1.data");
    unlink("wenni_intensity2.data");
    unlink("wenni_logratio.data");
    unlink("wenni_weight.data");
}

sub cmdline_options($) {
    my $option=$_[0];
    my $help=0;
    my $result=GetOptions("beta=f",\$$option{beta},
              "datatype=s",\$$option{datatype},
              "help",\$help,
              "neighbours=i",\$$option{neighbours},
              "nodelete",\$$option{nodelete});
    if (defined($$option{datatype}) && ($$option{datatype} ne "derived") &&
        ($$option{datatype} ne "raw")) {
        $help=1;
    }
    if (!$result || $help) {
      help();
      exit(!$result);
    }
}

sub help() {
  print << "HELPTXT";

Usage wenni.pl [options]

Options
Options are read from the input BASE file, but the settings may be
changed with the corresponding command line option:

  --beta float
    Set beta value to float.
  --datatype string
    Allowed string values for datatype is 'raw' or
    'derived' (without '' characters). When 'derived' is
    given, the user defined intensity1 and intensity2 will
    be used from the Basefile. When 'raw' is used,
    intensity1 will be calculated as FCh1Mean-BCh1Mean
    (intensity2 is defined correspondingly).
  --neighbours integer
    Set the number of neighbours to integer.
  --nodelete
    Do not clean up, i.e. temporary files will not be deleted.
  --help
    Print this message.

This script runs WeNNI on data read from STDIN, and write the result
to STDOUT.

HELPTXT
}

sub nni($) {
    my $neighbours=$_[0];
    my $command="$FindBin::Bin/nni";
    $command.=" -data wenni_logratio.data";
    $command.=" -neighbours $neighbours";
    $command.=" -nni_algorithm WeNNI";
    $command.=" -weight wenni_weight.data";
    $command.=" > wenni_imputed.data";
    return system("$command");
}

sub NNIFileConverter($$) {
    my $beta=$_[0];
    my $datatype=$_[1];
    my $command="$FindBin::Bin/NNIFileConverter";
    $command.=" -beta $beta";
    if ($datatype eq "derived") {
      $command.=" -fg1 wenni_intensity1.data";
      $command.=" -fg2 wenni_intensity2.data";
      $command.=" -datatype derived";
    }
    else {
        $command.=" -bg1 wenni_BCh1Mean.data";
        $command.=" -bg2 wenni_BCh2Mean.data";
        $command.=" -fg1 wenni_FCh1Mean.data";
        $command.=" -fg2 wenni_FCh2Mean.data";
        $command.=" -datatype raw";
    }
    $command.=" -bgstd1 wenni_BCh1SD.data";
    $command.=" -bgstd2 wenni_BCh2SD.data";
    $command.=" -logratio wenni_logratio.data";
    $command.=" -weight wenni_weight.data";
    return system("$command");
}

sub result($$) {
    my $basefile=$_[0];
    my $datatype=$_[1];
    print "BASEfile\n";
    print "section\tspots\n";
    print "columns\tposition\treporter\tassayData\n";
    print "assayFields\tl2ratio1_2\tl10intgmean1_2\n";

    open(BASEFILE,"<$basefile");
    my $assays_found=0;
    while (<BASEFILE>) {
        if (/^assays\t/) {
            $assays_found=1;
            last;
        }
    }
    if (!$assays_found) {
        print STDERR "Error: Cannot find assays in 'section spot'\n";
        exit(-1);
    }
    print;

    my $percent_found=0;
    while (<BASEFILE>) {
        if (/^%/) {
            $percent_found=1;
            last;
        }
    }
    if (!$percent_found) {
        print STDERR "Error: Cannot find end of 'section spot' header\n";
        exit(-1);
    }
    print;

    my $log10=log(10);
    open(IMPUTEFILE,"<wenni_imputed.data");
    if ($datatype eq "raw") {
        open(B1FILE,"<wenni_BCh1Mean.data");
        open(F1FILE,"<wenni_FCh1Mean.data");
        open(B2FILE,"<wenni_BCh2Mean.data");
        open(F2FILE,"<wenni_FCh2Mean.data");
    }
    else {
        open(F1FILE,"<wenni_intensity1.data");
        open(F2FILE,"<wenni_intensity2.data");
    }
    while (<IMPUTEFILE>) {
        chop;
        my @impute_values=split(/\s/);
        my $row=<F1FILE>;
        chop;
        my @f1_values=split(/\s/,$row);
        $row=<F2FILE>;
        chop;
        my @f2_values=split(/\s/,$row);
        my @b1_values;
        my @b2_values;
        if ($datatype eq "raw") {
            $row=<B1FILE>;
            chop;
            @b1_values=split(/\s/,$row);
            $row=<B2FILE>;
            chop;
            @b2_values=split(/\s/,$row);
        }
        my $basefile_row=<BASEFILE>;
        my ($position,$reporter)=split("\t",$basefile_row);
        print "$position\t$reporter";
        my $columns=@impute_values;
        for (my $i=0; $i<$columns; $i++) {
            my $int1=$f1_values[$i]-($datatype eq "raw" ? $b1_values[$i] : 0);
            my $int2=$f2_values[$i]-($datatype eq "raw" ? $b2_values[$i] : 0);
            $int1=($int1<0 ? 0 : $int1);
            $int2=($int2<0 ? 0 : $int2);
            my $A=( ($int1>0) && ($int2>0) ? log(sqrt($int1*$int2))/$log10 : 0);
            print "\t$impute_values[$i]\t$A";
        }
        print "\n";
    }
    if ($datatype eq "raw") {
        close(B1FILE);
        close(B2FILE);
    }
    close(F1FILE);
    close(F2FILE);
    close(IMPUTEFILE);
    close(BASEFILE);
}