annotate tools/regVariation/compute_q_values.pl @ 1:cdcb0ce84a1b

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