Mercurial > repos > jfrancoismartin > mixmodel4repeated_measures

comparison mixmodel_wrapper.R @ 0:1422de181204 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/mixmodel4repeated_measures commit 6ea32b3182383c19e5333201d2385a61d8da3d50
author jfrancoismartin
date Wed, 10 Oct 2018 05:18:42 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1422de181204
1 #!/usr/bin/env Rscript
2
3 library(batch) ## parseCommandArgs
4
5 library(lme4) ## mixed model computing
6 library(Matrix)
7 library(lmerTest) ## computing pvalue and lsmeans from results of lme4 package
8 library(multtest) ## multiple testing
9
10
11 source_local <- function(fname){
12 argv <- commandArgs(trailingOnly = FALSE)
13 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
14 source(paste(base_dir, fname, sep="/"))
15 }
16
17 #source_local("univariate_script.R")
18 source_local("mixmodel_script.R")
19 source_local("diagmfl.R")
20
21 argVc <- unlist(parseCommandArgs(evaluate=FALSE))
22
23 ##------------------------------
24 ## Initializing
25 ##------------------------------
26
27 ## options
28 ##--------
29
30 strAsFacL <- options()$stringsAsFactors
31 options(stringsAsFactors = FALSE)
32
33 ## constants
34 ##----------
35
36 modNamC <- "mixmodel" ## module name
37
38 topEnvC <- environment()
39 flagC <- "\n"
40
41 ## functions
42 ##----------
43
44 flgF <- function(tesC,
45 envC = topEnvC,
46 txtC = NA) { ## management of warning and error messages
47
48 tesL <- eval(parse(text = tesC), envir = envC)
49
50 if(!tesL) {
51
52 sink(NULL)
53 stpTxtC <- ifelse(is.na(txtC),
54 paste0(tesC, " is FALSE"),
55 txtC)
56
57 stop(stpTxtC,
58 call. = FALSE)
59
60 }
61
62 } ## flgF
63
64 ## log file
65 ##---------
66
67 sink(argVc["information"])
68
69 cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep="")
70
71 ## loading
72 ##--------
73
74 datMN <- t(as.matrix(read.table(argVc["dataMatrix_in"],
75 check.names = FALSE,
76 header = TRUE,
77 row.names = 1,
78 sep = "\t")))
79
80 samDF <- read.table(argVc["sampleMetadata_in"],
81 check.names = FALSE,
82 header = TRUE,
83 row.names = 1,
84 sep = "\t")
85
86 varDF <- read.table(argVc["variableMetadata_in"],
87 check.names = FALSE,
88 header = TRUE,
89 row.names = 1,
90 sep = "\t")
91
92
93 ## checking
94 ##---------
95
96 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")
97 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")
98
99 flgF("argVc['fixfact'] %in% colnames(samDF)", txtC = paste0("Required fixed factor '" , argVc['fixfact'], "' could not be found in the column names of the sampleMetadata"))
100 flgF("argVc['time'] %in% colnames(samDF)", txtC = paste0("Required time factor '" , argVc['time'] , "' could not be found in the column names of the sampleMetadata"))
101 flgF("argVc['subject'] %in% colnames(samDF)", txtC = paste0("Required subject factor '", argVc['subject'], "' could not be found in the column names of the sampleMetadata"))
102
103 flgF("mode(samDF[, argVc['fixfact']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc['fixfact'], "' column of the sampleMetadata should contain either number only, or character only"))
104 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"))
105 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"))
106
107 flgF("argVc['adjC'] %in% c('holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none')")
108 flgF("argVc['trf'] %in% c('none', 'log10', 'log2')")
109
110 flgF("0 <= as.numeric(argVc['thrN']) && as.numeric(argVc['thrN']) <= 1", txtC = "(corrected) p-value threshold must be between 0 and 1")
111 flgF("argVc['diaR'] %in% c('no', 'yes')")
112
113
114 ##------------------------------
115 ## Computation
116 ##------------------------------
117
118
119 varDF <- lmixedm(datMN = datMN,
120 samDF = samDF,
121 varDF = varDF,
122 fixfact = argVc["fixfact"],
123 time = argVc["time"],
124 subject = argVc["subject"],
125 logtr = argVc['trf'],
126 pvalCutof = argVc['thrN'],
127 pvalcorMeth= argVc["adjC"],
128 dffOption = "Satterthwaite",
129 visu = argVc["diaR"],
130 least.confounded = FALSE,
131 outlier.limit = 3,
132 pdfC = argVc["out_graph_pdf"],
133 pdfE = argVc["out_estim_pdf"]
134 )
135
136
137 ##------------------------------
138 ## Ending
139 ##------------------------------
140
141
142 ## saving
143 ##--------
144
145 varDF <- cbind.data.frame(variableMetadata = rownames(varDF),
146 varDF)
147
148 write.table(varDF,
149 file = argVc["variableMetadata_out"],
150 quote = FALSE,
151 row.names = FALSE,
152 sep = "\t")
153
154 ## closing
155 ##--------
156
157 cat("\nEnd of '", modNamC, "' Galaxy module call: ",
158 as.character(Sys.time()), "\n", sep = "")
159
160 sink()
161
162 options(stringsAsFactors = strAsFacL)
163
164 rm(list = ls())