Mercurial > repos > marie-tremblay-metatoul > asca_w4m

comparison galaxy/asca_wrapper.R @ 0:c5f11e6f8f99 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Mon, 30 Jul 2018 07:29:40 -0400
parents
children 20395c0079ae
comparison
equal deleted inserted replaced
-1:000000000000 0:c5f11e6f8f99
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
116 samDF <- read.table(listArguments[["sampleMetadata_in"]],
117 check.names = FALSE,
118 header = TRUE,
119 row.names = 1,
120 sep = "\t")
121
122 varDF <- read.table(listArguments[["variableMetadata_in"]],
123 check.names = FALSE,
124 header = TRUE,
125 row.names = 1,
126 sep = "\t")
127
128 result <- asca_w4m(xMN, samDF, c(listArguments[["factor1"]],listArguments[["factor2"]]), varDF, as.numeric(listArguments[["threshold"]]),
129 scaling=listArguments[["scaling"]], listArguments[["nPerm"]])
130
131
132 ##saving
133
134 if (exists("result")) {
135 ## writing output files
136 cat("\n\nWriting output files\n\n");
137 write.table(result[[4]],
138 file = listArguments$sampleMetadata_out,
139 quote = FALSE,
140 row.names = TRUE,
141 sep = "\t")
142
143 write.table(result[[5]],
144 file = listArguments$variableMetadata_out,
145 quote = FALSE,
146 row.names = TRUE,
147 sep = "\t")
148
149 # Graphical display for each significant parameter
150 print(result[[3]])
151 cat("\n p-value of Residuals must not be taken into account\n")
152
153 if (any(result[[2]] < as.numeric(listArguments[["threshold"]])))
154 {
155 data.asca.permutation <- result[[2]]
156 design <- data.matrix(samDF[, colnames(samDF) %in% c(listArguments[["factor1"]],listArguments[["factor2"]])])
157
158 pdf(listArguments$figure, onefile=TRUE)
159 par(mfrow=c(1,3))
160 if (data.asca.permutation[1] < as.numeric(listArguments[["threshold"]]))
161 {
162 eigenvalues <- data.frame(1:length(unique(design[,1])), result[[1]]$'1'$svd$var.explained[1:length(unique(design[,1]))])
163 colnames(eigenvalues) <- c("PC", "explainedVariance")
164 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
165 noms <- levels(as.factor(samDF[, listArguments$factor1]))
166 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="1", interaction=0, factorName=listArguments$factor1, factorModalite=noms)
167 Date.loadings <- data.matrix(result[[5]][,2:3])
168 Date.loadings.leverage <- diag(Date.loadings%*%t(Date.loadings))
169 names(Date.loadings.leverage) <- colnames(xMN)
170 Date.loadings.leverage <- sort(Date.loadings.leverage, decreasing=TRUE)
171 barplot(Date.loadings.leverage[Date.loadings.leverage > 0.001], main="PC1 loadings")
172 }
173 if (data.asca.permutation[2] < as.numeric(listArguments[["threshold"]]))
174 {
175 eigenvalues <- data.frame(1:length(unique(design[,2])), result[[1]]$'2'$svd$var.explained[1:length(unique(design[,2]))])
176 colnames(eigenvalues) <- c("PC", "explainedVariance")
177 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
178 noms <- levels(as.factor(samDF[, listArguments$factor2]))
179 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="2", interaction=0, factorName=listArguments$factor2, factorModalite=noms)
180 Date.loadings <- data.matrix(result[[5]][,4:5])
181 Date.loadings.leverage <- diag(Date.loadings%*%t(Date.loadings))
182 names(Date.loadings.leverage) <- colnames(xMN)
183 Date.loadings.leverage <- sort(Date.loadings.leverage, decreasing=TRUE)
184 barplot(Date.loadings.leverage[Date.loadings.leverage > 0.001], main="PC1 loadings")
185 }
186 if (data.asca.permutation[3] < as.numeric(listArguments[["threshold"]]))
187 {
188 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])))])
189 colnames(eigenvalues) <- c("PC", "explainedVariance")
190 barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component")
191 noms1 <- data.matrix(levels(as.factor(samDF[, listArguments$factor1])))
192 noms2 <- data.matrix(levels(as.factor(samDF[, listArguments$factor2])))
193 noms <- apply(noms1, 1, FUN=function(x){paste(x, "-", noms2, sep="")})
194 noms <- apply(noms, 1, FUN=function(x){c(noms)})
195 ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="12", interaction=1, factorModalite=noms[,1])
196 Date.loadings <- data.matrix(result[[5]][,6:7])
197 Date.loadings.leverage <- diag(Date.loadings%*%t(Date.loadings))
198 names(Date.loadings.leverage) <- colnames(xMN)
199 Date.loadings.leverage <- sort(Date.loadings.leverage, decreasing=TRUE)
200 barplot(Date.loadings.leverage[Date.loadings.leverage > 0.001], main="PC1 loadings")
201 }
202 dev.off()
203 }
204
205 tryCatch({
206 save(result, file="asca.RData");
207 }, warning = function(w) {
208 print(paste("Warning: ", w));
209 }, error = function(err) {
210 stop(paste("ERROR saving result RData object:", err));
211 });
212 }
213
214 ## ending
215 ##-------
216
217 cat("\nEnd of the '", modNamC, "' Galaxy module call: ",
218 format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "")
219
220 sink()
221
222 # options(stringsAsFactors = strAsFacL)
223
224
225 rm(list = ls())