Mercurial > repos > marie-tremblay-metatoul > asca

comparison asca_wrapper.R @ 0:93312041f1d5 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/workflow4metabolomics/ascaw4m commit 7ea9b0f8abc5a60c2c04fd2098788497f14766b6
author marie-tremblay-metatoul
date Fri, 21 Sep 2018 05:51:14 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:93312041f1d5
1 #!/usr/bin/env Rscript
2
3 ###################################################################################################
4 #
5 # MetStaT ASCA.calculate function
6 #
7 #
8 # R-Package: MetStaT
9 #
10 # Version: 1.0
11 #
12 # Author (asca.calculate): Tim Dorscheidt
13 # Author (wrapper & .r adaptation for workflow4metabolomics.org): M. Tremblay-Franco & Y. Guitton #
14 #
15 # Expected parameters from the commandline
16 # input files:
17 # dataMatrix
18 # sampleMetadata
19 # variableMetadata
20 # params:
21 # Factors (Factor1 & Factor2)
22 # scaling
23 # Number of permutations
24 # Significance threshold
25 # output files:
26 # sampleMetadata
27 # variableMetadata
28 # Graphical outputs
29 # Information text
30 ###################################################################################################
31 pkgs=c("MetStaT","batch","pcaMethods")
32 for(pkg in pkgs) {
33 suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
34 cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="")
35 }
36
37
38 listArguments = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects
39
40 #Redirect all stdout to the log file
41 sink(listArguments$information)
42
43 # ----- PACKAGE -----
44 cat("\tPACKAGE INFO\n")
45 sessionInfo()
46
47 source_local <- function(fname) {
48 argv <- commandArgs(trailingOnly = FALSE)
49 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
50 source(paste(base_dir, fname, sep="/"))
51 }
52
53 #load asca_w4m function
54 source_local("asca_w4m.R")
55 source_local("ASCA.Calculate_w4m.R")
56 source_local("ASCA.PlotScoresPerLevel_w4m.R")
57 print("first loadings OK")
58
59 ## libraries
60 ##----------
61
62 cat('\n\nRunning ASCA.calculate\n');
63 options(warn=-1);
64 #remove rgl warning
65 options(rgl.useNULL = TRUE);
66
67
68 ## constants
69 ##----------
70
71 modNamC <- "asca" ## module name
72
73 topEnvC <- environment()
74 flgC <- "\n"
75
76 ## functions
77 ##----------For manual input of function
78 ##--end function
79
80 flgF <- function(tesC,
81 envC = topEnvC,
82 txtC = NA) { ## management of warning and error messages
83
84 tesL <- eval(parse(text = tesC), envir = envC)
85
86 if(!tesL) {
87
88 sink(NULL)
89 stpTxtC <- ifelse(is.na(txtC),
90 paste0(tesC, " is FALSE"),
91 txtC)
92
93 stop(stpTxtC,
94 call. = FALSE)
95
96 }
97
98 } ## flgF
99
100
101 ## log file
102 ##---------
103 cat("\nStart of the '", modNamC, "' Galaxy module call: ",
104 format(Sys.time(), "%a %d %b %Y %X"), "\n", sep="")
105
106
107 ## arguments
108 ##----------
109 ## loading files and checks
110 xMN <- t(as.matrix(read.table(listArguments[["dataMatrix_in"]],
111 check.names = FALSE,
112 header = TRUE,
113 row.names = 1,
114 sep = "\t")))
115 varIdDM <- rownames(xMN)
116
117 samDF <- read.table(listArguments[["sampleMetadata_in"]],
118 check.names = FALSE,
119 header = TRUE,
120 row.names = 1,
121 sep = "\t")
122 obsIdSMD <- rownames(samDF)
123
124 varDF <- read.table(listArguments[["variableMetadata_in"]],
125 check.names = FALSE,
126 header = TRUE,
127 row.names = 1,
128 sep = "\t")
129 varIdVDM <- rownames(varDF)
130
131 result <- asca_w4m(xMN, samDF, c(listArguments[["factor1"]],listArguments[["factor2"]]), varDF, as.numeric(listArguments[["threshold"]]),
132 scaling=listArguments[["scaling"]], listArguments[["nPerm"]])
133
134
135 ##saving
136
137 if (exists("result")) {
138 ## writing output files
139 cat("\n\nWriting output files\n\n");
140 write.table(data.frame(cbind(obsIdSMD, result[[4]])),
141 file = listArguments$sampleMetadata_out,
142 quote = FALSE,
143 row.names = FALSE,
144 sep = "\t")
145
146 write.table(data.frame(cbind(varIdVDM, result[[5]])),
147 file = listArguments$variableMetadata_out,
148 quote = FALSE,
149 row.names = FALSE,
150 sep = "\t")
151
152 # Graphical display for each significant parameter
153 print(result[[3]])
154 cat("\n p-value of Residuals must not be taken into account\n")
155
156 if (any(result[[2]] < as.numeric(listArguments[["threshold"]])))
157 {
158 data.asca.permutation <- result[[2]]
159 design <- data.matrix(samDF[, colnames(samDF) %in% c(listArguments[["factor1"]],listArguments[["factor2"]])])
160
161 pdf(listArguments$figure, onefile=TRUE)
162 par(mfrow=c(1,3))
163 if (data.asca.permutation[1] < as.numeric(listArguments[["threshold"]]))
164 {
165 eigenvalues <- data.frame(1:length(unique(design[,1])), result[[1]]$'1'$svd$var.explained[1:length(unique(design[,1]))])
166 colnames(eigenvalues) <- c("PC", "explainedVariance")
167 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
168 noms <- levels(as.factor(samDF[, listArguments$factor1]))
169 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="1", interaction=0, factorName=listArguments$factor1, factorModalite=noms)
170
171 v1 <- paste(listArguments[["factor1"]],"_XLOAD-p1", sep="")
172 v2 <- paste(listArguments[["factor1"]],"_XLOAD-p2", sep="")
173 f1.loadings <- data.matrix(result[[5]][,c(v1, v2)])
174 f1.loadings.leverage <- diag(f1.loadings%*%t(f1.loadings))
175 names(f1.loadings.leverage) <- colnames(xMN)
176 f1.loadings.leverage <- sort(f1.loadings.leverage, decreasing=TRUE)
177 barplot(f1.loadings.leverage[f1.loadings.leverage > 0.001], main="Leverage values")
178 }
179 if (data.asca.permutation[2] < as.numeric(listArguments[["threshold"]]))
180 {
181 eigenvalues <- data.frame(1:length(unique(design[,2])), result[[1]]$'2'$svd$var.explained[1:length(unique(design[,2]))])
182 colnames(eigenvalues) <- c("PC", "explainedVariance")
183 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
184 noms <- levels(as.factor(samDF[, listArguments$factor2]))
185 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="2", interaction=0, factorName=listArguments$factor2, factorModalite=noms)
186
187 v1 <- paste(listArguments[["factor2"]],"_XLOAD-p1", sep="")
188 v2 <- paste(listArguments[["factor2"]],"_XLOAD-p2", sep="")
189 f2.loadings <- data.matrix(result[[5]][,c(v1, v2)])
190 f2.loadings.leverage <- diag(f2.loadings%*%t(f2.loadings))
191 names(f2.loadings.leverage) <- colnames(xMN)
192 f2.loadings.leverage <- sort(f2.loadings.leverage, decreasing=TRUE)
193 barplot(f2.loadings.leverage[f2.loadings.leverage > 0.001], main="Leverage values")
194 }
195 if (data.asca.permutation[3] < as.numeric(listArguments[["threshold"]]))
196 {
197 eigenvalues <- data.frame(1:(length(unique(design[,1]))*length(unique(design[,2]))), result[[1]]$'12'$svd$var.explained[1:(length(unique(design[,1]))*length(unique(design[,2])))])
198 colnames(eigenvalues) <- c("PC", "explainedVariance")
199 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
200 noms1 <- data.matrix(levels(as.factor(samDF[, listArguments$factor1])))
201 noms2 <- data.matrix(levels(as.factor(samDF[, listArguments$factor2])))
202 noms <- apply(noms1, 1, FUN=function(x){paste(x, "-", noms2, sep="")})
203 noms <- apply(noms, 1, FUN=function(x){c(noms)})
204 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="12", interaction=1, factorModalite=noms[,1])
205
206 v1 <- "Interact_XLOAD-p1"
207 v2 <- "Interact_XLOAD-p2"
208 f1f2.loadings <- data.matrix(result[[5]][, c(v1, v2)])
209 f1f2.loadings.leverage <- diag(f1f2.loadings%*%t(f1f2.loadings))
210 names(f1f2.loadings.leverage) <- colnames(xMN)
211 f1f2.loadings.leverage <- sort(f1f2.loadings.leverage, decreasing=TRUE)
212 barplot(f1f2.loadings.leverage[f1f2.loadings.leverage > 0.001], main="Leverage values")
213 }
214 dev.off()
215 }
216
217 tryCatch({
218 save(result, file="asca.RData");
219 }, warning = function(w) {
220 print(paste("Warning: ", w));
221 }, error = function(err) {
222 stop(paste("ERROR saving result RData object:", err));
223 });
224 }
225
226 ## ending
227 ##-------
228
229 cat("\nEnd of the '", modNamC, "' Galaxy module call: ",
230 format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "")
231
232 sink()
233
234 # options(stringsAsFactors = strAsFacL)
235
236
237 rm(list = ls())