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

Last change on this file since 115 was 115, checked in by Jari Häkkinen, 17 years ago

Well defined A values are now kept in BASE.

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 7.6 KB
Line 
1#!/usr/bin/perl -w
2
3# $Id: wenni.pl 115 2006-06-22 08:43:07Z jari $
4
5# Copyright (C) 2005, 2006 Jari Häkkinen
6#
7# This file is part of WeNNI,
8# http://lev.thep.lu.se/trac/baseplugins/wiki/WeNNI
9#
10# WeNNI is free software; you can redistribute it and/or modify it
11# under the terms of the GNU General Public License as published by
12# the Free Software Foundation; either version 2 of the License, or
13# (at your option) any later version.
14#
15# WeNNI is distributed in the hope that it will be useful, but WITHOUT
16# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
18# License for more details.
19#
20# You should have received a copy of the GNU General Public License
21# along with this program; if not, write to the Free Software
22# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
23# 02111-1307, USA.
24
25use strict;
26use FindBin;
27use Getopt::Long;
28
29sub basefile_options($$);
30sub BaseFileConverter($);
31sub cleanup($);
32sub cmdline_options($);
33sub nni($);
34sub NNIFileConverter($$);
35sub result($$);
36
37# read command line options
38my %cmdline_option;
39cmdline_options(\%cmdline_option);
40
41# write STDIN to a file
42my $basefile="BASEfile.data";
43open(BASEFILE,">$basefile");
44while(<>) {
45    print BASEFILE;
46}
47close(BASEFILE);
48
49# determine beta and number of neighbours. Presedence order is (low to
50# high priority): default value, basefile value, command line value.
51my %option=( beta => 0.6, datatype => "derived",
52             neighbours => 4, nodelete => 0 );
53basefile_options(\%option,$basefile);
54foreach my $key (keys %cmdline_option) {
55    # override any options if set at command line
56    if (defined($cmdline_option{$key})) {
57        $option{$key}=$cmdline_option{$key};
58    }
59}
60
61BaseFileConverter($option{'datatype'});
62NNIFileConverter($option{'beta'},$option{'datatype'});
63nni($option{'neighbours'});
64result($basefile,$option{'datatype'});
65if (!$option{'nodelete'}) {
66    cleanup($basefile);
67}
68
69sub basefile_options($$) {
70    my $option=$_[0];
71    open(BASEFILE,"<$_[1]");
72    my $parameter_section_found=0;
73    while(<BASEFILE>) {
74        chomp;
75        if (/section\tWeNNIParams/) {
76            $parameter_section_found=1;
77            last;
78        }
79    }
80    if ($parameter_section_found) {
81        while(<BASEFILE>) {
82            chomp;
83            if (/%/) {
84                last;
85            }
86            my @param=split "\t";
87            $$option{$param[0]}=$param[1];
88        }
89    }
90    close(BASEFILE);
91}
92
93sub BaseFileConverter($) {
94    my $datatype=$_[0];
95    my $command="$FindBin::Bin/BaseFileConverter";
96    $command.=" BASEfile.data";
97    $command.=" wenni_";
98    if ($datatype eq "derived") {
99      $command.=" -assayFields intensity1";
100      $command.=" -assayFields intensity2";
101    }
102    else {
103      $command.=" -assayFields FCh1Mean";
104      $command.=" -assayFields BCh1Mean";
105      $command.=" -assayFields FCh2Mean";
106      $command.=" -assayFields BCh2Mean";
107    }
108    $command.=" -assayFields BCh1SD";
109    $command.=" -assayFields BCh2SD";
110    return system("$command");
111}
112
113sub cleanup($) {
114    my $basefile=$_[0];
115    unlink($basefile);
116    unlink("wenni_BCh1Mean.data");
117    unlink("wenni_BCh1SD.data");
118    unlink("wenni_BCh2Mean.data");
119    unlink("wenni_BCh2SD.data");
120    unlink("wenni_FCh1Mean.data");
121    unlink("wenni_FCh2Mean.data");
122    unlink("wenni_imputed.data");
123    unlink("wenni_intensity1.data");
124    unlink("wenni_intensity2.data");
125    unlink("wenni_logratio.data");
126    unlink("wenni_weight.data");
127}
128
129sub cmdline_options($) {
130    my $option=$_[0];
131    my $help=0;
132    my $result=GetOptions("beta=f",\$$option{beta},
133              "datatype=s",\$$option{datatype},
134              "help",\$help,
135              "neighbours=i",\$$option{neighbours},
136              "nodelete",\$$option{nodelete});
137    if (defined($$option{datatype}) && ($$option{datatype} ne "derived") &&
138        ($$option{datatype} ne "raw")) {
139        $help=1;
140    }
141    if (!$result || $help) {
142      help();
143      exit(!$result);
144    }
145}
146
147sub help() {
148  print << "HELPTXT";
149
150Usage wenni.pl [options]
151
152Options
153Options are read from the input BASE file, but the settings may be
154changed with the corresponding command line option:
155
156  --beta float
157    Set beta value to float.
158  --datatype string
159    Allowed string values for datatype is 'raw' or
160    'derived' (without '' characters). When 'derived' is
161    given, the user defined intensity1 and intensity2 will
162    be used from the Basefile. When 'raw' is used,
163    intensity1 will be calculated as FCh1Mean-BCh1Mean
164    (intensity2 is defined correspondingly).
165  --neighbours integer
166    Set the number of neighbours to integer.
167  --nodelete
168    Do not clean up, i.e. temporary files will not be deleted.
169  --help
170    Print this message.
171
172This script runs WeNNI on data read from STDIN, and write the result
173to STDOUT.
174
175HELPTXT
176}
177
178sub nni($) {
179    my $neighbours=$_[0];
180    my $command="$FindBin::Bin/nni";
181    $command.=" -data wenni_logratio.data";
182    $command.=" -neighbours $neighbours";
183    $command.=" -nni_algorithm WeNNI";
184    $command.=" -weight wenni_weight.data";
185    $command.=" > wenni_imputed.data";
186    return system("$command");
187}
188
189sub NNIFileConverter($$) {
190    my $beta=$_[0];
191    my $datatype=$_[1];
192    my $command="$FindBin::Bin/NNIFileConverter";
193    $command.=" -beta $beta";
194    if ($datatype eq "derived") {
195      $command.=" -fg1 wenni_intensity1.data";
196      $command.=" -fg2 wenni_intensity2.data";
197      $command.=" -datatype derived";
198    }
199    else {
200        $command.=" -bg1 wenni_BCh1Mean.data";
201        $command.=" -bg2 wenni_BCh2Mean.data";
202        $command.=" -fg1 wenni_FCh1Mean.data";
203        $command.=" -fg2 wenni_FCh2Mean.data";
204        $command.=" -datatype raw";
205    }
206    $command.=" -bgstd1 wenni_BCh1SD.data";
207    $command.=" -bgstd2 wenni_BCh2SD.data";
208    $command.=" -logratio wenni_logratio.data";
209    $command.=" -weight wenni_weight.data";
210    return system("$command");
211}
212
213sub result($$) {
214    my $basefile=$_[0];
215    my $datatype=$_[1];
216    print "BASEfile\n";
217    print "section\tspots\n";
218    print "columns\tposition\treporter\tassayData\n";
219    print "assayFields\tl2ratio1_2\tl10intgmean1_2\n";
220
221    open(BASEFILE,"<$basefile");
222    my $assays_found=0;
223    while (<BASEFILE>) {
224        if (/^assays\t/) {
225            $assays_found=1;
226            last;
227        }
228    }
229    if (!$assays_found) {
230        print STDERR "Error: Cannot find assays in 'section spot'\n";
231        exit(-1);
232    }
233    print;
234
235    my $percent_found=0;
236    while (<BASEFILE>) {
237        if (/^%/) {
238            $percent_found=1;
239            last;
240        }
241    }
242    if (!$percent_found) {
243        print STDERR "Error: Cannot find end of 'section spot' header\n";
244        exit(-1);
245    }
246    print;
247
248    my $log10=log(10);
249    open(IMPUTEFILE,"<wenni_imputed.data");
250    if ($datatype eq "raw") {
251        open(B1FILE,"<wenni_BCh1Mean.data");
252        open(F1FILE,"<wenni_FCh1Mean.data");
253        open(B2FILE,"<wenni_BCh2Mean.data");
254        open(F2FILE,"<wenni_FCh2Mean.data");
255    }
256    else {
257        open(F1FILE,"<wenni_intensity1.data");
258        open(F2FILE,"<wenni_intensity2.data");
259    }
260    while (<IMPUTEFILE>) {
261        chop;
262        my @impute_values=split(/\s/);
263        my $row=<F1FILE>;
264        chop;
265        my @f1_values=split(/\s/,$row);
266        $row=<F2FILE>;
267        chop;
268        my @f2_values=split(/\s/,$row);
269        my @b1_values;
270        my @b2_values;
271        if ($datatype eq "raw") {
272            $row=<B1FILE>;
273            chop;
274            @b1_values=split(/\s/,$row);
275            $row=<B2FILE>;
276            chop;
277            @b2_values=split(/\s/,$row);
278        }
279        my $basefile_row=<BASEFILE>;
280        my ($position,$reporter)=split("\t",$basefile_row);
281        print "$position\t$reporter";
282        my $columns=@impute_values;
283        for (my $i=0; $i<$columns; $i++) {
284            my $int1=$f1_values[$i]-($datatype eq "raw" ? $b1_values[$i] : 0);
285            my $int2=$f2_values[$i]-($datatype eq "raw" ? $b2_values[$i] : 0);
286            $int1=($int1<0 ? 0 : $int1);
287            $int2=($int2<0 ? 0 : $int2);
288            my $A=( ($int1>0) && ($int2>0) ? log(sqrt($int1*$int2))/$log10 : 0);
289            print "\t$impute_values[$i]\t$A";
290        }
291        print "\n";
292    }
293    if ($datatype eq "raw") {
294        close(B1FILE);
295        close(B2FILE);
296    }
297    close(F1FILE);
298    close(F2FILE);
299    close(IMPUTEFILE);
300    close(BASEFILE);
301}
Note: See TracBrowser for help on using the repository browser.