0
|
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
|