Mercurial > repos > devteam > dwt_var_perclass

comparison execute_dwt_var_perClass.R @ 1:781e68074f84 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-devteam/tree/master/tools/dwt_var_perclass commit f929353ffb0623f2218d7dec459c7da62f3b0d24"
author devteam
date Mon, 06 Jul 2020 20:34:10 -0400
parents
children
comparison
equal deleted inserted replaced
0:cb422b6f49d2 1:781e68074f84
1 ######################################################################
2 ## plot power spectra, i.e. wavelet variance by class
3 ## add code to create null bands by permuting the original data series
4 ## get class of maximum significant variance per feature
5 ## generate plots and table matrix of variance including p-values
6 ######################################################################
7 library("wavethresh");
8 library("waveslim");
9
10 options(echo = FALSE)
11
12 ## normalize data
13 norm <- function(data) {
14 v <- (data - mean(data)) / sd(data);
15 if (sum(is.na(v)) >= 1) {
16 v <- data;
17 }
18 return(v);
19 }
20
21 dwt_var_permut_get_max <- function(data, names, outfile, filter = 4, bc = "symmetric", method = "kendall", wf = "haar", boundary = "reflection") {
22 max_var <- NULL;
23 matrix <- NULL;
24 title <- NULL;
25 final_pvalue <- NULL;
26 short_levels <- NULL;
27 scale <- NULL;
28
29 print(names);
30
31 par(mfcol = c(length(names), length(names)), mar = c(0, 0, 0, 0), oma = c(4, 3, 3, 2), xaxt = "s", cex = 1, las = 1);
32
33 short_levels <- wavethresh::wd(data[, 1], filter.number = filter, bc = bc)$nlevels;
34
35 title <- c("motif");
36 for (i in seq_len(short_levels)) {
37 title <- c(title, paste(i, "var", sep = "_"), paste(i, "pval", sep = "_"), paste(i, "test", sep = "_"));
38 }
39 print(title);
40
41 ## normalize the raw data
42 data <- apply(data, 2, norm);
43
44 for (i in seq_len(length(names))) {
45 for (j in seq_len(length(names))) {
46 temp <- NULL;
47 results <- NULL;
48 wave1_dwt <- NULL;
49 out <- NULL;
50
51 out <- vector(length = length(title));
52 temp <- vector(length = short_levels);
53
54 if (i != j) {
55 plot(temp, type = "n", axes = FALSE, xlab = NA, ylab = NA);
56 box(col = "grey");
57 grid(ny = 0, nx = NULL);
58 } else {
59
60 wave1_dwt <- waveslim::dwt(data[, i], wf = wf, short_levels, boundary = boundary);
61
62 temp_row <- (short_levels + 1) * -1;
63 temp_col <- 1;
64 temp <- waveslim::wave.variance(wave1_dwt)[temp_row, temp_col];
65
66 ##permutations code :
67 feature1 <- NULL;
68 null <- NULL;
69 var_25 <- NULL;
70 var_975 <- NULL;
71 med <- NULL;
72
73 feature1 <- data[, i];
74 for (k in seq_len(1000)) {
75 nk_1 <- NULL;
76 null_levels <- NULL;
77 var <- NULL;
78 null_wave1 <- NULL;
79
80 nk_1 <- sample(feature1, length(feature1), replace = FALSE);
81 null_levels <- wavethresh::wd(nk_1, filter.number = filter, bc = bc)$nlevels;
82 var <- vector(length = length(null_levels));
83 null_wave1 <- waveslim::dwt(nk_1, wf = wf, short_levels, boundary = boundary);
84 var <- waveslim::wave.variance(null_wave1)[-8, 1];
85 null <- rbind(null, var);
86 }
87 null <- apply(null, 2, sort, na.last = TRUE);
88 var_25 <- null[25, ];
89 var_975 <- null[975, ];
90 med <- (apply(null, 2, median, na.rm = TRUE));
91
92 ## plot
93 results <- cbind(temp, var_25, var_975);
94 matplot(results, type = "b", pch = "*", lty = 1, col = c(1, 2, 2), axes = F);
95
96 ## get pvalues by comparison to null distribution
97 out <- (names[i]);
98 for (m in seq_len(length(temp))) {
99 print(paste("scale", m, sep = " "));
100 print(paste("var", temp[m], sep = " "));
101 print(paste("med", med[m], sep = " "));
102 pv <- NULL;
103 tail <- NULL;
104 out <- c(out, format(temp[m], digits = 3));
105 if (temp[m] >= med[m]) {
106 ## R tail test
107 print("R");
108 tail <- "R";
109 pv <- (length(which(null[, m] >= temp[m]))) / (length(na.exclude(null[, m])));
110
111 } else {
112 ## L tail test
113 print("L");
114 tail <- "L";
115 pv <- (length(which(null[, m] <= temp[m]))) / (length(na.exclude(null[, m])));
116 }
117 out <- c(out, pv);
118 print(pv);
119 out <- c(out, tail);
120 ## get variances outside null bands by comparing temp to null
121 ### temp stores variance for each scale, and null stores permuted variances for null bands
122 if (temp[m] <= var_975[m]) {
123 temp[m] <- NA;
124 }
125 }
126 final_pvalue <- rbind(final_pvalue, out);
127 matrix <- rbind(matrix, temp)
128 }
129 ## labels
130 if (i == 1) {
131 mtext(names[j], side = 2, line = 0.5, las = 3, cex = 0.25);
132 }
133 if (j == 1) {
134 mtext(names[i], side = 3, line = 0.5, cex = 0.25);
135 }
136 if (j == length(names)) {
137 axis(1, at = (1:short_levels), las = 3, cex.axis = 0.5);
138 }
139 }
140 }
141 colnames(final_pvalue) <- title;
142
143 ## get maximum variance larger than expectation by comparison to null bands
144 varnames <- vector();
145 for (i in seq_len(length(names))) {
146 name1 <- paste(names[i], "var", sep = "_")
147 varnames <- c(varnames, name1)
148 }
149 rownames(matrix) <- varnames;
150 colnames(matrix) <- (1:short_levels);
151 max_var <- names;
152 scale <- vector(length = length(names));
153 for (x in seq_len(nrow(matrix))) {
154 if (length(which.max(matrix[x, ])) == 0) {
155 scale[x] <- NA;
156 } else{
157 scale[x] <- colnames(matrix)[which.max(matrix[x, ])];
158 }
159 }
160 max_var <- cbind(max_var, scale);
161 write.table(max_var, file = outfile, sep = "\t", quote = FALSE, row.names = FALSE, append = TRUE);
162 return(final_pvalue);
163 }
164
165 ## execute
166 ## read in data
167 args <- commandArgs(trailingOnly = TRUE)
168
169 data_test <- NULL;
170 data_test <- read.delim(args[1]);
171
172 count <- ncol(data_test)
173 print(paste("The number of columns in the input file is: ", count));
174
175 # check if the number of motifs is not a multiple of 12, and round up is so
176 if (count %% 12 != 0) {
177 print("the number of motifs is not a multiple of 12")
178 count2 <- ceiling(count / 12);
179 }else{
180 print("the number of motifs is a multiple of 12")
181 count2 <- count / 12
182 }
183 print(paste("There will be", count2, "subfiles"))
184
185 pdf(file = args[4], width = 11, height = 8);
186
187 ## loop to read and execute on all count2 subfiles
188 final <- NULL;
189 for (x in seq_len(count2)) {
190 sub <- NULL;
191 sub_names <- NULL;
192 a <- NULL;
193 b <- NULL;
194
195 a <- ((x - 1) * 12 + 1);
196 b <- x * 12;
197
198 if (x < count2) {
199 sub <- data_test[, +c(a:b)];
200 sub_names <- colnames(data_test)[a:b];
201 final <- rbind(final, dwt_var_permut_get_max(sub, sub_names, args[2]));
202 }
203 else{
204 sub <- data_test[, +c(a:ncol(data_test))];
205 sub_names <- colnames(data_test)[a:ncol(data_test)];
206 final <- rbind(final, dwt_var_permut_get_max(sub, sub_names, args[2]));
207 }
208 }
209
210 dev.off();
211
212 write.table(final, file = args[3], sep = "\t", quote = FALSE, row.names = FALSE);