Mercurial > repos > workflow4metabolomics > mixmodel4repeated_measures

comparison mixmodel_wrapper.R @ 0:a4d89d47646f 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/tools-metabolomics commit 8d2ca678d973501b60479a8dc3f212eecd56eab8
author workflow4metabolomics
date Mon, 16 May 2022 09:25:01 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a4d89d47646f
1 #!/usr/bin/env Rscript
2
3 library(lme4) ## mixed model computing
4 library(Matrix)
5 library(MASS)
6 library(lmerTest) ## computing pvalue and lsmeans from results of lme4 package
7 library(multtest) ## multiple testing
8 library(ggplot2)
9 library(gridExtra)
10 library(grid)
11
12 source_local <- function(fname) {
13 argv <- commandArgs(trailingOnly = FALSE)
14 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
15 source(paste(base_dir, fname, sep = "/"))
16 }
17
18 source_local("mixmodel_script.R")
19 source_local("diagmfl.R")
20
21 parse_args <- function() {
22 args <- commandArgs()
23 start <- which(args == "--args")[1] + 1
24 if (is.na(start)) {
25 return(list())
26 }
27 seq_by2 <- seq(start, length(args), by = 2)
28 result <- as.list(args[seq_by2 + 1])
29 names(result) <- args[seq_by2]
30 return(result)
31 }
32
33 argVc <- unlist(parse_args())
34
35 ##------------------------------
36 ## Initializing
37 ##------------------------------
38
39 ## options
40 ##--------
41
42 strAsFacL <- options()$stringsAsFactors
43 options(stringsAsFactors = FALSE)
44
45 ## constants
46 ##----------
47
48 modNamC <- "mixmodel" ## module name
49
50 topEnvC <- environment()
51 flagC <- "\n"
52
53 ## functions
54 ##----------
55
56 flgF <- function(tesC,
57 envC = topEnvC,
58 txtC = NA) { ## management of warning and error messages
59
60 tesL <- eval(parse(text = tesC), envir = envC)
61
62 if (!tesL) {
63
64 sink(NULL)
65 stpTxtC <- ifelse(is.na(txtC),
66 paste0(tesC, " is FALSE"),
67 txtC)
68
69 stop(stpTxtC,
70 call. = FALSE)
71
72 }
73
74 } ## flgF
75
76 ## log file
77 ##---------
78
79 sink(argVc["information"])
80
81 cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "")
82 cat("\nParameters used:\n\n")
83 print(argVc)
84 cat("\n\n")
85
86 ## loading
87 ##--------
88
89 datMN <- t(as.matrix(read.table(argVc["dataMatrix_in"],
90 check.names = FALSE,
91 header = TRUE,
92 row.names = 1,
93 sep = "\t")))
94
95 samDF <- read.table(argVc["sampleMetadata_in"],
96 check.names = FALSE,
97 header = TRUE,
98 row.names = 1,
99 sep = "\t")
100
101 varDF <- read.table(argVc["variableMetadata_in"],
102 check.names = FALSE,
103 header = TRUE,
104 row.names = 1,
105 sep = "\t")
106
107
108 ## checking
109 ##---------
110
111 flgF("identical(rownames(datMN), rownames(samDF))", txtC = "Column names of the dataMatrix are not identical to the row names of the sampleMetadata; check your data with the 'Check Format' module in the 'Quality Control' section")
112 flgF("identical(colnames(datMN), rownames(varDF))", txtC = "Row names of the dataMatrix are not identical to the row names of the variableMetadata; check your data with the 'Check Format' module in the 'Quality Control' section")
113
114 flgF("argVc['time'] %in% colnames(samDF)", txtC = paste0("Required time factor '", argVc["time"], "' could not be found in the column names of the sampleMetadata"))
115 flgF("argVc['subject'] %in% colnames(samDF)", txtC = paste0("Required subject factor '", argVc["subject"], "' could not be found in the column names of the sampleMetadata"))
116
117 flgF("mode(samDF[, argVc['time']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["time"], "' column of the sampleMetadata should contain either number only, or character only"))
118 flgF("mode(samDF[, argVc['subject']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["subject"], "' column of the sampleMetadata should contain either number only, or character only"))
119
120 flgF("argVc['adjC'] %in% c('holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none')")
121 flgF("argVc['trf'] %in% c('none', 'log10', 'log2')")
122
123 flgF("0 <= as.numeric(argVc['thrN']) && as.numeric(argVc['thrN']) <= 1", txtC = "(corrected) p-value threshold must be between 0 and 1")
124 flgF("argVc['diaR'] %in% c('no', 'yes')")
125
126
127 ##------------------------------
128 ## Formating
129 ##------------------------------
130
131 if (argVc["dff"] == "Satt") {
132 dffmeth <- "Satterthwaite"
133 } else {
134 dffmeth <- "Kenward-Roger"
135 }
136
137
138 ##------------------------------
139 ## Computation
140 ##------------------------------
141
142
143 varDFout <- lmixedm(datMN = datMN,
144 samDF = samDF,
145 varDF = varDF,
146 fixfact = argVc["fixfact"],
147 time = argVc["time"],
148 subject = argVc["subject"],
149 logtr = argVc["trf"],
150 pvalCutof = argVc["thrN"],
151 pvalcorMeth = argVc["adjC"],
152 dffOption = dffmeth,
153 visu = argVc["diaR"],
154 least.confounded = FALSE,
155 outlier.limit = 3,
156 pdfC = argVc["out_graph_pdf"],
157 pdfE = argVc["out_estim_pdf"]
158 )
159
160
161
162
163 ##------------------------------
164 ## Rounding
165 ##------------------------------
166
167 if (argVc["rounding"] == "yes") {
168 varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))] <- apply(varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))], 2, round, digits = as.numeric(argVc["decplaces"]))
169 }
170
171 ##------------------------------
172 ## Ending
173 ##------------------------------
174
175
176 ## saving
177 ##--------
178
179 varDFout <- cbind.data.frame(variableMetadata = rownames(varDFout),
180 varDFout)
181
182 write.table(varDFout,
183 file = argVc["variableMetadata_out"],
184 quote = FALSE,
185 row.names = FALSE,
186 sep = "\t")
187
188 ## closing
189 ##--------
190
191 cat("\n\nEnd of '", modNamC, "' Galaxy module call: ",
192 as.character(Sys.time()), "\n", sep = "")
193 cat("\nInformation about R (version, Operating System, attached or loaded packages):\n\n")
194 sessionInfo()
195
196 sink()
197
198 options(stringsAsFactors = strAsFacL)
199
200 rm(list = ls())