[2] | 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 | |
---|