/tools/discreteWavelet/execute_dwt_cor_aVb_all.pl

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