/tools/discreteWavelet/execute_dwt_cor_aVa_perClass.pl
Perl | 221 lines | 156 code | 36 blank | 29 comment | 22 complexity | d3ad4f072a2925c11263dfc8d541cd71 MD5 | raw file
1#!/usr/bin/perl -w 2 3use warnings; 4use IO::Handle; 5 6$usage = "execute_dwt_cor_aVa_perClass.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] [PDF.out] \n"; 7die $usage unless @ARGV == 4; 8 9#get the input arguments 10my $firstInputFile = $ARGV[0]; 11my $secondInputFile = $ARGV[1]; 12my $firstOutputFile = $ARGV[2]; 13my $secondOutputFile = $ARGV[3]; 14 15open (INPUT1, "<", $firstInputFile) || die("Could not open file $firstInputFile \n"); 16open (INPUT2, "<", $secondInputFile) || die("Could not open file $secondInputFile \n"); 17open (OUTPUT1, ">", $firstOutputFile) || die("Could not open file $firstOutputFile \n"); 18open (OUTPUT2, ">", $secondOutputFile) || die("Could not open file $secondOutputFile \n"); 19open (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 22STDERR -> fdopen( \*ERROR, "w" ) or die ("Could not direct errors to the error file error.txt \n"); 23 24print "There are two input data files: \n"; 25print "The input data file is: $firstInputFile \n"; 26print "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 32print "\n"; 33 34$r_script = "get_dwt_cor_aVa_test.r"; 35print "$r_script \n"; 36 37open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n"; 38print 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 188print 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"; 209close Rcmd; 210 211system("echo \"wavelet IvC test started on \`hostname\` at \`date\`\"\n"); 212system("R --no-restore --no-save --no-readline < $r_script > $r_script.out\n"); 213system("echo \"wavelet IvC test ended on \`hostname\` at \`date\`\"\n"); 214 215#close the input and output and error files 216close(ERROR); 217close(OUTPUT2); 218close(OUTPUT1); 219close(INPUT2); 220close(INPUT1); 221