annotate tools/discreteWavelet/execute_dwt_IvC_all.pl @ 2:c2a356708570

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:45:42 -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 #!/usr/bin/perl -w
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 use warnings;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 use IO::Handle;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 $usage = "execute_dwt_IvC_all.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] [PDF.out] \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 die $usage unless @ARGV == 4;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 #get the input arguments
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 my $firstInputFile = $ARGV[0];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 my $secondInputFile = $ARGV[1];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11 my $firstOutputFile = $ARGV[2];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 my $secondOutputFile = $ARGV[3];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 open (INPUT1, "<", $firstInputFile) || die("Could not open file $firstInputFile \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 open (INPUT2, "<", $secondInputFile) || die("Could not open file $secondInputFile \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 open (OUTPUT1, ">", $firstOutputFile) || die("Could not open file $firstOutputFile \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 open (OUTPUT2, ">", $secondOutputFile) || die("Could not open file $secondOutputFile \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 open (ERROR, ">", "error.txt") or die ("Could not open file error.txt \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 #save all error messages into the error file $errorFile using the error file handle ERROR
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 STDERR -> fdopen( \*ERROR, "w" ) or die ("Could not direct errors to the error file error.txt \n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 print "There are two input data files: \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25 print "The input data file is: $firstInputFile \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 print "The control data file is: $secondInputFile \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 # IvC test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29 $test = "IvC";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 # construct an R script to implement the IvC test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 print "\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 $r_script = "get_dwt_IvC_test.r";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 print "$r_script \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 # R script
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 print Rcmd "
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40 ###########################################################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 # code to do wavelet Indel vs. Control
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 # signal is the difference I-C; function is second moment i.e. variance from zero not mean
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 # to perform wavelet transf. of signal, scale-by-scale analysis of the function
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 # create null bands by permuting the original data series
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 # generate plots and table matrix of correlation coefficients including p-values
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 ############################################################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 library(\"Rwave\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 library(\"wavethresh\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 library(\"waveslim\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 options(echo = FALSE)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53 # normalize data
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 norm <- function(data){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 v <- (data - mean(data))/sd(data);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 if(sum(is.na(v)) >= 1){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 v <- data;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 return(v);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 dwt_cor <- function(data.short, names.short, data.long, names.long, test, pdf, table, filter = 4, bc = \"symmetric\", wf = \"haar\", boundary = \"reflection\") {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 print(test);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 print(pdf);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65 print(table);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 pdf(file = pdf);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 final_pvalue = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 title = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 short.levels <- wd(data.short[, 1], filter.number = filter, bc = bc)\$nlevels;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 title <- c(\"motif\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 for (i in 1:short.levels){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 title <- c(title, paste(i, \"moment2\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"), paste(i, \"test\", sep = \"_\"));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 print(title);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 # loop to compare a vs a
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 for(i in 1:length(names.short)){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 wave1.dwt = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 m2.dwt = diff = var.dwt = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 out = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 out <- vector(length = length(title));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 print(names.short[i]);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 print(names.long[i]);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 # need exit if not comparing motif(a) vs motif(a)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 if (names.short[i] != names.long[i]){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 stop(paste(\"motif\", names.short[i], \"is not the same as\", names.long[i], sep = \" \"));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 else {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 # signal is the difference I-C data sets
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 diff<-data.short[,i]-data.long[,i];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 # normalize the signal
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 diff<-norm(diff);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99 # function is 2nd moment
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 # 2nd moment m_j = 1/N[sum_N(W_j + V_J)^2] = 1/N sum_N(W_j)^2 + (X_bar)^2
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 wave1.dwt <- dwt(diff, wf = wf, short.levels, boundary = boundary);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 var.dwt <- wave.variance(wave1.dwt);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 m2.dwt <- vector(length = short.levels)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 for(level in 1:short.levels){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 m2.dwt[level] <- var.dwt[level, 1] + (mean(diff)^2);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 # CI bands by permutation of time series
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 feature1 = feature2 = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 feature1 = data.short[, i];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111 feature2 = data.long[, i];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112 null = results = med = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 m2_25 = m2_975 = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115 for (k in 1:1000) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 nk_1 = nk_2 = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117 m2_null = var_null = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118 null.levels = null_wave1 = null_diff = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119 nk_1 <- sample(feature1, length(feature1), replace = FALSE);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 nk_2 <- sample(feature2, length(feature2), replace = FALSE);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 null.levels <- wd(nk_1, filter.number = filter, bc = bc)\$nlevels;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 null_diff <- nk_1-nk_2;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 null_diff <- norm(null_diff);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124 null_wave1 <- dwt(null_diff, wf = wf, short.levels, boundary = boundary);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125 var_null <- wave.variance(null_wave1);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 m2_null <- vector(length = null.levels);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 for(level in 1:null.levels){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 m2_null[level] <- var_null[level, 1] + (mean(null_diff)^2);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 null= rbind(null, m2_null);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 null <- apply(null, 2, sort, na.last = TRUE);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 m2_25 <- null[25,];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 m2_975 <- null[975,];
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 med <- apply(null, 2, median, na.rm = TRUE);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 # plot
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 results <- cbind(m2.dwt, m2_25, m2_975);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2), xlab = \"Wavelet Scale\", ylab = c(\"Wavelet 2nd Moment\", test), main = (names.short[i]), cex.main = 0.75);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 abline(h = 1);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 # get pvalues by comparison to null distribution
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 out <- c(names.short[i]);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 for (m in 1:length(m2.dwt)){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 print(paste(\"scale\", m, sep = \" \"));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 print(paste(\"m2\", m2.dwt[m], sep = \" \"));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 print(paste(\"median\", med[m], sep = \" \"));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 out <- c(out, format(m2.dwt[m], digits = 4));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150 pv = NULL;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151 if(is.na(m2.dwt[m])){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152 pv <- \"NA\";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 else {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 if (m2.dwt[m] >= med[m]){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156 # R tail test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157 tail <- \"R\";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158 pv <- (length(which(null[, m] >= m2.dwt[m])))/(length(na.exclude(null[, m])));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160 else{
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
161 if (m2.dwt[m] < med[m]){
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
162 # L tail test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
163 tail <- \"L\";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
164 pv <- (length(which(null[, m] <= m2.dwt[m])))/(length(na.exclude(null[, m])));
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
165 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
166 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
167 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
168 out <- c(out, pv);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
169 print(pv);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
170 out <- c(out, tail);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
171 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
172 final_pvalue <-rbind(final_pvalue, out);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
173 print(out);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
174 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
175 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
176
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
177 colnames(final_pvalue) <- title;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
178 write.table(final_pvalue, file = table, sep = \"\\t\", quote = FALSE, row.names = FALSE);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
179 dev.off();
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
180 }\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
181
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
182 print Rcmd "
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
183 # execute
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
184 # read in data
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
185
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
186 inputData <- read.delim(\"$firstInputFile\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
187 inputDataNames <- colnames(inputData);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
188
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
189 controlData <- read.delim(\"$secondInputFile\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
190 controlDataNames <- colnames(controlData);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
191
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
192 # call the test function to implement IvC test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
193 dwt_cor(inputData, inputDataNames, controlData, controlDataNames, test = \"$test\", pdf = \"$secondOutputFile\", table = \"$firstOutputFile\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
194 print (\"done with the correlation test\");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
195 \n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
196
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
197 print Rcmd "#eof\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
198
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
199 close Rcmd;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
200
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
201 system("echo \"wavelet IvC test started on \`hostname\` at \`date\`\"\n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
202 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out\n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
203 system("echo \"wavelet IvC test ended on \`hostname\` at \`date\`\"\n");
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
204
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
205 #close the input and output and error files
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
206 close(ERROR);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
207 close(OUTPUT2);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
208 close(OUTPUT1);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
209 close(INPUT2);
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
210 close(INPUT1);