Mercurial > repos > devteam > dwt_ivc_all

comparison execute_dwt_IvC_all.R @ 1:506ae7b0d85d 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_ivc_all commit f929353ffb0623f2218d7dec459c7da62f3b0d24"
author devteam
date Mon, 06 Jul 2020 20:31:56 -0400
parents
children
comparison
equal deleted inserted replaced
0:0b89b03ad760 1:506ae7b0d85d
1 ###########################################################################################
2 ## code to do wavelet Indel vs. Control
3 ## signal is the difference I-C; function is second moment i.e. variance from zero not mean
4 ## to perform wavelet transf. of signal, scale-by-scale analysis of the function
5 ## create null bands by permuting the original data series
6 ## generate plots and table matrix of correlation coefficients including p-values
7 ############################################################################################
8 library("wavethresh");
9 library("waveslim");
10
11 options(echo = FALSE)
12
13 ## normalize data
14 norm <- function(data) {
15 v <- (data - mean(data)) / sd(data);
16 if (sum(is.na(v)) >= 1) {
17 v <- data;
18 }
19 return(v);
20 }
21
22 dwt_cor <- function(data_short, names_short, data_long, names_long, test, pdf, table, filter = 4, bc = "symmetric", wf = "haar", boundary = "reflection") {
23 print(test);
24 print(pdf);
25 print(table);
26
27 pdf(file = pdf);
28 final_pvalue <- NULL;
29 title <- NULL;
30
31 short_levels <- wavethresh::wd(data_short[, 1], filter.number = filter, bc = bc)$nlevels;
32 title <- c("motif");
33 for (i in 1:short_levels) {
34 title <- c(title, paste(i, "moment2", sep = "_"), paste(i, "pval", sep = "_"), paste(i, "test", sep = "_"));
35 }
36 print(title);
37
38 ## loop to compare a vs a
39 for (i in seq_len(length(names_short))) {
40 wave1_dwt <- NULL;
41 m2_dwt <- NULL;
42 diff <- NULL;
43 var_dwt <- NULL;
44 out <- NULL;
45 out <- vector(length = length(title));
46
47 print(names_short[i]);
48 print(names_long[i]);
49
50 ## need exit if not comparing motif(a) vs motif(a)
51 if (names_short[i] != names_long[i]) {
52 stop(paste("motif", names_short[i], "is not the same as", names_long[i], sep = " "));
53 }
54 else {
55 ## signal is the difference I-C data sets
56 diff <- data_short[, i] - data_long[, i];
57
58 ## normalize the signal
59 diff <- norm(diff);
60
61 ## function is 2nd moment
62 ## 2nd moment m_j = 1/N[sum_N(W_j + V_J)^2] = 1/N sum_N(W_j)^2 + (X_bar)^2
63 wave1_dwt <- waveslim::dwt(diff, wf = wf, short_levels, boundary = boundary);
64 var_dwt <- waveslim::wave.variance(wave1_dwt);
65 m2_dwt <- vector(length = short_levels)
66 for (level in 1:short_levels) {
67 m2_dwt[level] <- var_dwt[level, 1] + (mean(diff)^2);
68 }
69
70 ## CI bands by permutation of time series
71 feature1 <- NULL;
72 feature2 <- NULL;
73 feature1 <- data_short[, i];
74 feature2 <- data_long[, i];
75 null <- NULL;
76 results <- NULL;
77 med <- NULL;
78 m2_25 <- NULL;
79 m2_975 <- NULL;
80
81 for (k in 1:1000) {
82 nk_1 <- NULL;
83 nk_2 <- NULL;
84 m2_null <- NULL;
85 var_null <- NULL;
86 null_levels <- NULL;
87 null_wave1 <- NULL;
88 null_diff <- NULL;
89 nk_1 <- sample(feature1, length(feature1), replace = FALSE);
90 nk_2 <- sample(feature2, length(feature2), replace = FALSE);
91 null_levels <- wavethresh::wd(nk_1, filter.number = filter, bc = bc)$nlevels;
92 null_diff <- nk_1 - nk_2;
93 null_diff <- norm(null_diff);
94 null_wave1 <- waveslim::dwt(null_diff, wf = wf, short_levels, boundary = boundary);
95 var_null <- waveslim::wave.variance(null_wave1);
96 m2_null <- vector(length = null_levels);
97 for (level in 1:null_levels) {
98 m2_null[level] <- var_null[level, 1] + (mean(null_diff)^2);
99 }
100 null <- rbind(null, m2_null);
101 }
102
103 null <- apply(null, 2, sort, na.last = TRUE);
104 m2_25 <- null[25, ];
105 m2_975 <- null[975, ];
106 med <- apply(null, 2, median, na.rm = TRUE);
107
108 ## plot
109 results <- cbind(m2_dwt, m2_25, m2_975);
110 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);
111 abline(h = 1);
112
113 ## get pvalues by comparison to null distribution
114 out <- c(names_short[i]);
115 for (m in seq_len(length(m2_dwt))) {
116 print(paste("scale", m, sep = " "));
117 print(paste("m2", m2_dwt[m], sep = " "));
118 print(paste("median", med[m], sep = " "));
119 out <- c(out, format(m2_dwt[m], digits = 4));
120 pv <- NULL;
121 if (is.na(m2_dwt[m])) {
122 pv <- "NA";
123 }
124 else {
125 if (m2_dwt[m] >= med[m]) {
126 ## R tail test
127 tail <- "R";
128 pv <- (length(which(null[, m] >= m2_dwt[m]))) / (length(na.exclude(null[, m])));
129 }
130 else{
131 if (m2_dwt[m] < med[m]) {
132 ## L tail test
133 tail <- "L";
134 pv <- (length(which(null[, m] <= m2_dwt[m]))) / (length(na.exclude(null[, m])));
135 }
136 }
137 }
138 out <- c(out, pv);
139 print(pv);
140 out <- c(out, tail);
141 }
142 final_pvalue <- rbind(final_pvalue, out);
143 print(out);
144 }
145 }
146
147 colnames(final_pvalue) <- title;
148 write.table(final_pvalue, file = table, sep = "\t", quote = FALSE, row.names = FALSE);
149 dev.off();
150 }
151 ## execute
152 ## read in data
153 args <- commandArgs(trailingOnly = TRUE)
154
155 input_data <- read.delim(args[1]);
156 input_data_names <- colnames(input_data);
157
158 control_data <- read.delim(args[2]);
159 control_data_names <- colnames(control_data);
160
161 ## call the test function to implement IvC test
162 dwt_cor(input_data, input_data_names, control_data, control_data_names, test = "IvC", pdf = args[3], table = args[4]);
163 print("done with the correlation test");