|
0
|
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())
|
