|
5
|
1 ###############################################################################
|
|
|
2 # In summary this script calculates the z prime factor for one or more plates
|
|
|
3 # of data.
|
|
|
4 #
|
|
|
5 # Args:
|
|
|
6 # input.data.frame contains one or more columns of plate data. The first column
|
|
|
7 # must contain the well annotation ie A01, A02 etc. Each subsequent column represents
|
|
|
8 # a plate of data where the column name is the plate name.
|
|
|
9 #
|
|
|
10 # plate.conf.data.frame is a file of the following format:
|
|
|
11 # well, type(syntax must be either "poscontr" or "negcontr"), name (name of control)
|
|
|
12 #
|
|
|
13 # Returns:
|
|
|
14 # z'factor for each plate for each combination of each type of positive vs negative control.
|
|
|
15 #
|
|
|
16 # Author: gouldkate
|
|
|
17 ###############################################################################
|
|
|
18
|
|
|
19 use strict;
|
|
|
20 use warnings;
|
|
|
21 use IO::Handle;
|
|
|
22 use File::Temp qw/ tempfile tempdir /;
|
|
|
23 my $tdir = tempdir( CLEANUP => 0 );
|
|
|
24
|
|
|
25 # check to make sure having correct input and output files
|
|
|
26 my $usage = "usage: calculateZPrimeFactor.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] \n";
|
|
|
27 die $usage unless @ARGV == 3;
|
|
|
28
|
|
|
29 #get the input arguments
|
|
|
30 my $plateData = $ARGV[0];
|
|
|
31 my $plateConfig = $ARGV[1];
|
|
|
32 my $zPrimeFactor = $ARGV[2];
|
|
|
33
|
|
|
34 #open the input files
|
|
|
35 open (INPUT1, "<", $plateData) || die("Could not open file $plateData \n");
|
|
|
36 open (INPUT2, "<", $plateConfig) || die("Could not open file $plateConfig \n");
|
|
|
37 open (OUTPUT1, ">", $zPrimeFactor) || die("Could not open file $zPrimeFactor \n");
|
|
|
38
|
|
|
39 #variable to store the name of the R script file
|
|
|
40 my $r_script;
|
|
|
41
|
|
|
42 # R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values
|
|
|
43 # construct an R script file and save it in a temp directory
|
|
|
44 chdir $tdir;
|
|
|
45 $r_script = "calculateZPrimeFactor.r";
|
|
|
46
|
|
|
47 open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
|
|
|
48 print Rcmd "
|
|
|
49 #options(show.error.messages = FALSE);
|
|
|
50 input.data.frame <- read.table(\"$plateData\", head=T, sep=\"\\t\", comment=\"\");
|
|
|
51 plate.conf.data.frame <- read.table(\"$plateConfig\", head=T, sep=\"\\t\", comment=\"\");
|
|
|
52
|
|
|
53 # assumed second column is type ie negcontr or poscontr and third column is name of control
|
|
|
54 negative.controls.list <- unique(plate.conf.data.frame[which (plate.conf.data.frame[,2] == \"negcontr\"),][,3])
|
|
|
55 positive.controls.list <- unique(plate.conf.data.frame[which (plate.conf.data.frame[,2] == \"poscontr\"),][,3])
|
|
|
56
|
|
|
57 z.prime.factor.report.data.frame <- data.frame()
|
|
|
58
|
|
|
59
|
|
|
60 for (negative.control in negative.controls.list){
|
|
|
61 for (positive.control in positive.controls.list){
|
|
|
62 for (i in 2:length(colnames(input.data.frame))){
|
|
|
63
|
|
|
64 negative.control.wells <- plate.conf.data.frame[which (plate.conf.data.frame[3] == negative.control),][,1]
|
|
|
65 positive.control.wells <- plate.conf.data.frame[which (plate.conf.data.frame[3] == positive.control),][,1]
|
|
|
66
|
|
|
67 control.duo <- paste(negative.control, positive.control, sep=\"/\")
|
|
|
68
|
|
|
69 negative.control.values <- input.data.frame[((input.data.frame[,1] %in% negative.control.wells)&(!(is.na(input.data.frame[,i])))),][,i]
|
|
|
70 positive.control.values <- input.data.frame[((input.data.frame[,1] %in% positive.control.wells)&(!(is.na(input.data.frame[,i])))),][,i]
|
|
|
71
|
|
|
72 if ((length(negative.control.values)==0)|(length(positive.control.values)==0)){
|
|
|
73 z.prime.factor <- NA
|
|
|
74 } else {
|
|
|
75 if ((sum(negative.control.values)) < (sum(positive.control.values))){
|
|
|
76 low.value.controls <- negative.control.values
|
|
|
77 high.value.controls <- positive.control.values
|
|
|
78 } else {
|
|
|
79 low.value.controls <- positive.control.values
|
|
|
80 high.value.controls <- negative.control.values
|
|
|
81 }
|
|
|
82 z.prime.factor <- round(1 - ((3 * (sd(high.value.controls, na.rm = TRUE) + sd(low.value.controls, na.rm = TRUE))/(abs(mean(high.value.controls, na.rm = TRUE) - mean(low.value.controls, na.rm = TRUE))))), 2)
|
|
|
83 }
|
|
|
84
|
|
|
85 plate.name <- colnames(input.data.frame)[i]
|
|
|
86 z.prime.factor.current.row.data.frame <- data.frame(plate=plate.name, control.duo=control.duo, z.prime.factor=z.prime.factor, stringsAsFactors=FALSE)
|
|
|
87 z.prime.factor.report.data.frame <-rbind(z.prime.factor.report.data.frame, z.prime.factor.current.row.data.frame)
|
|
|
88 }
|
|
|
89 }
|
|
|
90 }
|
|
|
91 z.prime.factor.report.data.frame <- z.prime.factor.report.data.frame[order(z.prime.factor.report.data.frame\$control.duo),]
|
|
|
92 write.table(z.prime.factor.report.data.frame, file=\"$zPrimeFactor\", quote=F, sep=\"\\t\", row.names=F);
|
|
|
93 #eof\n";
|
|
|
94
|
|
|
95 close Rcmd;
|
|
|
96
|
|
|
97 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
|
|
|
98
|
|
|
99 #close the input and output files
|
|
|
100 close(OUTPUT1);
|
|
|
101 close(INPUT1);
|
|
|
102 close(INPUT2);
|
