/tools/discreteWavelet/execute_dwt_IvC_all.pl

https://bitbucket.org/cistrome/cistrome-harvard/ · Perl · 210 lines · 188 code · 15 blank · 7 comment · 5 complexity · d6f92d43152dc5ebe684296c4b3065cc MD5 · raw file

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