PageRenderTime 29ms CodeModel.GetById 17ms app.highlight 9ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/discreteWavelet/execute_dwt_cor_aVa_perClass.pl

https://bitbucket.org/cistrome/cistrome-harvard/
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