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