root/galaxy-central/tools/discreteWavelet/execute_dwt_cor_aVb_all.pl

リビジョン 2, 7.4 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
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);
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。