|
0
|
1
|
|
|
2
|
|
19
|
3 assign("TOOL_NAME", "XSeekerPreparator", envir = globalenv())
|
|
|
4 lockBinding("TOOL_NAME", globalenv())
|
|
|
5 assign("VERSION", "1.3.0", envir = globalenv())
|
|
|
6 lockBinding("VERSION", globalenv())
|
|
|
7 assign("DEBUG_FAST", FALSE, envir = globalenv())
|
|
|
8 lockBinding("DEBUG_FAST", globalenv())
|
|
|
9 assign("DEBUG_FAST_IGNORE_SLOW_OP", DEBUG_FAST, envir = globalenv())
|
|
|
10 lockBinding("DEBUG_FAST_IGNORE_SLOW_OP", globalenv())
|
|
|
11 assign("PROCESS_SMOL_BATCH", DEBUG_FAST, envir = globalenv())
|
|
|
12 lockBinding("PROCESS_SMOL_BATCH", globalenv())
|
|
|
13 assign("FAST_FEATURE_RATIO", 10, envir = globalenv())
|
|
|
14 lockBinding("FAST_FEATURE_RATIO", globalenv())
|
|
|
15 assign("OUTPUT_SPECIFIC_TOOL", "XSeeker_Galaxy", envir = globalenv())
|
|
|
16 lockBinding("OUTPUT_SPECIFIC_TOOL", globalenv())
|
|
14
|
17
|
|
19
|
18 assign(
|
|
|
19 "ENRICHED_RDATA_VERSION",
|
|
|
20 paste(VERSION, OUTPUT_SPECIFIC_TOOL, sep = "-"),
|
|
|
21 envir = globalenv()
|
|
|
22 )
|
|
|
23 lockBinding("ENRICHED_RDATA_VERSION", globalenv())
|
|
|
24 assign("ENRICHED_RDATA_DOC", sprintf("
|
|
0
|
25 Welcome to the enriched <Version %s> of the output of CAMERA/xcms.
|
|
|
26 This doc was generated by the tool: %s - Version %s
|
|
|
27 To show the different variables contained in this rdata, type:
|
|
|
28 - `load('this_rdata.rdata', rdata_env <- new.env())`
|
|
|
29 - `names(rdata_env)`
|
|
|
30
|
|
|
31 Sections
|
|
|
32 ######
|
|
|
33
|
|
|
34
|
|
|
35 This tools helpers
|
|
|
36 ------
|
|
|
37 The version number is somewhat special because the evolution of the
|
|
|
38 rdata's format is non-linear.
|
|
|
39 There may be different branches, each evolving separatly.
|
|
|
40 To reflect these branches's diversions, there may be a prepended
|
|
|
41 branch name following this format:
|
|
|
42 major.minor.patch-branch_name
|
|
|
43 Like this, we can process rdata with the same tool, and output
|
|
|
44 rdata formated differently, for each tool.
|
|
|
45
|
|
|
46
|
|
|
47 - enriched_rdata:
|
|
|
48 - Description: flag created by that tool to tell it was enriched.
|
|
|
49 - Retrieval method: enriched_rdata <- TRUE
|
|
|
50
|
|
|
51 - enriched_rdata_version:
|
|
|
52 - Description: A flag created by that tool to tell which version of
|
|
|
53 this tool has enriched the rdata.
|
|
19
|
54 - Retrieval method:
|
|
|
55 enriched_rdata_version <- sprintf(
|
|
|
56 \"%s\",
|
|
|
57 ENRICHED_RDATA_VERSION
|
|
|
58 )
|
|
0
|
59
|
|
|
60 - enriched_rdata_doc:
|
|
|
61 - Description: Contains the documentation string.
|
|
|
62
|
|
|
63 Data from original mzxml file
|
|
|
64 ------
|
|
|
65 - tic:
|
|
|
66 - Description: Those are the tic values from the original mzxml
|
|
|
67 file, extracted using xcms 2.
|
|
|
68 - Retrieval method: xcms::xcmsRaw('original_file.mzxml')@tic
|
|
|
69 - xcms version: 2.0
|
|
|
70
|
|
|
71 - mz:
|
|
|
72 - Description: Those are the m/z values from the original mzxml
|
|
|
73 file, extracted using xcms 2.
|
|
|
74 - Retrieval method: xcms::xcmsRaw('original_file.mzxml')@env$mz
|
|
|
75 - xcms version: 2.0
|
|
|
76
|
|
|
77 - scanindex:
|
|
|
78 - Description: Those are the scanindex values from the original mzxml
|
|
|
79 file, extracted using xcms 2.
|
|
|
80 - Retrieval method: xcms::xcmsRaw('original_file.mzxml')@scanindex
|
|
|
81 - xcms version: 2.0
|
|
|
82
|
|
|
83 - scantime:
|
|
|
84 - Description: Those are the scantime values from the original mzxml
|
|
|
85 file, extracted using xcms 2.
|
|
|
86 - Retrieval method: xcms::xcmsRaw('original_file.mzxml')@scantime
|
|
|
87 - xcms version: 2.0
|
|
|
88
|
|
|
89 - intensity:
|
|
|
90 - Description: Those are the intensity values from the original mzxml
|
|
|
91 file, extracted using xcms 2.
|
|
|
92 - Retrieval method: xcms::xcmsRaw('original_file.mzxml')@env$intensity
|
|
|
93 - xcms version: 2.0
|
|
|
94
|
|
|
95 - polarity:
|
|
|
96 - Description: Those are the polarity values from the original mzxml
|
|
|
97 file, extracted using xcms 2.
|
|
19
|
98 - Retrieval method:
|
|
|
99 as.character(xcms::xcmsRaw(
|
|
|
100 'original_file.mzxml'
|
|
|
101 )@polarity[[1]])
|
|
0
|
102 - xcms version: 2.0
|
|
|
103
|
|
|
104 Data taken from incoming rdata
|
|
|
105 ------
|
|
|
106 - variableMetadata:
|
|
|
107 - Description: Unmodified copy of variableMetadata from incoming rdata.
|
|
|
108 - Retrieval method: rdata_file$variableMetadata
|
|
|
109
|
|
|
110 - process_params:
|
|
|
111 - Description: Those are the processing parameters values from the
|
|
|
112 curent rdata. They have been simplified to allow easy access like:
|
|
|
113 for (params in process_params) {
|
|
|
114 if (params[[\"xfunction\"]] == \"annotatediff\") {
|
|
|
115 process_peak_picking_params(params)
|
|
|
116 }
|
|
|
117 }
|
|
|
118 - Retrieval method:
|
|
|
119 ## just he same list, but simplified
|
|
|
120 process_params <- list()
|
|
|
121 for (list_name in names(rdata_file$listOFlistArguments)) {
|
|
|
122 param_list <- list()
|
|
19
|
123 for (param_name in names(
|
|
|
124 rdata_file$listOFlistArguments[[list_name]]
|
|
|
125 )) {
|
|
|
126 param_list[[param_name]] <- rdata_file$listOFlistArguments[[
|
|
|
127 list_name
|
|
|
128 ]][[param_name]]
|
|
0
|
129 }
|
|
|
130 process_params[[length(process_params)+1]] <- param_list
|
|
|
131 }
|
|
19
|
132 ", ENRICHED_RDATA_VERSION, TOOL_NAME, VERSION, ENRICHED_RDATA_VERSION),
|
|
|
133 envir = globalenv())
|
|
|
134 lockBinding("ENRICHED_RDATA_DOC", globalenv())
|
|
0
|
135
|
|
|
136
|
|
|
137
|
|
|
138 get_models <- function(path) {
|
|
|
139 if (is.null(path)) {
|
|
|
140 stop("No models to define the database schema")
|
|
|
141 } else {
|
|
|
142 message(sprintf("Loading models from %s", path))
|
|
|
143 }
|
|
|
144 ## galaxy mangles the "@" to a "__at__"
|
|
|
145 if (substr(path, 1, 9) == "git__at__") {
|
|
19
|
146 path <- sub("^git__at__", "git@", path, perl = TRUE)
|
|
0
|
147 }
|
|
|
148 if (
|
|
|
149 substr(path, 1, 4) == "git@"
|
|
19
|
150 || substr(path, length(path) - 4, 4) == ".git"
|
|
0
|
151 ) {
|
|
19
|
152 return(get_models_from_git(path))
|
|
0
|
153 }
|
|
|
154 if (substr(path, 1, 4) == "http") {
|
|
19
|
155 return(get_models_from_url(path))
|
|
0
|
156 }
|
|
19
|
157 return(source(path)$value)
|
|
0
|
158 }
|
|
|
159
|
|
19
|
160 get_models_from_git <- function(url, target_file = "models.R", rm = TRUE) {
|
|
0
|
161 tmp <- tempdir()
|
|
|
162 message(sprintf("Cloning %s", url))
|
|
|
163 system2("git", c("clone", url, tmp))
|
|
|
164 result <- search_tree(file.path(tmp, dir), target_file)
|
|
|
165 if (!is.null(result)) {
|
|
|
166 models <- source(result)$value
|
|
|
167 if (rm) {
|
|
19
|
168 unlink(tmp, recursive = TRUE)
|
|
0
|
169 }
|
|
19
|
170 return(models)
|
|
0
|
171 }
|
|
|
172 if (rm) {
|
|
19
|
173 unlink(tmp, recursive = TRUE)
|
|
0
|
174 }
|
|
|
175 stop(sprintf(
|
|
|
176 "Could not find any file named \"%s\" in this repo",
|
|
|
177 target_file
|
|
|
178 ))
|
|
|
179 }
|
|
|
180
|
|
19
|
181 get_models_from_url <- function(url, target_file = "models.R", rm = TRUE) {
|
|
0
|
182 tmp <- tempdir()
|
|
|
183 message(sprintf("Downloading %s", url))
|
|
|
184 result <- file.path(tmp, target_file)
|
|
19
|
185 if (download.file(url, destfile = result) == 0) {
|
|
0
|
186 models <- source(result)$value
|
|
|
187 if (rm) {
|
|
19
|
188 unlink(tmp, recursive = TRUE)
|
|
0
|
189 }
|
|
19
|
190 return(models)
|
|
0
|
191 }
|
|
|
192 if (rm) {
|
|
19
|
193 unlink(tmp, recursive = TRUE)
|
|
0
|
194 }
|
|
|
195 stop("Could not download any file at this adress.")
|
|
|
196 }
|
|
|
197
|
|
|
198 search_tree <- function(path, target) {
|
|
|
199 target <- tolower(target)
|
|
|
200 for (file in list.files(path)) {
|
|
|
201 if (is.dir(file)) {
|
|
|
202 result <- search_tree(file.path(path, file), target)
|
|
|
203 if (!is.null(result)) {
|
|
19
|
204 return(result)
|
|
0
|
205 }
|
|
|
206 } else if (tolower(file) == target) {
|
|
19
|
207 return(file.path(path, file))
|
|
0
|
208 }
|
|
|
209 }
|
|
19
|
210 return(NULL)
|
|
0
|
211 }
|
|
|
212
|
|
|
213 create_database <- function(orm) {
|
|
19
|
214 orm$recreate_database(no_exists = FALSE)
|
|
0
|
215 set_database_version(orm, "created")
|
|
|
216 }
|
|
|
217
|
|
|
218 insert_adducts <- function(orm) {
|
|
|
219 message("Creating adducts...")
|
|
|
220 adducts <- list(
|
|
19
|
221 list("[M-H2O-H]-", 1, -1, -48.992020312000001069, 1, 0, 0.5, "H0", "H1O3"),
|
|
|
222 list("[M-H-Cl+O]-", 1, -1, -19.981214542000000022, 2, 0, 0.5, "O1", "H1Cl1"),
|
|
|
223 list("[M-Cl+O]-", 1, -1, -18.973389510000000512, 3, 0, 0.5, "O1", "Cl1"),
|
|
|
224 list("[M-3H]3-", 1, -3, -3.0218293560000000219, 4, 0, 1.0, "H0", "H3"),
|
|
|
225 list("[2M-3H]3-", 2, -3, -3.0218293560000000219, 4, 0, 0.5, "H0", "H3"),
|
|
|
226 list("[3M-3H]3-", 3, -3, -3.0218293560000000219, 4, 0, 0.5, "H0", "H3"),
|
|
|
227 list("[M-2H]2-", 1, -2, -2.0145529039999998666, 5, 0, 1.0, "H0", "H2"),
|
|
|
228 list("[2M-2H]2-", 2, -2, -2.0145529039999998666, 5, 0, 0.5, "H0", "H2"),
|
|
|
229 list("[3M-2H]2-", 3, -2, -2.0145529039999998666, 5, 0, 0.5, "H0", "H2"),
|
|
|
230 list("[M-H]-", 1, -1, -1.0072764519999999333, 6, 1, 1.0, "H0", "H1"),
|
|
|
231 list("[2M-H]-", 2, -1, -1.0072764519999999333, 6, 0, 0.5, "H0", "H1"),
|
|
|
232 list("[3M-H]-", 3, -1, -1.0072764519999999333, 6, 0, 0.5, "H0", "H1"),
|
|
|
233 list("[M]+", 1, 1, -0.00054858000000000000945, 7, 1, 1.0, "H0", "H0"),
|
|
|
234 list("[M]-", 1, -1, 0.00054858000000000000945, 8, 1, 1.0, "H0", "H0"),
|
|
|
235 list("[M+H]+", 1, 1, 1.0072764519999999333, 9, 1, 1.0, "H1", "H0"),
|
|
|
236 list("[2M+H]+", 2, 1, 1.0072764519999999333, 9, 0, 0.5, "H1", "H0"),
|
|
|
237 list("[3M+H]+", 3, 1, 1.0072764519999999333, 9, 0, 0.25, "H1", "H0"),
|
|
|
238 list("[M+2H]2+", 1, 2, 2.0145529039999998666, 10, 0, 0.75, "H2", "H0"),
|
|
|
239 list("[2M+2H]2+", 2, 2, 2.0145529039999998666, 10, 0, 0.5, "H2", "H0"),
|
|
|
240 list("[3M+2H]2+", 3, 2, 2.0145529039999998666, 10, 0, 0.25, "H2", "H0"),
|
|
|
241 list("[M+3H]3+", 1, 3, 3.0218293560000000219, 11, 0, 0.75, "H3", "H0"),
|
|
|
242 list("[2M+3H]3+", 2, 3, 3.0218293560000000219, 11, 0, 0.5, "H3", "H0"),
|
|
|
243 list("[3M+3H]3+", 3, 3, 3.0218293560000000219, 11, 0, 0.25, "H3", "H0"),
|
|
|
244 list("[M-2H+NH4]-", 1, -1, 16.019272654000001665, 12, 0, 0.25, "N1H4", "H2"),
|
|
|
245 list("[2M-2H+NH4]-", 2, -1, 16.019272654000001665, 12, 0, 0.0, "N1H4", "H2"),
|
|
|
246 list("[3M-2H+NH4]-", 3, -1, 16.019272654000001665, 12, 0, 0.25, "N1H4", "H2"),
|
|
|
247 list("[M+NH4]+", 1, 1, 18.033825558000000199, 13, 1, 1.0, "N1H4", "H0"),
|
|
|
248 list("[2M+NH4]+", 2, 1, 18.033825558000000199, 13, 0, 0.5, "N1H4", "H0"),
|
|
|
249 list("[3M+NH4]+", 3, 1, 18.033825558000000199, 13, 0, 0.25, "N1H4", "H0"),
|
|
|
250 list("[M+H+NH4]2+", 1, 2, 19.041102009999999467, 14, 0, 0.5, "N1H5", "H0"),
|
|
|
251 list("[2M+H+NH4]2+", 2, 2, 19.041102009999999467, 14, 0, 0.5, "N1H5", "H0"),
|
|
|
252 list("[3M+H+NH4]2+", 3, 2, 19.041102009999999467, 14, 0, 0.25, "N1H5", "H0"),
|
|
|
253 list("[M+Na-2H]-", 1, -1, 20.974668176000001551, 15, 0, 0.75, "Na1", "H2"),
|
|
|
254 list("[2M-2H+Na]-", 2, -1, 20.974668176000001551, 15, 0, 0.25, "Na1", "H2"),
|
|
|
255 list("[3M-2H+Na]-", 3, -1, 20.974668176000001551, 15, 0, 0.25, "Na1", "H2"),
|
|
|
256 list("[M+Na]+", 1, 1, 22.989221080000000086, 16, 1, 1.0, "Na1", "H0"),
|
|
|
257 list("[2M+Na]+", 2, 1, 22.989221080000000086, 16, 0, 0.5, "Na1", "H0"),
|
|
|
258 list("[3M+Na]+", 3, 1, 22.989221080000000086, 16, 0, 0.25, "Na1", "H0"),
|
|
|
259 list("[M+H+Na]2+", 1, 2, 23.996497531999999353, 17, 0, 0.5, "Na1H1", "H0"),
|
|
|
260 list("[2M+H+Na]2+", 2, 2, 23.996497531999999353, 17, 0, 0.5, "Na1H1", "H0"),
|
|
|
261 list("[3M+H+Na]2+", 3, 2, 23.996497531999999353, 17, 0, 0.25, "Na1H1", "H0"),
|
|
|
262 list("[M+2H+Na]3+", 1, 3, 25.003773983999998619, 18, 0, 0.25, "H2Na1", "H0"),
|
|
|
263 list("[M+CH3OH+H]+", 1, 1, 33.033491200000000276, 19, 0, 0.25, "C1O1H5", "H0"),
|
|
|
264 list("[M-H+Cl]2-", 1, -2, 33.962124838000001148, 20, 0, 1.0, "Cl1", "H1"),
|
|
|
265 list("[2M-H+Cl]2-", 2, -2, 33.962124838000001148, 20, 0, 0.5, "Cl1", "H1"),
|
|
|
266 list("[3M-H+Cl]2-", 3, -2, 33.962124838000001148, 20, 0, 0.5, "Cl1", "H1"),
|
|
|
267 list("[M+Cl]-", 1, -1, 34.969401290000000416, 21, 1, 1.0, "Cl1", "H0"),
|
|
|
268 list("[2M+Cl]-", 2, -1, 34.969401290000000416, 21, 0, 0.5, "Cl1", "H0"),
|
|
|
269 list("[3M+Cl]-", 3, -1, 34.969401290000000416, 21, 0, 0.5, "Cl1", "H0"),
|
|
|
270 list("[M+K-2H]-", 1, -1, 36.948605415999999479, 22, 0, 0.5, "K1", "H2"),
|
|
|
271 list("[2M-2H+K]-", 2, -1, 36.948605415999999479, 22, 0, 0.0, "K1", "H2"),
|
|
|
272 list("[3M-2H+K]-", 3, -1, 36.948605415999999479, 22, 0, 0.0, "K1", "H2"),
|
|
|
273 list("[M+K]+", 1, 1, 38.963158319999998013, 23, 1, 1.0, "K1", "H0"),
|
|
|
274 list("[2M+K]+", 2, 1, 38.963158319999998013, 23, 0, 0.5, "K1", "H0"),
|
|
|
275 list("[3M+K]+", 3, 1, 38.963158319999998013, 23, 0, 0.25, "K1", "H0"),
|
|
|
276 list("[M+H+K]2+", 1, 2, 39.970434771999997281, 24, 0, 0.5, "K1H1", "H0"),
|
|
|
277 list("[2M+H+K]2+", 2, 2, 39.970434771999997281, 24, 0, 0.5, "K1H1", "H0"),
|
|
|
278 list("[3M+H+K]2+", 3, 2, 39.970434771999997281, 24, 0, 0.25, "K1H1", "H0"),
|
|
|
279 list("[M+ACN+H]+", 1, 1, 42.033825557999996646, 25, 0, 0.25, "C2H4N1", "H0"),
|
|
|
280 list("[2M+ACN+H]+", 2, 1, 42.033825557999996646, 25, 0, 0.25, "C2H4N1", "H0"),
|
|
|
281 list("[M+2Na-H]+", 1, 1, 44.971165708000000902, 26, 0, 0.5, "Na2", "H1"),
|
|
|
282 list("[2M+2Na-H]+", 2, 1, 44.971165708000000902, 26, 0, 0.25, "Na2", "H1"),
|
|
|
283 list("[3M+2Na-H]+", 3, 1, 44.971165708000000902, 26, 0, 0.25, "Na2", "H1"),
|
|
|
284 list("[2M+FA-H]-", 2, -1, 44.998202851999998586, 27, 0, 0.25, "C1O2H2", "H1"),
|
|
|
285 list("[M+FA-H]-", 1, -1, 44.998202851999998586, 27, 0, 0.5, "C1O2H2", "H1"),
|
|
|
286 list("[M+2Na]2+", 1, 2, 45.978442160000000172, 28, 0, 0.5, "Na2", "H0"),
|
|
|
287 list("[2M+2Na]2+", 2, 2, 45.978442160000000172, 28, 0, 0.5, "Na2", "H0"),
|
|
|
288 list("[3M+2Na]2+", 3, 2, 45.978442160000000172, 28, 0, 0.25, "Na2", "H0"),
|
|
|
289 list("[M+H+2Na]3+", 1, 3, 46.985718611999999438, 29, 0, 0.25, "H1Na2", "H0"),
|
|
|
290 list("[M+H+FA]+", 1, 1, 47.012755755999997122, 30, 0, 0.25, "C1O2H3", "H0"),
|
|
|
291 list("[M+Hac-H]-", 1, -1, 59.013852915999997607, 31, 0, 0.25, "C2O2H4", "H1"),
|
|
|
292 list("[2M+Hac-H]-", 2, -1, 59.013852915999997607, 31, 0, 0.25, "C2O2H4", "H1"),
|
|
|
293 list("[M+IsoProp+H]+", 1, 1, 61.064791327999998317, 32, 0, 0.25, "C3H9O1", "H0"),
|
|
|
294 list("[M+Na+K]2+", 1, 2, 61.9523793999999981, 33, 0, 0.5, "Na1K1", "H0"),
|
|
|
295 list("[2M+Na+K]2+", 2, 2, 61.9523793999999981, 33, 0, 0.5, "Na1K1", "H0"),
|
|
|
296 list("[3M+Na+K]2+", 3, 2, 61.9523793999999981, 33, 0, 0.25, "Na1K1", "H0"),
|
|
|
297 list("[M+NO3]-", 1, -1, 61.988366450000000895, 34, 0, 0.5, "N1O3", "H0"),
|
|
|
298 list("[M+ACN+Na]+", 1, 1, 64.015770185999997464, 35, 0, 0.25, "C2H3N1Na1", "H0"),
|
|
|
299 list("[2M+ACN+Na]+", 2, 1, 64.015770185999997464, 35, 0, 0.25, "C2H3N1Na1", "H0"),
|
|
|
300 list("[M+NH4+FA]+", 1, 1, 64.039304861999994502, 36, 0, 0.25, "N1C1O2H6", "H0"),
|
|
|
301 list("[M-2H+Na+FA]-", 1, -1, 66.980147479999999405, 37, 0, 0.5, "NaC1O2H2", "H2"),
|
|
|
302 list("[M+3Na]3+", 1, 3, 68.967663239999993153, 38, 0, 0.25, "Na3", "H0"),
|
|
|
303 list("[M+Na+FA]+", 1, 1, 68.99470038399999794, 39, 0, 0.25, "Na1C1O2H2", "H0"),
|
|
|
304 list("[M+2Cl]2-", 1, -2, 69.938802580000000832, 40, 0, 1.0, "Cl2", "H0"),
|
|
|
305 list("[2M+2Cl]2-", 2, -2, 69.938802580000000832, 40, 0, 0.5, "Cl2", "H0"),
|
|
|
306 list("[3M+2Cl]2-", 3, -2, 69.938802580000000832, 40, 0, 0.5, "Cl2", "H0"),
|
|
|
307 list("[M+2K-H]+", 1, 1, 76.919040187999996758, 41, 0, 0.5, "K2", "H1"),
|
|
|
308 list("[2M+2K-H]+", 2, 1, 76.919040187999996758, 41, 0, 0.25, "K2", "H1"),
|
|
|
309 list("[3M+2K-H]+", 3, 1, 76.919040187999996758, 41, 0, 0.25, "K2", "H1"),
|
|
|
310 list("[M+2K]2+", 1, 2, 77.926316639999996028, 42, 0, 0.5, "K2", "H0"),
|
|
|
311 list("[2M+2K]2+", 2, 2, 77.926316639999996028, 42, 0, 0.5, "K2", "H0"),
|
|
|
312 list("[3M+2K]2+", 3, 2, 77.926316639999996028, 42, 0, 0.25, "K2", "H0"),
|
|
|
313 list("[M+Br]-", 1, -1, 78.918886479999997619, 43, 1, 1.0, "Br1", "H0"),
|
|
|
314 list("[M+Cl+FA]-", 1, -1, 80.974880593999998268, 44, 0, 0.5, "Cl1C1O2H2", "H0"),
|
|
|
315 list("[M+AcNa-H]-", 1, -1, 80.995797543999998426, 45, 0, 0.25, "C2H3Na1O2", "H1"),
|
|
|
316 list("[M+2ACN+2H]2+", 1, 2, 84.067651115999993292, 46, 0, 0.25, "C4H8N2", "H0"),
|
|
|
317 list("[M+K+FA]+", 1, 1, 84.968637623999995868, 47, 0, 0.25, "K1C1O2H2", "H0"),
|
|
|
318 list("[M+Cl+Na+FA-H]-", 1, -1, 102.95682522200000619, 48, 0, 0.5, "Cl1Na1C1O2H2", "H1"),
|
|
|
319 list("[2M+3H2O+2H]+", 2, 1, 104.03153939599999944, 49, 0, 0.25, "H8O6", "H0"),
|
|
|
320 list("[M+TFA-H]-", 1, -1, 112.98558742000000165, 50, 0, 0.5, "C2F3O2H1", "H1"),
|
|
|
321 list("[M+H+TFA]+", 1, 1, 115.00014032400000019, 51, 0, 0.25, "C2F3O2H2", "H0"),
|
|
|
322 list("[M+3ACN+2H]2+", 1, 2, 125.09420022199999778, 52, 0, 0.25, "C6H11N3", "H0"),
|
|
|
323 list("[M+NH4+TFA]+", 1, 1, 132.02668943000000468, 53, 0, 0.25, "N1C2F3O2H5", "H0"),
|
|
|
324 list("[M+Na+TFA]+", 1, 1, 136.98208495200000811, 54, 0, 0.25, "Na1C2F3O2H1", "H0"),
|
|
|
325 list("[M+Cl+TFA]-", 1, -1, 148.96226516199999423, 55, 0, 0.5, "Cl1C2F3O2H1", "H0"),
|
|
|
326 list("[M+K+TFA]+", 1, 1, 152.95602219200000604, 56, 0, 0.25, "K1C2F3O2H1","H0")
|
|
0
|
327 )
|
|
|
328 dummy_adduct <- orm$adduct()
|
|
|
329 for (adduct in adducts) {
|
|
|
330 i <- 0
|
|
19
|
331 dummy_adduct$set_name(adduct[[i <- i + 1]])
|
|
|
332 dummy_adduct$set_multi(adduct[[i <- i + 1]])
|
|
|
333 dummy_adduct$set_charge(adduct[[i <- i + 1]])
|
|
|
334 dummy_adduct$set_mass(adduct[[i <- i + 1]])
|
|
|
335 dummy_adduct$set_oidscore(adduct[[i <- i + 1]])
|
|
|
336 dummy_adduct$set_quasi(adduct[[i <- i + 1]])
|
|
|
337 dummy_adduct$set_ips(adduct[[i <- i + 1]])
|
|
|
338 dummy_adduct$set_formula_add(adduct[[i <- i + 1]])
|
|
|
339 dummy_adduct$set_formula_ded(adduct[[i <- i + 1]])
|
|
11
|
340 invisible(dummy_adduct$save())
|
|
19
|
341 dummy_adduct$clear(unset_id = TRUE)
|
|
0
|
342 }
|
|
|
343 message("Adducts created")
|
|
|
344 }
|
|
|
345
|
|
19
|
346 insert_base_data <- function(orm, path, archetype = FALSE) {
|
|
0
|
347 if (archetype) {
|
|
|
348 ## not implemented yet
|
|
19
|
349 return()
|
|
0
|
350 }
|
|
|
351 base_data <- readLines(path)
|
|
19
|
352 for (sql in strsplit(paste(base_data, collapse = " "), ";")[[1]]) {
|
|
0
|
353 orm$execute(sql)
|
|
|
354 }
|
|
|
355 set_database_version(orm, "enriched")
|
|
|
356 }
|
|
|
357
|
|
|
358 insert_compounds <- function(orm, compounds_path) {
|
|
19
|
359 compounds <- read.csv(file = compounds_path, sep = "\t")
|
|
0
|
360 if (is.null(compounds <- translate_compounds(compounds))) {
|
|
|
361 stop("Could not find asked compound's attributes in csv file.")
|
|
|
362 }
|
|
|
363 dummy_compound <- orm$compound()
|
|
|
364 compound_list <- list()
|
|
|
365 for (i in seq_len(nrow(compounds))) {
|
|
|
366 dummy_compound$set_mz(compounds[i, "mz"])
|
|
|
367 dummy_compound$set_name(compounds[i, "name"])
|
|
|
368 dummy_compound$set_common_name(compounds[i, "common_name"])
|
|
|
369 dummy_compound$set_formula(compounds[i, "formula"])
|
|
19
|
370 compound_list[[length(compound_list) + 1]] <- as.list(
|
|
0
|
371 dummy_compound,
|
|
|
372 c("mz", "name", "common_name", "formula")
|
|
|
373 )
|
|
19
|
374 dummy_compound$clear(unset_id = TRUE)
|
|
0
|
375 }
|
|
19
|
376 invisible(dummy_compound$save(bulk = compound_list))
|
|
0
|
377 }
|
|
|
378
|
|
|
379 translate_compounds <- function(compounds) {
|
|
|
380 recognized_headers <- list(
|
|
19
|
381 c(
|
|
|
382 "HMDB_ID", "MzBank", "X.M.H..", "X.M.H...1",
|
|
|
383 "MetName", "ChemFormula", "INChIkey"
|
|
|
384 )
|
|
0
|
385 )
|
|
|
386 header_translators <- list(
|
|
|
387 hmdb_header_translator
|
|
|
388 )
|
|
|
389 for (index in seq_along(recognized_headers)) {
|
|
|
390 headers <- recognized_headers[[index]]
|
|
|
391 if (identical(colnames(compounds), headers)) {
|
|
19
|
392 return(header_translators[[index]](compounds))
|
|
0
|
393 }
|
|
|
394 }
|
|
|
395 if (is.null(translator <- guess_translator(colnames(compounds)))) {
|
|
19
|
396 return(NULL)
|
|
0
|
397 }
|
|
19
|
398 return(csv_header_translator(translator, compounds))
|
|
0
|
399 }
|
|
|
400
|
|
|
401 guess_translator <- function(header) {
|
|
|
402 result <- list(
|
|
19
|
403 # HMDB_ID = NULL,
|
|
|
404 mz = NULL,
|
|
|
405 name = NULL,
|
|
|
406 common_name = NULL,
|
|
|
407 formula = NULL,
|
|
|
408 # inchi_key = NULL
|
|
0
|
409 )
|
|
|
410 asked_cols <- names(result)
|
|
|
411 for (asked_col in asked_cols) {
|
|
|
412 for (col in header) {
|
|
|
413 if ((twisted <- tolower(col)) == asked_col
|
|
|
414 || gsub("-", "_", twisted) == asked_col
|
|
|
415 || gsub(" ", "_", twisted) == asked_col
|
|
|
416 || tolower(gsub("(.)([A-Z])", "\\1_\\2", col)) == asked_col
|
|
|
417 ) {
|
|
|
418 result[[asked_col]] <- col
|
|
|
419 next
|
|
|
420 }
|
|
|
421 }
|
|
|
422 }
|
|
|
423 if (any(mapply(is.null, result))) {
|
|
19
|
424 return(NULL)
|
|
0
|
425 }
|
|
19
|
426 return(result)
|
|
0
|
427 }
|
|
|
428
|
|
|
429 hmdb_header_translator <- function(compounds) {
|
|
19
|
430 return(csv_header_translator(
|
|
0
|
431 list(
|
|
19
|
432 HMDB_ID = "HMDB_ID",
|
|
|
433 mz = "MzBank",
|
|
|
434 name = "MetName",
|
|
|
435 common_name = "MetName",
|
|
|
436 formula = "ChemFormula",
|
|
|
437 inchi_key = "INChIkey"
|
|
0
|
438 ), compounds
|
|
|
439 ))
|
|
|
440 }
|
|
|
441
|
|
|
442 csv_header_translator <- function(translation_table, csv) {
|
|
|
443 header_names <- names(translation_table)
|
|
19
|
444 result <- data.frame(seq_len(nrow(csv)))
|
|
0
|
445 for (i in seq_along(header_names)) {
|
|
|
446 result[, header_names[[i]]] <- csv[, translation_table[[i]]]
|
|
|
447 }
|
|
|
448 result[, "mz"] <- as.numeric(result[, "mz"])
|
|
19
|
449 return(result)
|
|
0
|
450 }
|
|
|
451
|
|
|
452 set_database_version <- function(orm, version) {
|
|
|
453 orm$set_tag(
|
|
|
454 version,
|
|
19
|
455 tag_name = "database_version",
|
|
|
456 tag_table_name = "XSeeker_tagging_table"
|
|
0
|
457 )
|
|
|
458 }
|
|
|
459
|
|
|
460 process_rdata <- function(orm, rdata, options) {
|
|
|
461 mzml_tmp_dir <- gather_mzml_files(rdata)
|
|
|
462 samples <- names(rdata$singlefile)
|
|
|
463 if (!is.null(options$samples)) {
|
|
|
464 samples <- samples[options$samples %in% samples]
|
|
|
465 }
|
|
|
466 show_percent <- (
|
|
|
467 is.null(options$`not-show-percent`)
|
|
|
468 || options$`not-show-percent` == FALSE
|
|
|
469 )
|
|
|
470 error <- tryCatch({
|
|
|
471 process_sample_list(
|
|
|
472 orm, rdata, samples,
|
|
19
|
473 show_percent = show_percent,
|
|
|
474 file_grouping_var = options$class
|
|
0
|
475 )
|
|
|
476 NULL
|
|
19
|
477 }, error = function(e) {
|
|
0
|
478 message(e)
|
|
|
479 e
|
|
|
480 })
|
|
|
481 if (!is.null(mzml_tmp_dir)) {
|
|
19
|
482 unlink(mzml_tmp_dir, recursive = TRUE)
|
|
0
|
483 }
|
|
|
484 if (!is.null(error)) {
|
|
|
485 stop(error)
|
|
|
486 }
|
|
|
487 }
|
|
|
488
|
|
|
489 gather_mzml_files <- function(rdata) {
|
|
|
490 if (is.null(rdata$singlefile)) {
|
|
|
491 message("Extracting mxml files")
|
|
|
492 tmp <- tempdir()
|
|
19
|
493 rdata$singlefile <- utils::unzip(rdata$zipfile, exdir = tmp)
|
|
|
494 names(rdata$singlefile) <- tools::file_path_sans_ext(
|
|
|
495 basename(rdata$singlefile)
|
|
|
496 )
|
|
0
|
497 message("Extracted")
|
|
19
|
498 return(tmp)
|
|
0
|
499 } else {
|
|
19
|
500 message(sprintf(
|
|
|
501 "Not a zip file, loading files directly from path: %s",
|
|
|
502 paste(rdata$singlefile, collapse = " ; ")
|
|
|
503 ))
|
|
0
|
504 }
|
|
19
|
505 return(NULL)
|
|
0
|
506 }
|
|
|
507
|
|
19
|
508 process_sample_list <- function(
|
|
|
509 orm,
|
|
|
510 rdata,
|
|
|
511 sample_names,
|
|
|
512 show_percent,
|
|
|
513 file_grouping_var = NULL
|
|
|
514 ) {
|
|
|
515 if (is.null(file_grouping_var)) {
|
|
|
516 file_grouping_var <- find_grouping_var(rdata$variableMetadata)
|
|
|
517 if (is.null(file_grouping_var)) {
|
|
|
518 stop("Malformed variableMetada.")
|
|
|
519 }
|
|
|
520 }
|
|
|
521 tryCatch({
|
|
|
522 headers <- colnames(rdata$variableMetadata)
|
|
|
523 file_grouping_var <- headers[[as.numeric(file_grouping_var)]]
|
|
|
524 }, error = function(e) NULL)
|
|
|
525 if (
|
|
|
526 is.null(file_grouping_var)
|
|
|
527 || !(file_grouping_var %in% colnames(rdata$variableMetadata))
|
|
|
528 ) {
|
|
|
529 stop(sprintf(
|
|
|
530 "Could not find grouping variable %s in var meta file.",
|
|
|
531 file_grouping_var
|
|
|
532 ))
|
|
|
533 }
|
|
0
|
534 message("Processing samples.")
|
|
|
535 message(sprintf("File grouping variable: %s", file_grouping_var))
|
|
|
536
|
|
11
|
537 context <- new.env()
|
|
|
538 context$samples <- list()
|
|
|
539 context$peaks <- rdata$xa@xcmsSet@peaks
|
|
|
540 context$groupidx <- rdata$xa@xcmsSet@groupidx
|
|
|
541 xcms_set <- rdata$xa@xcmsSet
|
|
|
542 singlefile <- rdata$singlefile
|
|
0
|
543 process_arg_list <- rdata$listOFlistArguments
|
|
11
|
544 var_meta <- rdata$variableMetadata
|
|
|
545
|
|
0
|
546 process_params <- list()
|
|
12
|
547 if (is.null(process_arg_list)) {
|
|
|
548 for (history in xcms_set@.processHistory) {
|
|
|
549 if (
|
|
|
550 class(history@param) == "CentWaveParam"
|
|
|
551 && history@type == "Peak detection"
|
|
|
552 ) {
|
|
|
553 params <- history@param
|
|
|
554 process_params <- list(list(
|
|
19
|
555 xfunction = "annotatediff",
|
|
|
556 ppm = params@ppm,
|
|
|
557 peakwidth = sprintf(
|
|
|
558 "%s - %s",
|
|
|
559 params@peakwidth[[1]],
|
|
|
560 params@peakwidth[[2]]
|
|
|
561 ),
|
|
|
562 snthresh = params@snthresh,
|
|
|
563 prefilterStep = params@prefilter[[1]],
|
|
|
564 prefilterLevel = params@prefilter[[2]],
|
|
|
565 mzdiff = params@mzdiff,
|
|
|
566 fitgauss = params@fitgauss,
|
|
|
567 noise = params@noise,
|
|
|
568 mzCenterFun = params@mzCenterFun,
|
|
|
569 integrate = params@integrate,
|
|
|
570 firstBaselineCheck = params@firstBaselineCheck,
|
|
|
571 snthreshIsoROIs = !identical(params@roiScales, numeric(0))
|
|
12
|
572 ))
|
|
|
573 break
|
|
|
574 }
|
|
0
|
575 }
|
|
12
|
576 } else {
|
|
|
577 for (list_name in names(process_arg_list)) {
|
|
|
578 param_list <- list()
|
|
|
579 for (param_name in names(process_arg_list[[list_name]])) {
|
|
19
|
580 param_list[[param_name]] <- process_arg_list[[
|
|
|
581 list_name
|
|
|
582 ]][[param_name]]
|
|
12
|
583 }
|
|
19
|
584 process_params[[length(process_params) + 1]] <- param_list
|
|
12
|
585 }
|
|
0
|
586 }
|
|
12
|
587
|
|
0
|
588 message("Parameters from previous processes extracted.")
|
|
|
589
|
|
|
590
|
|
|
591 indices <- as.numeric(unique(var_meta[, file_grouping_var]))
|
|
19
|
592 if (any(is.null(names(singlefile)[indices]))) {
|
|
|
593 stop(sprintf(
|
|
|
594 paste(
|
|
|
595 "Indices defined by grouping variable %s are not all present",
|
|
|
596 "in singlefile names (%s).\nCannot continue. Indices: %s"
|
|
|
597 ),
|
|
|
598 file_grouping_var,
|
|
|
599 paste(names(singlefile), collapse = ", "),
|
|
|
600 paste(indices, collapse = ", ")
|
|
|
601 ))
|
|
|
602 }
|
|
0
|
603 smol_xcms_set <- orm$smol_xcms_set()
|
|
|
604 mz_tab_info <- new.env()
|
|
|
605 g <- xcms::groups(xcms_set)
|
|
|
606 mz_tab_info$group_length <- nrow(g)
|
|
|
607 mz_tab_info$dataset_path <- xcms::filepaths(xcms_set)
|
|
|
608 mz_tab_info$sampnames <- xcms::sampnames(xcms_set)
|
|
|
609 mz_tab_info$sampclass <- xcms::sampclass(xcms_set)
|
|
19
|
610 mz_tab_info$rtmed <- g[, "rtmed"]
|
|
|
611 mz_tab_info$mzmed <- g[, "mzmed"]
|
|
|
612 mz_tab_info$smallmolecule_abundance_assay <- xcms::groupval(
|
|
|
613 xcms_set,
|
|
|
614 value = "into"
|
|
|
615 )
|
|
|
616 blogified <- blob::blob(fst::compress_fst(
|
|
|
617 serialize(mz_tab_info, NULL),
|
|
|
618 compression = 100
|
|
|
619 ))
|
|
11
|
620 rm(mz_tab_info)
|
|
|
621
|
|
|
622 invisible(smol_xcms_set$set_raw(blogified)$save())
|
|
|
623 smol_xcms_set_id <- smol_xcms_set$get_id()
|
|
|
624 rm(smol_xcms_set)
|
|
|
625
|
|
0
|
626 for (no in indices) {
|
|
11
|
627 sample_name <- names(singlefile)[[no]]
|
|
|
628 sample_path <- singlefile[[no]]
|
|
0
|
629 if (
|
|
|
630 is.na(no)
|
|
|
631 || is.null(sample_path)
|
|
|
632 || !(sample_name %in% sample_names)
|
|
|
633 ) {
|
|
|
634 next
|
|
|
635 }
|
|
|
636 env <- new.env()
|
|
11
|
637
|
|
14
|
638 if (!DEBUG_FAST_IGNORE_SLOW_OP) {
|
|
|
639 ms_file <- xcms::xcmsRaw(sample_path)
|
|
|
640 env$tic <- ms_file@tic
|
|
|
641 env$mz <- ms_file@env$mz
|
|
|
642 env$scanindex <- ms_file@scanindex
|
|
|
643 env$scantime <- ms_file@scantime * 60
|
|
|
644 env$intensity <- ms_file@env$intensity
|
|
|
645 env$polarity <- as.character(ms_file@polarity[[1]])
|
|
11
|
646
|
|
14
|
647 ## Again, ms file is huge, so we get rid of it quickly.
|
|
|
648 rm(ms_file)
|
|
|
649
|
|
|
650 }
|
|
0
|
651 env$sample_name <- sample_name
|
|
|
652 env$dataset_path <- sample_path
|
|
|
653 env$process_params <- process_params
|
|
|
654 env$enriched_rdata <- TRUE
|
|
|
655 env$enriched_rdata_version <- ENRICHED_RDATA_VERSION
|
|
|
656 env$tool_name <- TOOL_NAME
|
|
|
657 env$enriched_rdata_doc <- ENRICHED_RDATA_DOC
|
|
14
|
658
|
|
11
|
659 sample <- add_sample_to_database(orm, env, context, smol_xcms_set_id)
|
|
19
|
660 rm(env)
|
|
11
|
661 context$samples[no] <- sample$get_id()
|
|
19
|
662 rm(sample)
|
|
0
|
663 }
|
|
11
|
664 context$clusters <- list()
|
|
|
665 context$show_percent <- show_percent
|
|
|
666 context$cluster_mean_rt_abundance <- list()
|
|
|
667 context$central_feature <- list()
|
|
12
|
668 context$adducts <- list()
|
|
11
|
669 load_variable_metadata(orm, var_meta, context)
|
|
|
670 clusters <- context$clusters
|
|
|
671 rm(context)
|
|
0
|
672 message("Features enrichment")
|
|
11
|
673 complete_features(orm, clusters, show_percent)
|
|
0
|
674 message("Features enrichment done.")
|
|
19
|
675 return(NULL)
|
|
0
|
676 }
|
|
|
677
|
|
|
678 find_grouping_var <- function(var_meta) {
|
|
19
|
679 known_colnames <- c(
|
|
6
|
680 "name", "namecustom", "mz", "mzmin", "mzmax",
|
|
19
|
681 "rt", "rtmin", "rtmax", "npeaks", "isotopes", "adduct",
|
|
|
682 "pcgroup", "ms_level"
|
|
6
|
683 )
|
|
|
684 col_names <- colnames(var_meta)
|
|
19
|
685 classes <- list()
|
|
6
|
686 for (name in col_names) {
|
|
|
687 if (!(name %in% known_colnames)) {
|
|
19
|
688 classes[[length(classes) + 1]] <- name
|
|
0
|
689 }
|
|
|
690 }
|
|
6
|
691 if (length(classes) > 1) {
|
|
19
|
692 stop(sprintf(
|
|
|
693 "Only one class expected in the variable metadata. Found %d .",
|
|
|
694 length(classes)
|
|
|
695 ))
|
|
6
|
696 }
|
|
7
|
697 if (length(classes) == 0) {
|
|
6
|
698 stop("Could not find any class column in your variableMetadata.")
|
|
|
699 }
|
|
19
|
700 return(classes[[1]])
|
|
0
|
701 }
|
|
|
702
|
|
11
|
703 add_sample_to_database <- function(orm, env, context, smol_xcms_set_id) {
|
|
0
|
704 message(sprintf("Processing sample %s", env$sample_name))
|
|
|
705 sample <- (
|
|
|
706 orm$sample()
|
|
|
707 $set_name(env$sample_name)
|
|
|
708 $set_path(env$dataset_path)
|
|
|
709 $set_kind("enriched_rdata")
|
|
|
710 $set_polarity(
|
|
19
|
711 if (
|
|
|
712 is.null(env$polarity)
|
|
|
713 || identical(env$polarity, character(0))
|
|
|
714 ) ""
|
|
0
|
715 else env$polarity
|
|
|
716 )
|
|
|
717 $set_raw(blob::blob(fst::compress_fst(
|
|
|
718 serialize(env, NULL),
|
|
19
|
719 compression = 100
|
|
0
|
720 )))
|
|
|
721 )
|
|
11
|
722 sample[["smol_xcms_set_id"]] <- smol_xcms_set_id
|
|
|
723 sample$modified__[["smol_xcms_set_id"]] <- smol_xcms_set_id
|
|
|
724 sample <- sample$save()
|
|
0
|
725 load_process_params(orm, sample, env$process_params)
|
|
|
726 message(sprintf("Sample %s inserted.", env$sample_name))
|
|
19
|
727 return(sample)
|
|
0
|
728 }
|
|
|
729
|
|
|
730
|
|
11
|
731 load_variable_metadata <- function(orm, var_meta, context) {
|
|
0
|
732 all_clusters <- orm$cluster()$all()
|
|
|
733
|
|
11
|
734 next_feature_id <- get_next_id(orm$feature()$all(), "featureID") + 1
|
|
|
735 next_cluster_id <- 0
|
|
0
|
736 next_pc_group <- get_next_id(all_clusters, "pc_group")
|
|
11
|
737 next_align_group <- get_next_id(all_clusters, "align_group") + 1
|
|
0
|
738 message("Extracting features")
|
|
|
739 invisible(create_features(
|
|
11
|
740 orm, var_meta, context,
|
|
0
|
741 next_feature_id, next_cluster_id,
|
|
|
742 next_pc_group, next_align_group
|
|
|
743 ))
|
|
|
744 message("Extracting features done.")
|
|
19
|
745 return(NULL)
|
|
0
|
746 }
|
|
|
747
|
|
|
748 get_next_id <- function(models, attribute) {
|
|
|
749 if ((id <- models$max(attribute)) == Inf || id == -Inf) {
|
|
19
|
750 return(0)
|
|
0
|
751 }
|
|
19
|
752 return(id)
|
|
0
|
753 }
|
|
|
754
|
|
|
755 create_features <- function(
|
|
11
|
756 orm, var_meta, context,
|
|
0
|
757 next_feature_id, next_cluster_id,
|
|
|
758 next_pc_group, next_align_group
|
|
|
759 ) {
|
|
|
760 field_names <- as.list(names(orm$feature()$fields__))
|
|
19
|
761 field_names[field_names == "id"] <- NULL
|
|
0
|
762
|
|
|
763 features <- list()
|
|
|
764 dummy_feature <- orm$feature()
|
|
|
765
|
|
|
766 if (show_percent <- context$show_percent) {
|
|
|
767 percent <- -1
|
|
|
768 total <- nrow(var_meta)
|
|
|
769 }
|
|
14
|
770 rows <- seq_len(nrow(var_meta))
|
|
|
771 if (PROCESS_SMOL_BATCH) {
|
|
|
772
|
|
19
|
773 rows <- rows[1:as.integer(FAST_FEATURE_RATIO / 100.0 * length(rows))]
|
|
14
|
774 }
|
|
|
775 cluster_row <- list()
|
|
|
776 for (row in rows) {
|
|
0
|
777 if (show_percent && (row / total) * 100 > percent) {
|
|
|
778 percent <- percent + 1
|
|
19
|
779 message("\r", sprintf("\r%d %%", percent), appendLF = FALSE)
|
|
0
|
780 }
|
|
|
781
|
|
11
|
782 dummy_feature$set_featureID(next_feature_id)
|
|
|
783 next_feature_id <- next_feature_id + 1
|
|
14
|
784
|
|
|
785 curent_var_meta <- var_meta[row, ]
|
|
|
786 set_feature_fields_from_var_meta(dummy_feature, curent_var_meta)
|
|
11
|
787 fake_iso <- dummy_feature$get_iso()
|
|
|
788 iso <- extract_iso(fake_iso)
|
|
|
789 clusterID <- extract_clusterID(fake_iso, next_cluster_id)
|
|
|
790 context$clusterID <- clusterID
|
|
|
791 dummy_feature$set_iso(iso)
|
|
|
792
|
|
0
|
793 peak_list <- context$peaks[context$groupidx[[row]], ]
|
|
10
|
794 if (! ("matrix" %in% class(peak_list))) {
|
|
19
|
795 peak_list <- matrix(
|
|
|
796 peak_list,
|
|
|
797 nrow = 1,
|
|
|
798 ncol = length(peak_list),
|
|
|
799 dimnames = list(c(), names(peak_list))
|
|
|
800 )
|
|
10
|
801 }
|
|
11
|
802
|
|
|
803 clusterID <- as.character(clusterID)
|
|
|
804 if (is.null(context$central_feature[[clusterID]])) {
|
|
|
805 int_o <- extract_peak_var(peak_list, "into")
|
|
|
806 context$central_feature[[clusterID]] <- (
|
|
19
|
807 peak_list[peak_list[, "into"] == int_o, ]["sample"]
|
|
11
|
808 )
|
|
|
809 }
|
|
|
810
|
|
14
|
811 if (!DEBUG_FAST_IGNORE_SLOW_OP) {
|
|
19
|
812 central_feature <- context$central_feature[[clusterID]]
|
|
|
813 sample_peak_list <- peak_list[
|
|
|
814 as.integer(peak_list[, "sample"]) == central_feature,
|
|
|
815 ,
|
|
|
816 drop = FALSE
|
|
|
817 ]
|
|
|
818 if (
|
|
|
819 !identical(sample_peak_list, numeric(0))
|
|
|
820 && !is.null(nrow(sample_peak_list))
|
|
|
821 && nrow(sample_peak_list) != 0
|
|
|
822 ) {
|
|
|
823 int_o <- extract_peak_var(sample_peak_list, "into")
|
|
|
824 if (!is.na(int_o)) {
|
|
14
|
825 dummy_feature$set_int_o(int_o)
|
|
|
826 }
|
|
19
|
827 int_b <- extract_peak_var(sample_peak_list, "intb")
|
|
|
828 if (!is.na(int_b)) {
|
|
14
|
829 dummy_feature$set_int_b(int_b)
|
|
|
830 }
|
|
19
|
831 max_o <- extract_peak_var(sample_peak_list, "maxo")
|
|
|
832 if (!is.na(max_o)) {
|
|
14
|
833 dummy_feature$set_max_o(max_o)
|
|
|
834 }
|
|
0
|
835 }
|
|
|
836 }
|
|
14
|
837
|
|
|
838 cluster_row[[row]] <- create_associated_cluster(
|
|
12
|
839 orm,
|
|
14
|
840 context$samples[context$central_feature[[clusterID]]][[1]],
|
|
11
|
841 dummy_feature, clusterID,
|
|
0
|
842 context, curent_var_meta, next_pc_group,
|
|
|
843 next_align_group
|
|
|
844 )
|
|
|
845 next_align_group <- next_align_group + 1
|
|
19
|
846 features[[length(features) + 1]] <- as.list(dummy_feature, field_names)
|
|
0
|
847 dummy_feature$clear()
|
|
|
848 }
|
|
14
|
849 rm(var_meta)
|
|
|
850 message("")
|
|
0
|
851 message("Saving features")
|
|
19
|
852 invisible(dummy_feature$save(bulk = features))
|
|
14
|
853
|
|
|
854 ## We link manually clusters to the sample they're in.
|
|
|
855 link_cache <- list()
|
|
|
856 for (row in rows) {
|
|
|
857 sample_nos <- unique(context$peaks[context$groupidx[[row]], "sample"])
|
|
|
858 for (sample_id in context$samples[sample_nos]) {
|
|
|
859 cluster_id <- cluster_row[[row]]$get_id()
|
|
19
|
860 id <- paste(sample_id, cluster_id, sep = ";")
|
|
|
861 if (is.null(link_cache[[id]])) {
|
|
14
|
862 link_cache[[id]] <- 1
|
|
|
863 orm$cluster_sample(
|
|
19
|
864 sample_id = sample_id,
|
|
|
865 cluster_id = cluster_id
|
|
14
|
866 )$save()
|
|
|
867 }
|
|
|
868 }
|
|
|
869 }
|
|
|
870
|
|
0
|
871 message("Saved.")
|
|
19
|
872 return(context$clusters)
|
|
0
|
873 }
|
|
|
874
|
|
19
|
875 extract_peak_var <- function(peak_list, var_name, selector = max) {
|
|
0
|
876 value <- peak_list[, var_name]
|
|
|
877 names(value) <- NULL
|
|
19
|
878 return(selector(value))
|
|
0
|
879 }
|
|
|
880
|
|
|
881 set_feature_fields_from_var_meta <- function(feature, var_meta) {
|
|
|
882 if (!is.null(mz <- var_meta[["mz"]]) && !is.na(mz)) {
|
|
|
883 feature$set_mz(mz)
|
|
|
884 }
|
|
|
885 if (!is.null(mzmin <- var_meta[["mzmin"]]) && !is.na(mzmin)) {
|
|
|
886 feature$set_mz_min(mzmin)
|
|
|
887 }
|
|
|
888 if (!is.null(mzmax <- var_meta[["mzmax"]]) && !is.na(mzmax)) {
|
|
|
889 feature$set_mz_max(mzmax)
|
|
|
890 }
|
|
|
891 if (!is.null(rt <- var_meta[["rt"]]) && !is.na(rt)) {
|
|
|
892 feature$set_rt(rt)
|
|
|
893 }
|
|
|
894 if (!is.null(rtmin <- var_meta[["rtmin"]]) && !is.na(rtmin)) {
|
|
|
895 feature$set_rt_min(rtmin)
|
|
|
896 }
|
|
|
897 if (!is.null(rtmax <- var_meta[["rtmax"]]) && !is.na(rtmax)) {
|
|
|
898 feature$set_rt_max(rtmax)
|
|
|
899 }
|
|
|
900 if (!is.null(isotopes <- var_meta[["isotopes"]]) && !is.na(isotopes)) {
|
|
|
901 feature$set_iso(isotopes)
|
|
|
902 }
|
|
19
|
903 return(feature)
|
|
0
|
904 }
|
|
|
905
|
|
|
906 extract_iso <- function(weird_data) {
|
|
|
907 if (grepl("^\\[\\d+\\]", weird_data)[[1]]) {
|
|
19
|
908 return(sub("^\\[\\d+\\]", "", weird_data, perl = TRUE))
|
|
0
|
909 }
|
|
19
|
910 return(weird_data)
|
|
0
|
911 }
|
|
|
912
|
|
19
|
913 extract_clusterID <- function(weird_data, next_cluster_id) {
|
|
0
|
914 if (grepl("^\\[\\d+\\]", weird_data)[[1]]) {
|
|
|
915 clusterID <- stringr::str_extract(weird_data, "^\\[\\d+\\]")
|
|
|
916 clusterID <- as.numeric(stringr::str_extract(clusterID, "\\d+"))
|
|
|
917 } else {
|
|
|
918 clusterID <- 0
|
|
|
919 }
|
|
19
|
920 return(clusterID + next_cluster_id)
|
|
0
|
921 }
|
|
|
922
|
|
|
923 create_associated_cluster <- function(
|
|
12
|
924 orm,
|
|
14
|
925 main_sample_id, feature, clusterID,
|
|
0
|
926 context, curent_var_meta, next_pc_group, next_align_group
|
|
|
927 ) {
|
|
11
|
928 clusterID <- as.character(clusterID)
|
|
|
929 if (is.null(cluster <- context$clusters[[clusterID]])) {
|
|
12
|
930 pcgroup <- as.numeric(curent_var_meta[["pcgroup"]])
|
|
|
931 adduct_name <- as.character(curent_var_meta[["adduct"]])
|
|
|
932 annotation <- curent_var_meta[["isotopes"]]
|
|
11
|
933 cluster <- context$clusters[[clusterID]] <- orm$cluster(
|
|
19
|
934 pc_group = pcgroup + next_pc_group,
|
|
12
|
935 # adduct=adduct,
|
|
19
|
936 align_group = next_align_group,
|
|
0
|
937 # curent_group=curent_group,
|
|
19
|
938 clusterID = context$clusterID,
|
|
|
939 annotation = annotation
|
|
11
|
940 )
|
|
12
|
941 if (is.null(adduct <- context$adducts[[adduct_name]])) {
|
|
19
|
942 context$adducts[[adduct_name]] <- orm$adduct()$load_by(
|
|
|
943 name = adduct_name
|
|
|
944 )$first()
|
|
12
|
945 if (is.null(adduct <- context$adducts[[adduct_name]])) {
|
|
19
|
946 adduct <- context$adducts[[adduct_name]] <- orm$adduct(
|
|
|
947 name = adduct_name,
|
|
|
948 charge = 0
|
|
|
949 )
|
|
12
|
950 adduct$save()
|
|
|
951 }
|
|
|
952 }
|
|
|
953 cluster$set_adduct(adduct)
|
|
19
|
954 ## Crappy hack to assign sample id to cluster without loading the
|
|
|
955 ## sample. Samples are too big (their sample$env) and slows the
|
|
|
956 ## process, and eat all the menory so we dont't want to load them.
|
|
14
|
957 cluster[["sample_id"]] <- main_sample_id
|
|
|
958 cluster$modified__[["sample_id"]] <- main_sample_id
|
|
0
|
959 } else {
|
|
|
960 if (context$clusterID != 0 && cluster$get_clusterID() == 0) {
|
|
|
961 cluster$set_clusterID(context$clusterID)
|
|
|
962 }
|
|
|
963 }
|
|
|
964 cluster$save()
|
|
|
965 feature$set_cluster(cluster)
|
|
19
|
966 return(cluster)
|
|
0
|
967 }
|
|
|
968
|
|
11
|
969 complete_features <- function(orm, clusters, show_percent) {
|
|
|
970 total <- length(clusters)
|
|
|
971 percent <- -1
|
|
|
972 i <- 0
|
|
|
973 for (cluster in clusters) {
|
|
19
|
974 i <- i + 1
|
|
11
|
975 if (show_percent && (i / total) * 100 > percent) {
|
|
|
976 percent <- percent + 1
|
|
19
|
977 message("\r", sprintf("\r%d %%", percent), appendLF = FALSE)
|
|
11
|
978 }
|
|
19
|
979 features <- orm$feature()$load_by(cluster_id = cluster$get_id())
|
|
0
|
980 if (features$any()) {
|
|
|
981 if (!is.null(rt <- features$mean("rt"))) {
|
|
|
982 cluster$set_mean_rt(rt)$save()
|
|
|
983 }
|
|
|
984 features_df <- as.data.frame(features)
|
|
19
|
985 central_feature <- features_df[
|
|
|
986 grepl("^\\[M\\]", features_df[, "iso"]),
|
|
|
987 ]
|
|
0
|
988 central_feature_into <- central_feature[["int_o"]]
|
|
19
|
989 if (
|
|
|
990 !identical(central_feature_into, numeric(0))
|
|
|
991 && central_feature_into != 0
|
|
|
992 ) {
|
|
0
|
993 for (feature in as.vector(features)) {
|
|
|
994 feature$set_abundance(
|
|
|
995 feature$get_int_o() / central_feature_into * 100
|
|
|
996 )$save()
|
|
|
997 }
|
|
|
998 }
|
|
|
999 }
|
|
|
1000 }
|
|
19
|
1001 return(NULL)
|
|
0
|
1002 }
|
|
|
1003
|
|
|
1004 load_process_params <- function(orm, sample, params) {
|
|
|
1005 for (param_list in params) {
|
|
|
1006 if (is.null(param_list[["xfunction"]])) {
|
|
|
1007 next
|
|
|
1008 }
|
|
|
1009 if (param_list[["xfunction"]] == "annotatediff") {
|
|
|
1010 load_process_params_peak_picking(orm, sample, param_list)
|
|
|
1011 }
|
|
|
1012 }
|
|
19
|
1013 return(sample)
|
|
0
|
1014 }
|
|
|
1015
|
|
19
|
1016 load_process_params_peak_picking <- function(
|
|
|
1017 orm,
|
|
|
1018 sample,
|
|
|
1019 peak_picking_params
|
|
|
1020 ) {
|
|
|
1021 return(add_sample_process_parameters(
|
|
|
1022 params = peak_picking_params,
|
|
|
1023 params_translation = list(
|
|
|
1024 ppm = "ppm",
|
|
|
1025 maxcharge = "maxCharge",
|
|
|
1026 maxiso = "maxIso"
|
|
0
|
1027 ),
|
|
19
|
1028 param_model_generator = orm$peak_picking_parameters,
|
|
|
1029 sample_param_setter = sample$set_peak_picking_parameters
|
|
0
|
1030 ))
|
|
|
1031 }
|
|
|
1032
|
|
|
1033 add_sample_process_parameters <- function(
|
|
|
1034 params,
|
|
|
1035 params_translation,
|
|
|
1036 param_model_generator,
|
|
|
1037 sample_param_setter
|
|
|
1038 ) {
|
|
|
1039 model_params <- list()
|
|
|
1040 for (rdata_param_name in names(params_translation)) {
|
|
|
1041 database_param_name <- params_translation[[rdata_param_name]]
|
|
|
1042 if (is.null(rdata_param <- params[[rdata_param_name]])) {
|
|
|
1043 next
|
|
|
1044 }
|
|
|
1045 model_params[[database_param_name]] <- rdata_param
|
|
|
1046 }
|
|
|
1047 params_models <- do.call(param_model_generator()$load_by, model_params)
|
|
|
1048 if (params_models$any()) {
|
|
|
1049 params_model <- params_models$first()
|
|
|
1050 } else {
|
|
|
1051 params_model <- do.call(param_model_generator, model_params)
|
|
|
1052 params_model$save()
|
|
|
1053 }
|
|
19
|
1054 return(sample_param_setter(params_model)$save())
|
|
0
|
1055 }
|
|
|
1056
|
|
|
1057
|
|
|
1058 library(optparse)
|
|
|
1059
|
|
|
1060 option_list <- list(
|
|
|
1061 optparse::make_option(
|
|
|
1062 c("-v", "--version"),
|
|
19
|
1063 action = "store_true",
|
|
|
1064 help = "Display this tool's version and exits"
|
|
0
|
1065 ),
|
|
|
1066 optparse::make_option(
|
|
|
1067 c("-i", "--input"),
|
|
19
|
1068 type = "character",
|
|
|
1069 help = "The rdata path to import in XSeeker"
|
|
0
|
1070 ),
|
|
|
1071 optparse::make_option(
|
|
|
1072 c("-s", "--samples"),
|
|
19
|
1073 type = "character",
|
|
|
1074 help = "Samples to visualise in XSeeker"
|
|
0
|
1075 ),
|
|
|
1076 optparse::make_option(
|
|
|
1077 c("-B", "--archetype"),
|
|
19
|
1078 type = "character",
|
|
|
1079 help = "The name of the base database"
|
|
0
|
1080 ),
|
|
|
1081 optparse::make_option(
|
|
|
1082 c("-b", "--database"),
|
|
19
|
1083 type = "character",
|
|
|
1084 help = "The base database's path"
|
|
0
|
1085 ),
|
|
|
1086 optparse::make_option(
|
|
|
1087 c("-c", "--compounds-csv"),
|
|
19
|
1088 type = "character",
|
|
|
1089 help = "The csv containing compounds"
|
|
0
|
1090 ),
|
|
|
1091 optparse::make_option(
|
|
|
1092 c("-m", "--models"),
|
|
19
|
1093 type = "character",
|
|
|
1094 help = paste(
|
|
|
1095 "The path or url (must begin with http[s]:// or git@) to",
|
|
|
1096 "the database's models"
|
|
|
1097 )
|
|
0
|
1098 ),
|
|
|
1099 optparse::make_option(
|
|
19
|
1100 c("-k", "--class"),
|
|
|
1101 type = "character",
|
|
|
1102 help = "The name of the column containing the classes"
|
|
|
1103 ),
|
|
|
1104 optparse::make_option(
|
|
0
|
1105 c("-o", "--output"),
|
|
19
|
1106 type = "character",
|
|
|
1107 help = "The path where to output sqlite"
|
|
0
|
1108 ),
|
|
|
1109 optparse::make_option(
|
|
|
1110 c("-P", "--not-show-percent"),
|
|
19
|
1111 action = "store_true",
|
|
|
1112 help = "Flag not to show the percents",
|
|
|
1113 default = FALSE
|
|
0
|
1114 )
|
|
|
1115 )
|
|
|
1116
|
|
19
|
1117 options(error = function(){traceback(3)})
|
|
0
|
1118
|
|
19
|
1119 parser <- OptionParser(
|
|
|
1120 usage = "%prog [options] file",
|
|
|
1121 option_list = option_list
|
|
|
1122 )
|
|
|
1123 args <- parse_args(parser, positional_arguments = 0)
|
|
0
|
1124
|
|
|
1125 err_code <- 0
|
|
|
1126
|
|
|
1127 if (!is.null(args$options$version)) {
|
|
|
1128 message(sprintf("%s %s", TOOL_NAME, VERSION))
|
|
|
1129 quit()
|
|
|
1130 }
|
|
|
1131
|
|
|
1132 models <- get_models(args$options$models)
|
|
|
1133 orm <- DBModelR::ORM(
|
|
19
|
1134 connection_params = list(dbname=args$options$output),
|
|
|
1135 dbms = "SQLite"
|
|
0
|
1136 )
|
|
|
1137
|
|
|
1138 invisible(orm$models(models))
|
|
|
1139 invisible(create_database(orm))
|
|
|
1140
|
|
|
1141 message("Database model created")
|
|
|
1142
|
|
|
1143 insert_adducts(orm)
|
|
|
1144
|
|
|
1145 if (!is.null(args$options$database)) {
|
|
|
1146 insert_base_data(orm, args$options$database)
|
|
|
1147 }
|
|
|
1148 message(sprintf("Base data inserted using %s.", args$options$database))
|
|
|
1149
|
|
|
1150 if (!is.null(args$options$archetype)) {
|
|
19
|
1151 insert_base_data(orm, args$options$archetype, archetype = TRUE)
|
|
0
|
1152 }
|
|
|
1153 if (!is.null(args$options$`compounds-csv`)) {
|
|
|
1154 insert_compounds(orm, args$options$`compounds-csv`)
|
|
|
1155 }
|
|
|
1156
|
|
|
1157 # if (!is.null(args$options$rdata)) {
|
|
|
1158 # load_rdata_in_base(args$options$rdata, args$options$samples, args$options$`not-show-percent`)
|
|
|
1159 # }
|
|
|
1160
|
|
|
1161
|
|
|
1162 load(args$options$input, rdata <- new.env())
|
|
|
1163
|
|
|
1164 process_rdata(orm, rdata, args$options)
|
|
|
1165
|
|
19
|
1166 quit(status = err_code)
|
