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