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

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

import galaxy-central

行番号 
1#!/usr/bin/perl -w
2use warnings;
3use IO::Handle;
4
5$usage = "execute_dwt_IvC_all.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] [PDF.out]  \n";
6die $usage unless @ARGV == 4;
7
8#get the input arguments
9my $firstInputFile = $ARGV[0];
10my $secondInputFile = $ARGV[1];
11my $firstOutputFile = $ARGV[2];
12my $secondOutputFile = $ARGV[3];
13
14open (INPUT1, "<", $firstInputFile) || die("Could not open file $firstInputFile \n");
15open (INPUT2, "<", $secondInputFile) || die("Could not open file $secondInputFile \n");
16open (OUTPUT1, ">", $firstOutputFile) || die("Could not open file $firstOutputFile \n");
17open (OUTPUT2, ">", $secondOutputFile) || die("Could not open file $secondOutputFile \n");
18open (ERROR,  ">", "error.txt")  or die ("Could not open file error.txt \n");
19
20#save all error messages into the error file $errorFile using the error file handle ERROR
21STDERR -> fdopen( \*ERROR,  "w" ) or die ("Could not direct errors to the error file error.txt \n");
22
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 = "IvC";
30
31# construct an R script to implement the IvC test
32print "\n";
33
34$r_script = "get_dwt_IvC_test.r";
35print "$r_script \n";
36
37# R script
38open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
39print Rcmd "
40        ###########################################################################################
41        # code to do wavelet Indel vs. Control
42        # signal is the difference I-C; function is second moment i.e. variance from zero not mean
43        # to perform wavelet transf. of signal, scale-by-scale analysis of the function
44        # create null bands by permuting the original data series
45        # generate plots and table matrix of correlation coefficients including p-values
46        ############################################################################################
47        library(\"Rwave\");
48        library(\"wavethresh\");
49        library(\"waveslim\");
50       
51        options(echo = FALSE)
52       
53        # normalize data
54        norm <- function(data){
55            v <- (data - mean(data))/sd(data);
56            if(sum(is.na(v)) >= 1){
57                v <- data;
58            }
59            return(v);
60        }
61       
62        dwt_cor <- function(data.short, names.short, data.long, names.long, test, pdf, table, filter = 4, bc = \"symmetric\", wf = \"haar\", boundary = \"reflection\") {
63            print(test);
64            print(pdf);
65            print(table);
66           
67            pdf(file = pdf);
68            final_pvalue = NULL;
69            title = NULL;
70               
71            short.levels <- wd(data.short[, 1], filter.number = filter, bc = bc)\$nlevels;
72            title <- c(\"motif\");
73            for (i in 1:short.levels){
74                title <- c(title, paste(i, \"moment2\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"), paste(i, \"test\", sep = \"_\"));
75            }
76            print(title);
77       
78            # loop to compare a vs a
79            for(i in 1:length(names.short)){
80                        wave1.dwt = NULL;
81                        m2.dwt = diff = var.dwt = NULL;
82                        out = NULL;
83                out <- vector(length = length(title));
84       
85                        print(names.short[i]);
86                        print(names.long[i]);
87                       
88                        # need exit if not comparing motif(a) vs motif(a)
89                        if (names.short[i] != names.long[i]){
90                        stop(paste(\"motif\", names.short[i], \"is not the same as\", names.long[i], sep = \" \"));
91                        }
92                        else {
93                        # signal is the difference I-C data sets
94                    diff<-data.short[,i]-data.long[,i];
95       
96                    # normalize the signal
97                    diff<-norm(diff);
98       
99                    # function is 2nd moment
100                    # 2nd moment m_j = 1/N[sum_N(W_j + V_J)^2] = 1/N sum_N(W_j)^2 + (X_bar)^2
101                        wave1.dwt <- dwt(diff, wf = wf, short.levels, boundary = boundary);
102                        var.dwt <- wave.variance(wave1.dwt);
103                        m2.dwt <- vector(length = short.levels)
104                    for(level in 1:short.levels){
105                        m2.dwt[level] <- var.dwt[level, 1] + (mean(diff)^2);
106                    }
107                               
108                        # CI bands by permutation of time series
109                        feature1 = feature2 = NULL;
110                        feature1 = data.short[, i];
111                        feature2 = data.long[, i];
112                        null = results = med = NULL;
113                        m2_25 = m2_975 = NULL;
114           
115                        for (k in 1:1000) {
116                                nk_1 = nk_2 = NULL;
117                                m2_null = var_null = NULL;
118                                null.levels = null_wave1 = null_diff = NULL;
119                                nk_1 <- sample(feature1, length(feature1), replace = FALSE);
120                                nk_2 <- sample(feature2, length(feature2), replace = FALSE);
121                                null.levels <- wd(nk_1, filter.number = filter, bc = bc)\$nlevels;
122                                null_diff <- nk_1-nk_2;
123                                null_diff <- norm(null_diff);
124                                null_wave1 <- dwt(null_diff, wf = wf, short.levels, boundary = boundary);
125                        var_null <- wave.variance(null_wave1);
126                                m2_null <- vector(length = null.levels);
127                                for(level in 1:null.levels){
128                                m2_null[level] <- var_null[level, 1] + (mean(null_diff)^2);
129                                }
130                                null= rbind(null, m2_null);
131                        }
132               
133                        null <- apply(null, 2, sort, na.last = TRUE);
134                        m2_25 <- null[25,];
135                        m2_975 <- null[975,];
136                        med <- apply(null, 2, median, na.rm = TRUE);
137
138                        # plot
139                        results <- cbind(m2.dwt, m2_25, m2_975);
140                        matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2), xlab = \"Wavelet Scale\", ylab = c(\"Wavelet 2nd Moment\", test), main = (names.short[i]), cex.main = 0.75);
141                        abline(h = 1);
142
143                        # get pvalues by comparison to null distribution
144                        out <- c(names.short[i]);
145                        for (m in 1:length(m2.dwt)){
146                        print(paste(\"scale\", m, sep = \" \"));
147                        print(paste(\"m2\", m2.dwt[m], sep = \" \"));
148                        print(paste(\"median\", med[m], sep = \" \"));
149                        out <- c(out, format(m2.dwt[m], digits = 4));   
150                        pv = NULL;
151                        if(is.na(m2.dwt[m])){
152                                pv <- \"NA\";
153                        }
154                        else {
155                                if (m2.dwt[m] >= med[m]){
156                                # R tail test
157                                tail <- \"R\";
158                                pv <- (length(which(null[, m] >= m2.dwt[m])))/(length(na.exclude(null[, m])));
159                            }
160                            else{
161                                if (m2.dwt[m] < med[m]){
162                                        # L tail test
163                                    tail <- \"L\";
164                                    pv <- (length(which(null[, m] <= m2.dwt[m])))/(length(na.exclude(null[, m])));
165                                }
166                            }
167                        }
168                        out <- c(out, pv);
169                        print(pv); 
170                        out <- c(out, tail);
171                    }
172                    final_pvalue <-rbind(final_pvalue, out);
173                    print(out);
174                }
175            }
176           
177            colnames(final_pvalue) <- title;
178            write.table(final_pvalue, file = table, sep = \"\\t\", quote = FALSE, row.names = FALSE);
179            dev.off();
180        }\n";
181
182print Rcmd "
183        # execute
184        # read in data
185       
186        inputData <- read.delim(\"$firstInputFile\");
187        inputDataNames <- colnames(inputData);
188       
189        controlData <- read.delim(\"$secondInputFile\");
190        controlDataNames <- colnames(controlData);
191       
192        # call the test function to implement IvC test
193        dwt_cor(inputData, inputDataNames, controlData, controlDataNames, test = \"$test\", pdf = \"$secondOutputFile\", table = \"$firstOutputFile\");
194        print (\"done with the correlation test\");
195\n";
196
197print Rcmd "#eof\n";
198
199close Rcmd;
200
201system("echo \"wavelet IvC test started on \`hostname\` at \`date\`\"\n");
202system("R --no-restore --no-save --no-readline < $r_script > $r_script.out\n");
203system("echo \"wavelet IvC test ended on \`hostname\` at \`date\`\"\n");
204
205#close the input and output and error files
206close(ERROR);
207close(OUTPUT2);
208close(OUTPUT1);
209close(INPUT2);
210close(INPUT1);
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。