/tools/discreteWavelet/execute_dwt_cor_aVb_all.pl

https://bitbucket.org/cistrome/cistrome-harvard/ · Perl · 223 lines · 202 code · 14 blank · 7 comment · 5 complexity · 15f4c23579861176aa256c51a9e48f83 MD5 · raw file

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