Mercurial > repos > devteam > compute_q_values

comparison compute_q_values.pl @ 0:63abfc350814

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:53:54 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:63abfc350814
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);