Mercurial > repos > xuebing > sharplabtool
comparison tools/regVariation/compute_q_values.pl @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:9071e359b9a3 |
---|---|
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); |