root/galaxy-central/tools/regVariation/compute_q_values.pl

リビジョン 2, 3.7 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
1# A program to compute the q-values based on the p-values of multiple simultaneous tests.
2# The q-valules are computed using a specific R package created by John Storey called "qvalue".
3# The input is a TABULAR format file consisting of one column only that represents the p-values
4# of multiple simultaneous tests, one line for every p-value.
5# The first output is a TABULAR format file consisting of one column only that represents the q-values
6# corresponding to p-values, one line for every q-value.
7# the second output is a TABULAR format file consisting of three pages: the first page represents
8# the p-values histogram, the second page represents the q-values histogram, and the third page represents
9# the four Q-plots as introduced in the "qvalue" package manual.
10
11use strict;
12use warnings;
13use IO::Handle;
14
15# check to make sure having correct input and output files
16my $usage = "usage: compute_q_values.pl [TABULAR.in] [lambda] [pi0_method] [fdr_level] [robust] [TABULAR.out] [PDF.out] \n";
17die $usage unless @ARGV == 7;
18
19#get the input arguments
20my $p_valuesInputFile = $ARGV[0];
21my $lambdaValue =  $ARGV[1];
22my $pi0_method =  $ARGV[2];
23my $fdr_level =  $ARGV[3];
24my $robustValue =  $ARGV[4];
25my $q_valuesOutputFile = $ARGV[5];
26my $p_q_values_histograms_QPlotsFile = $ARGV[6];
27
28if($lambdaValue =~ /sequence/){
29        $lambdaValue = "seq(0, 0.95, 0.05)";
30}
31
32#open the input files
33open (INPUT, "<", $p_valuesInputFile) || die("Could not open file $p_valuesInputFile \n");
34open (OUTPUT1, ">", $q_valuesOutputFile) || die("Could not open file $q_valuesOutputFile \n");
35open (OUTPUT2, ">", $p_q_values_histograms_QPlotsFile) || die("Could not open file $p_q_values_histograms_QPlotsFile \n");
36#open (ERROR,  ">", "error.txt")  or die ("Could not open file error.txt \n");
37
38#save all error messages into the error file $errorFile using the error file handle ERROR
39#STDERR -> fdopen( \*ERROR,  "w" ) or die ("Could not direct errors to the error file error.txt \n");
40
41#warn "Hello Error File \n";
42
43#variable to store the name of the R script file
44my $r_script;
45
46# R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values       
47# construct an R script file and save it in the same directory where the perl file is located
48$r_script = "q_values_computation.r";
49
50open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
51print Rcmd "
52        #options(show.error.messages = FALSE);
53       
54        #load necessary packages
55        suppressPackageStartupMessages(library(tcltk));
56        library(qvalue);
57       
58        #read the p-values of the multiple simultaneous tests from the input file $p_valuesInputFile
59        p <- scan(\"$p_valuesInputFile\", quiet = TRUE);
60       
61        #compute the q-values that correspond to the p-values of the multiple simultaneous tests
62        qobj <- qvalue(p, pi0.meth = \"$pi0_method\", lambda = $lambdaValue, fdr.level = $fdr_level, robust = $robustValue);
63        #qobj <- qvalue(p, pi0.meth = \"smoother\", lambda = seq(0, 0.95, 0.05), fdr.level = 0.05);
64        #qobj <- qvalue(p, pi0.meth = \"bootstrap\", fdr.level = 0.05);
65       
66        #draw the p-values histogram, the q-values histogram, and the four Q-plots
67        # and save them on multiple pages of the output file $p_q_values_histograms_QPlotsFile
68        pdf(file = \"$p_q_values_histograms_QPlotsFile\", width = 6.25, height = 6, family = \"Times\", pointsize = 12, onefile = TRUE)
69        hist(qobj\$pvalues);
70        #dev.off();
71       
72        hist(qobj\$qvalues);
73        #dev.off();
74       
75        qplot(qobj); 
76        dev.off();
77       
78        #save the q-values in the output file $q_valuesOutputFile
79        qwrite(qobj, filename=\"$q_valuesOutputFile\");
80
81        #options(show.error.messages = TRUE);
82        #eof\n";
83close Rcmd;     
84
85system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
86
87#close the input and output and error files
88#close(ERROR);
89close(OUTPUT2);
90close(OUTPUT1);
91close(INPUT);
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。