|
0
|
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
|
|
|
11 use strict;
|
|
|
12 use warnings;
|
|
|
13 use IO::Handle;
|
|
|
14 use File::Temp qw/ tempfile tempdir /;
|
|
|
15 my $tdir = tempdir( CLEANUP => 0 );
|
|
|
16
|
|
|
17 # check to make sure having correct input and output files
|
|
|
18 my $usage = "usage: compute_q_values.pl [TABULAR.in] [lambda] [pi0_method] [fdr_level] [robust] [TABULAR.out] [PDF.out] \n";
|
|
|
19 die $usage unless @ARGV == 7;
|
|
|
20
|
|
|
21 #get the input arguments
|
|
|
22 my $p_valuesInputFile = $ARGV[0];
|
|
|
23 my $lambdaValue = $ARGV[1];
|
|
|
24 my $pi0_method = $ARGV[2];
|
|
|
25 my $fdr_level = $ARGV[3];
|
|
|
26 my $robustValue = $ARGV[4];
|
|
|
27 my $q_valuesOutputFile = $ARGV[5];
|
|
|
28 my $p_q_values_histograms_QPlotsFile = $ARGV[6];
|
|
|
29
|
|
|
30 if($lambdaValue =~ /sequence/){
|
|
|
31 $lambdaValue = "seq(0, 0.95, 0.05)";
|
|
|
32 }
|
|
|
33
|
|
|
34 #open the input files
|
|
|
35 open (INPUT, "<", $p_valuesInputFile) || die("Could not open file $p_valuesInputFile \n");
|
|
|
36 open (OUTPUT1, ">", $q_valuesOutputFile) || die("Could not open file $q_valuesOutputFile \n");
|
|
|
37 open (OUTPUT2, ">", $p_q_values_histograms_QPlotsFile) || die("Could not open file $p_q_values_histograms_QPlotsFile \n");
|
|
|
38 #open (ERROR, ">", "error.txt") or die ("Could not open file error.txt \n");
|
|
|
39
|
|
|
40 #save all error messages into the error file $errorFile using the error file handle ERROR
|
|
|
41 #STDERR -> fdopen( \*ERROR, "w" ) or die ("Could not direct errors to the error file error.txt \n");
|
|
|
42
|
|
|
43 #warn "Hello Error File \n";
|
|
|
44
|
|
|
45 #variable to store the name of the R script file
|
|
|
46 my $r_script;
|
|
|
47
|
|
|
48 # R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values
|
|
|
49 # construct an R script file and save it in a temp directory
|
|
|
50 chdir $tdir;
|
|
|
51 $r_script = "q_values_computation.r";
|
|
|
52
|
|
|
53 open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
|
|
|
54 print Rcmd "
|
|
|
55 #options(show.error.messages = FALSE);
|
|
|
56
|
|
|
57 #load necessary packages
|
|
|
58 suppressPackageStartupMessages(library(tcltk));
|
|
|
59 library(qvalue);
|
|
|
60
|
|
|
61 #read the p-values of the multiple simultaneous tests from the input file $p_valuesInputFile
|
|
|
62 p <- scan(\"$p_valuesInputFile\", quiet = TRUE);
|
|
|
63
|
|
|
64 #compute the q-values that correspond to the p-values of the multiple simultaneous tests
|
|
|
65 qobj <- qvalue(p, pi0.meth = \"$pi0_method\", lambda = $lambdaValue, fdr.level = $fdr_level, robust = $robustValue);
|
|
|
66 #qobj <- qvalue(p, pi0.meth = \"smoother\", lambda = seq(0, 0.95, 0.05), fdr.level = 0.05);
|
|
|
67 #qobj <- qvalue(p, pi0.meth = \"bootstrap\", fdr.level = 0.05);
|
|
|
68
|
|
|
69 #draw the p-values histogram, the q-values histogram, and the four Q-plots
|
|
|
70 # and save them on multiple pages of the output file $p_q_values_histograms_QPlotsFile
|
|
|
71 pdf(file = \"$p_q_values_histograms_QPlotsFile\", width = 6.25, height = 6, family = \"Times\", pointsize = 12, onefile = TRUE)
|
|
|
72 hist(qobj\$pvalues);
|
|
|
73 #dev.off();
|
|
|
74
|
|
|
75 hist(qobj\$qvalues);
|
|
|
76 #dev.off();
|
|
|
77
|
|
|
78 qplot(qobj);
|
|
|
79 dev.off();
|
|
|
80
|
|
|
81 #save the q-values in the output file $q_valuesOutputFile
|
|
|
82 qobj\$pi0 <- signif(qobj\$pi0,digits=6)
|
|
|
83 qwrite(qobj, filename=\"$q_valuesOutputFile\");
|
|
|
84
|
|
|
85 #options(show.error.messages = TRUE);
|
|
|
86 #eof\n";
|
|
|
87 close Rcmd;
|
|
|
88
|
|
|
89 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
|
|
|
90
|
|
|
91 #close the input and output and error files
|
|
|
92 #close(ERROR);
|
|
|
93 close(OUTPUT2);
|
|
|
94 close(OUTPUT1);
|
|
|
95 close(INPUT);
|
