[2] | 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 | |
---|
| 15 | # check to make sure having correct input and output files |
---|
| 16 | my $usage = "usage: compute_q_values.pl [TABULAR.in] [lambda] [pi0_method] [fdr_level] [robust] [TABULAR.out] [PDF.out] \n"; |
---|
| 17 | die $usage unless @ARGV == 7; |
---|
| 18 | |
---|
| 19 | #get the input arguments |
---|
| 20 | my $p_valuesInputFile = $ARGV[0]; |
---|
| 21 | my $lambdaValue = $ARGV[1]; |
---|
| 22 | my $pi0_method = $ARGV[2]; |
---|
| 23 | my $fdr_level = $ARGV[3]; |
---|
| 24 | my $robustValue = $ARGV[4]; |
---|
| 25 | my $q_valuesOutputFile = $ARGV[5]; |
---|
| 26 | my $p_q_values_histograms_QPlotsFile = $ARGV[6]; |
---|
| 27 | |
---|
| 28 | if($lambdaValue =~ /sequence/){ |
---|
| 29 | $lambdaValue = "seq(0, 0.95, 0.05)"; |
---|
| 30 | } |
---|
| 31 | |
---|
| 32 | #open the input files |
---|
| 33 | open (INPUT, "<", $p_valuesInputFile) || die("Could not open file $p_valuesInputFile \n"); |
---|
| 34 | open (OUTPUT1, ">", $q_valuesOutputFile) || die("Could not open file $q_valuesOutputFile \n"); |
---|
| 35 | open (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 |
---|
| 44 | my $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 | |
---|
| 50 | open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n"; |
---|
| 51 | print 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"; |
---|
| 83 | close Rcmd; |
---|
| 84 | |
---|
| 85 | system("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); |
---|
| 89 | close(OUTPUT2); |
---|
| 90 | close(OUTPUT1); |
---|
| 91 | close(INPUT); |
---|