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