Mercurial > repos > workflow4metabolomics > camera_combinexsannos

comparison lib.r @ 6:2443150161b4 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit cfad09eb4dd6b1439b7de6a0852cd8fa22210f58
author workflow4metabolomics
date Mon, 11 Sep 2023 22:41:22 +0000
parents 4c4945bce185
children
comparison
equal deleted inserted replaced
5:617cf0705775 6:2443150161b4
1 # lib.r 1 # lib.r
2 2
3 #@author G. Le Corguille 3 # @author G. Le Corguille
4 # solve an issue with batch if arguments are logical TRUE/FALSE 4 # solve an issue with batch if arguments are logical TRUE/FALSE
5 parseCommandArgs <- function(...) { 5 parseCommandArgs <- function(...) {
6 args <- batch::parseCommandArgs(...) 6 args <- batch::parseCommandArgs(...)
7 for (key in names(args)) { 7 for (key in names(args)) {
8 if (args[key] %in% c("TRUE", "FALSE")) 8 if (args[key] %in% c("TRUE", "FALSE")) {
9 args[key] <- as.logical(args[key]) 9 args[key] <- as.logical(args[key])
10 } 10 }
11 return(args) 11 }
12 } 12 return(args)
13 13 }
14 #@author G. Le Corguille 14
15 # @author G. Le Corguille
15 # This function will 16 # This function will
16 # - load the packages 17 # - load the packages
17 # - display the sessionInfo 18 # - display the sessionInfo
18 loadAndDisplayPackages <- function(pkgs) { 19 loadAndDisplayPackages <- function(pkgs) {
19 for (pkg in pkgs) suppressPackageStartupMessages(stopifnot(library(pkg, quietly = TRUE, logical.return = TRUE, character.only = TRUE))) 20 for (pkg in pkgs) suppressPackageStartupMessages(stopifnot(library(pkg, quietly = TRUE, logical.return = TRUE, character.only = TRUE)))
20 21
21 sessioninfo <- sessionInfo() 22 sessioninfo <- sessionInfo()
22 cat(sessioninfo$R.version$version.string, "\n") 23 cat(sessioninfo$R.version$version.string, "\n")
23 cat("Main packages:\n") 24 cat("Main packages:\n")
24 for (pkg in names(sessioninfo$otherPkgs)) { 25 for (pkg in names(sessioninfo$otherPkgs)) {
25 cat(paste(pkg, packageVersion(pkg)), "\t") 26 cat(paste(pkg, packageVersion(pkg)), "\t")
26 } 27 }
27 cat("\n") 28 cat("\n")
28 cat("Other loaded packages:\n") 29 cat("Other loaded packages:\n")
29 for (pkg in names(sessioninfo$loadedOnly)) { 30 for (pkg in names(sessioninfo$loadedOnly)) {
30 cat(paste(pkg, packageVersion(pkg)), "\t") 31 cat(paste(pkg, packageVersion(pkg)), "\t")
31 } 32 }
32 cat("\n") 33 cat("\n")
33 } 34 }
34 35
35 # This function retrieve a xset like object 36 # This function retrieve a xset like object
36 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr 37 # @author Gildas Le Corguille lecorguille@sb-roscoff.fr
37 getxcmsSetObject <- function(xobject) { 38 getxcmsSetObject <- function(xobject) {
38 # XCMS 1.x 39 # XCMS 1.x
39 if (class(xobject) == "xcmsSet") 40 if (class(xobject) == "xcmsSet") {
40 return(xobject) 41 return(xobject)
41 # XCMS 3.x 42 }
42 if (class(xobject) == "XCMSnExp") { 43 # XCMS 3.x
43 # Get the legacy xcmsSet object 44 if (class(xobject) == "XCMSnExp") {
44 suppressWarnings(xset <- as(xobject, "xcmsSet")) 45 # Get the legacy xcmsSet object
45 if (is.null(xset@phenoData$sample_group)) 46 suppressWarnings(xset <- as(xobject, "xcmsSet"))
46 sampclass(xset) <- "." 47 if (is.null(xset@phenoData$sample_group)) {
47 else 48 sampclass(xset) <- "."
48 sampclass(xset) <- xset@phenoData$sample_group 49 } else {
49 if (!is.null(xset@phenoData$sample_name)) 50 sampclass(xset) <- xset@phenoData$sample_group
50 rownames(xset@phenoData) <- xset@phenoData$sample_name 51 }
51 return(xset) 52 if (!is.null(xset@phenoData$sample_name)) {
52 } 53 rownames(xset@phenoData) <- xset@phenoData$sample_name
53 } 54 }
54 55 return(xset)
55 #@author G. Le Corguille 56 }
56 #The function create a pdf from the different png generated by diffreport 57 }
58
59 # @author G. Le Corguille
60 # The function create a pdf from the different png generated by diffreport
57 diffreport_png2pdf <- function(filebase) { 61 diffreport_png2pdf <- function(filebase) {
58 dir.create("pdf") 62 dir.create("pdf")
59 63
60 pdfEicOutput <- paste0("pdf/", filebase, "-eic_pdf.pdf") 64 pdfEicOutput <- paste0("pdf/", filebase, "-eic_pdf.pdf")
61 pdfBoxOutput <- paste0("pdf/", filebase, "-box_pdf.pdf") 65 pdfBoxOutput <- paste0("pdf/", filebase, "-box_pdf.pdf")
62 66
63 system(paste0("gm convert ", filebase, "_eic/*.png ", pdfEicOutput)) 67 system(paste0("gm convert ", filebase, "_eic/*.png ", pdfEicOutput))
64 system(paste0("gm convert ", filebase, "_box/*.png ", pdfBoxOutput)) 68 system(paste0("gm convert ", filebase, "_box/*.png ", pdfBoxOutput))
65 69 }
66 } 70
67 71 # @author G. Le Corguille
68 #@author G. Le Corguille 72 # The function create a zip archive from the different png generated by diffreport
69 #The function create a zip archive from the different png generated by diffreport
70 diffreport_png2zip <- function() { 73 diffreport_png2zip <- function() {
71 zip("eic.zip", dir(pattern = "_eic"), zip = Sys.which("zip")) 74 zip("eic.zip", dir(pattern = "_eic"), zip = Sys.which("zip"))
72 zip("box.zip", dir(pattern = "_box"), zip = Sys.which("zip")) 75 zip("box.zip", dir(pattern = "_box"), zip = Sys.which("zip"))
73 } 76 }
74 77
75 #The function create a zip archive from the different tabular generated by diffreport 78 # The function create a zip archive from the different tabular generated by diffreport
76 diffreport_tabular2zip <- function() { 79 diffreport_tabular2zip <- function() {
77 zip("tabular.zip", dir(pattern = "tabular/*"), zip = Sys.which("zip")) 80 zip("tabular.zip", dir(pattern = "tabular/*"), zip = Sys.which("zip"))
78 } 81 }
79 82
80 #@author G. Le Corguille 83 # @author G. Le Corguille
81 #This function convert if it is required the Retention Time in minutes 84 # This function convert if it is required the Retention Time in minutes
82 RTSecondToMinute <- function(variableMetadata, convertRTMinute) { 85 RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
83 if (convertRTMinute) { 86 if (convertRTMinute) {
84 #converting the retention times (seconds) into minutes 87 # converting the retention times (seconds) into minutes
85 print("converting the retention times into minutes in the variableMetadata") 88 print("converting the retention times into minutes in the variableMetadata")
86 variableMetadata[, "rt"] <- variableMetadata[, "rt"] / 60 89 variableMetadata[, "rt"] <- variableMetadata[, "rt"] / 60
87 variableMetadata[, "rtmin"] <- variableMetadata[, "rtmin"] / 60 90 variableMetadata[, "rtmin"] <- variableMetadata[, "rtmin"] / 60
88 variableMetadata[, "rtmax"] <- variableMetadata[, "rtmax"] / 60 91 variableMetadata[, "rtmax"] <- variableMetadata[, "rtmax"] / 60
89 } 92 }
90 return(variableMetadata) 93 return(variableMetadata)
91 } 94 }
92 95
93 #@author G. Le Corguille 96 # @author G. Le Corguille
94 #This function format ions identifiers 97 # This function format ions identifiers
95 formatIonIdentifiers <- function(variableMetadata, numDigitsRT = 0, numDigitsMZ = 0) { 98 formatIonIdentifiers <- function(variableMetadata, numDigitsRT = 0, numDigitsMZ = 0) {
96 splitDeco <- strsplit(as.character(variableMetadata$name), "_") 99 splitDeco <- strsplit(as.character(variableMetadata$name), "_")
97 idsDeco <- sapply(splitDeco, function(x) { 100 idsDeco <- sapply(splitDeco, function(x) {
98 deco <- unlist(x)[2]; if (is.na(deco)) return("") else return(paste0("_", deco)) 101 deco <- unlist(x)[2]
99 }) 102 if (is.na(deco)) {
100 namecustom <- make.unique(paste0("M", round(variableMetadata[, "mz"], numDigitsMZ), "T", round(variableMetadata[, "rt"], numDigitsRT), idsDeco)) 103 return("")
101 variableMetadata <- cbind(name = variableMetadata$name, namecustom = namecustom, variableMetadata[, !(colnames(variableMetadata) %in% c("name"))]) 104 } else {
102 return(variableMetadata) 105 return(paste0("_", deco))
103 } 106 }
104 107 })
105 #The function annotateDiffreport without the corr function which bugs 108 namecustom <- make.unique(paste0("M", round(variableMetadata[, "mz"], numDigitsMZ), "T", round(variableMetadata[, "rt"], numDigitsRT), idsDeco))
109 variableMetadata <- cbind(name = variableMetadata$name, namecustom = namecustom, variableMetadata[, !(colnames(variableMetadata) %in% c("name"))])
110 return(variableMetadata)
111 }
112
113 # The function annotateDiffreport without the corr function which bugs
106 annotatediff <- function(xset = xset, args = args, variableMetadataOutput = "variableMetadata.tsv") { 114 annotatediff <- function(xset = xset, args = args, variableMetadataOutput = "variableMetadata.tsv") {
107 # Resolve the bug with x11, with the function png 115 # Resolve the bug with x11, with the function png
108 options(bitmapType = "cairo") 116 options(bitmapType = "cairo")
109 117
110 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped. 118 # Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
119 res <- try(is.null(xset@filled))
120
121 # ------ annot -------
122 args$calcCiS <- as.logical(args$calcCiS)
123 args$calcIso <- as.logical(args$calcIso)
124 args$calcCaS <- as.logical(args$calcCaS)
125
126 # common parameters
127 args4annotate <- list(
128 object = xset,
129 nSlaves = args$nSlaves, sigma = args$sigma, perfwhm = args$perfwhm,
130 maxcharge = args$maxcharge, maxiso = args$maxiso, minfrac = args$minfrac,
131 ppm = args$ppm, mzabs = args$mzabs, quick = args$quick,
132 polarity = args$polarity, max_peaks = args$max_peaks, intval = args$intval
133 )
134
135 if (args$quick == FALSE) {
136 args4annotate <- append(
137 args4annotate,
138 list(
139 graphMethod = args$graphMethod, cor_eic_th = args$cor_eic_th, pval = args$pval,
140 calcCiS = args$calcCiS, calcIso = args$calcIso, calcCaS = args$calcCaS
141 )
142 )
143 # no ruleset
144 if (!is.null(args$multiplier)) {
145 args4annotate <- append(
146 args4annotate,
147 list(multiplier = args$multiplier)
148 )
149 } else { # ruleset
150 rulset <- read.table(args$rules, h = TRUE, sep = ";")
151 if (ncol(rulset) < 4) rulset <- read.table(args$rules, h = TRUE, sep = "\t")
152 if (ncol(rulset) < 4) rulset <- read.table(args$rules, h = TRUE, sep = ",")
153 if (ncol(rulset) < 4) {
154 error_message <- "Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation"
155 print(error_message)
156 stop(error_message)
157 }
158
159 args4annotate <- append(
160 args4annotate,
161 list(rules = rulset)
162 )
163 }
164 }
165
166
167 # launch annotate
168 xa <- do.call("annotate", args4annotate)
169 peakList <- getPeaklist(xa, intval = args$intval)
170 peakList <- cbind(groupnames(xa@xcmsSet), peakList)
171 colnames(peakList)[1] <- c("name")
172
173 # --- Multi condition : diffreport ---
174 diffrepOri <- NULL
175 if (!is.null(args$runDiffreport) && nlevels(sampclass(xset)) >= 2) {
176 # Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
111 res <- try(is.null(xset@filled)) 177 res <- try(is.null(xset@filled))
112 178 classes <- levels(sampclass(xset))
113 # ------ annot ------- 179 for (i in seq_len(length(classes) - 1)) {
114 args$calcCiS <- as.logical(args$calcCiS) 180 for (n in seq_len(length(classes))) {
115 args$calcIso <- as.logical(args$calcIso) 181 if (i + n <= length(classes)) {
116 args$calcCaS <- as.logical(args$calcCaS) 182 filebase <- paste(classes[i], class2 = classes[i + n], sep = "-vs-")
117 183
118 # common parameters 184 diffrep <- diffreport(
119 args4annotate <- list(object = xset, 185 object = xset, class1 = classes[i], class2 = classes[i + n],
120 nSlaves = args$nSlaves, sigma = args$sigma, perfwhm = args$perfwhm, 186 filebase = filebase, eicmax = args$eicmax, eicwidth = args$eicwidth,
121 maxcharge = args$maxcharge, maxiso = args$maxiso, minfrac = args$minfrac, 187 sortpval = TRUE, value = args$value, h = args$h, w = args$w, mzdec = args$mzdec, missing = 0
122 ppm = args$ppm, mzabs = args$mzabs, quick = args$quick, 188 )
123 polarity = args$polarity, max_peaks = args$max_peaks, intval = args$intval) 189
124 190 diffrepOri <- diffrep
125 if (args$quick == FALSE) { 191
126 args4annotate <- append(args4annotate, 192 # renamming of the column rtmed to rt to fit with camera peaklist function output
127 list(graphMethod = args$graphMethod, cor_eic_th = args$cor_eic_th, pval = args$pval, 193 colnames(diffrep)[colnames(diffrep) == "rtmed"] <- "rt"
128 calcCiS = args$calcCiS, calcIso = args$calcIso, calcCaS = args$calcCaS)) 194 colnames(diffrep)[colnames(diffrep) == "mzmed"] <- "mz"
129 # no ruleset 195
130 if (!is.null(args$multiplier)) { 196 # combines results and reorder columns
131 args4annotate <- append(args4annotate, 197 diffrep <- merge(peakList, diffrep[, c("name", "fold", "tstat", "pvalue")], by.x = "name", by.y = "name", sort = FALSE)
132 list(multiplier = args$multiplier)) 198 diffrep <- cbind(diffrep[, !(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))], diffrep[, (colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))])
199
200 diffrep <- RTSecondToMinute(diffrep, args$convertRTMinute)
201 diffrep <- formatIonIdentifiers(diffrep, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
202
203 if (args$sortpval) {
204 diffrep <- diffrep[order(diffrep$pvalue), ]
205 }
206
207 dir.create("tabular", showWarnings = FALSE)
208 write.table(diffrep, sep = "\t", quote = FALSE, row.names = FALSE, file = paste("tabular/", filebase, "_tsv.tabular", sep = ""))
209
210 if (args$eicmax != 0) {
211 if (args$png2 == "pdf") {
212 diffreport_png2pdf(filebase)
213 }
214 if (args$png2 == "zip") {
215 diffreport_png2zip()
216 }
217 }
133 } 218 }
134 # ruleset 219 }
135 else { 220 }
136 rulset <- read.table(args$rules, h = T, sep = ";") 221 if (args$tabular2 == "zip") {
137 if (ncol(rulset) < 4) rulset <- read.table(args$rules, h = T, sep = "\t") 222 diffreport_tabular2zip()
138 if (ncol(rulset) < 4) rulset <- read.table(args$rules, h = T, sep = ",") 223 }
139 if (ncol(rulset) < 4) { 224 }
140 error_message <- "Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation" 225
141 print(error_message) 226 # --- variableMetadata ---
142 stop(error_message) 227 variableMetadata <- peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
143 } 228 variableMetadata <- RTSecondToMinute(variableMetadata, args$convertRTMinute)
144 229 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
145 args4annotate <- append(args4annotate, 230 # if we have 2 conditions, we keep stat of diffrep
146 list(rules = rulset)) 231 if (!is.null(args$runDiffreport) && nlevels(sampclass(xset)) == 2) {
147 } 232 variableMetadata <- merge(variableMetadata, diffrep[, c("name", "fold", "tstat", "pvalue")], by.x = "name", by.y = "name", sort = FALSE)
148 } 233 if (exists("args[[\"sortpval\"]]")) {
149 234 variableMetadata <- variableMetadata[order(variableMetadata$pvalue), ]
150 235 }
151 # launch annotate 236 }
152 xa <- do.call("annotate", args4annotate) 237
153 peakList <- getPeaklist(xa, intval = args$intval) 238 variableMetadataOri <- variableMetadata
154 peakList <- cbind(groupnames(xa@xcmsSet), peakList); colnames(peakList)[1] <- c("name"); 239 write.table(variableMetadata, sep = "\t", quote = FALSE, row.names = FALSE, file = variableMetadataOutput)
155 240
156 # --- Multi condition : diffreport --- 241 return(list("xa" = xa, "diffrep" = diffrepOri, "variableMetadata" = variableMetadataOri))
157 diffrepOri <- NULL 242 }
158 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset)) >= 2) { 243
159 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped. 244
160 res <- try(is.null(xset@filled)) 245 combinexsAnnos_function <- function(
161 classes <- levels(sampclass(xset)) 246 xaP, xaN, diffrepP = NULL, diffrepN = NULL,
162 x <- 1:(length(classes) - 1) 247 pos = TRUE, tol = 2, ruleset = NULL, keep_meta = TRUE, convertRTMinute = FALSE, numDigitsMZ = 0,
163 for (i in seq(along = x)) {
164 y <- 1:(length(classes))
165 for (n in seq(along = y)) {
166 if (i + n <= length(classes)) {
167 filebase <- paste(classes[i], class2 = classes[i + n], sep = "-vs-")
168
169 diffrep <- diffreport(
170 object = xset, class1 = classes[i], class2 = classes[i + n],
171 filebase = filebase, eicmax = args$eicmax, eicwidth = args$eicwidth,
172 sortpval = TRUE, value = args$value, h = args$h, w = args$w, mzdec = args$mzdec, missing = 0)
173
174 diffrepOri <- diffrep
175
176 # renamming of the column rtmed to rt to fit with camera peaklist function output
177 colnames(diffrep)[colnames(diffrep) == "rtmed"] <- "rt"
178 colnames(diffrep)[colnames(diffrep) == "mzmed"] <- "mz"
179
180 # combines results and reorder columns
181 diffrep <- merge(peakList, diffrep[, c("name", "fold", "tstat", "pvalue")], by.x = "name", by.y = "name", sort = F)
182 diffrep <- cbind(diffrep[, !(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))], diffrep[, (colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))])
183
184 diffrep <- RTSecondToMinute(diffrep, args$convertRTMinute)
185 diffrep <- formatIonIdentifiers(diffrep, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
186
187 if (args$sortpval) {
188 diffrep <- diffrep[order(diffrep$pvalue), ]
189 }
190
191 dir.create("tabular", showWarnings = FALSE)
192 write.table(diffrep, sep = "\t", quote = FALSE, row.names = FALSE, file = paste("tabular/", filebase, "_tsv.tabular", sep = ""))
193
194 if (args$eicmax != 0) {
195 if (args$png2 == "pdf")
196 diffreport_png2pdf(filebase)
197 if (args$png2 == "zip")
198 diffreport_png2zip()
199 }
200 }
201 }
202 }
203 if (args$tabular2 == "zip")
204 diffreport_tabular2zip()
205 }
206
207 # --- variableMetadata ---
208 variableMetadata <- peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
209 variableMetadata <- RTSecondToMinute(variableMetadata, args$convertRTMinute)
210 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
211 # if we have 2 conditions, we keep stat of diffrep
212 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset)) == 2) {
213 variableMetadata <- merge(variableMetadata, diffrep[, c("name", "fold", "tstat", "pvalue")], by.x = "name", by.y = "name", sort = F)
214 if (exists("args[[\"sortpval\"]]")) {
215 variableMetadata <- variableMetadata[order(variableMetadata$pvalue), ]
216 }
217 }
218
219 variableMetadataOri <- variableMetadata
220 write.table(variableMetadata, sep = "\t", quote = FALSE, row.names = FALSE, file = variableMetadataOutput)
221
222 return(list("xa" = xa, "diffrep" = diffrepOri, "variableMetadata" = variableMetadataOri));
223
224 }
225
226
227 combinexsAnnos_function <- function(xaP, xaN, diffrepP = NULL, diffrepN = NULL,
228 pos = TRUE, tol = 2, ruleset = NULL, keep_meta = TRUE, convertRTMinute = F, numDigitsMZ = 0,
229 numDigitsRT = 0, variableMetadataOutput = "variableMetadata.tsv") { 248 numDigitsRT = 0, variableMetadataOutput = "variableMetadata.tsv") {
230 249 # Load the two Rdata to extract the xset objects from positive and negative mode
231 #Load the two Rdata to extract the xset objects from positive and negative mode 250 cat("\tObject xset from positive mode\n")
232 cat("\tObject xset from positive mode\n") 251 print(xaP)
233 print(xaP) 252 cat("\n")
234 cat("\n") 253
235 254 cat("\tObject xset from negative mode\n")
236 cat("\tObject xset from negative mode\n") 255 print(xaN)
237 print(xaN) 256 cat("\n")
238 cat("\n") 257
239 258 cat("\n")
240 cat("\n") 259 cat("\tCombining...\n")
241 cat("\tCombining...\n") 260 # Convert the string to numeric for creating matrix
242 #Convert the string to numeric for creating matrix 261 row <- as.numeric(strsplit(ruleset, ",")[[1]][1])
243 row <- as.numeric(strsplit(ruleset, ",")[[1]][1]) 262 column <- as.numeric(strsplit(ruleset, ",")[[1]][2])
244 column <- as.numeric(strsplit(ruleset, ",")[[1]][2]) 263 ruleset <- cbind(row, column)
245 ruleset <- cbind(row, column) 264 # Test if the file comes from an older version tool
246 #Test if the file comes from an older version tool 265 if ((!is.null(xaP)) && (!is.null(xaN))) {
247 if ((!is.null(xaP)) & (!is.null(xaN))) { 266 # Launch the combinexsannos function from CAMERA
248 #Launch the combinexsannos function from CAMERA 267 cAnnot <- combinexsAnnos(xaP, xaN, pos = pos, tol = tol, ruleset = ruleset)
249 cAnnot <- combinexsAnnos(xaP, xaN, pos = pos, tol = tol, ruleset = ruleset) 268 } else {
250 } else { 269 stop("You must relauch the CAMERA.annotate step with the lastest version.")
251 stop("You must relauch the CAMERA.annotate step with the lastest version.") 270 }
252 } 271
253 272 if (pos) {
254 if (pos) { 273 xa <- xaP
255 xa <- xaP 274 mode <- "neg. Mode"
256 mode <- "neg. Mode" 275 } else {
257 } else { 276 xa <- xaN
258 xa <- xaN 277 mode <- "pos. Mode"
259 mode <- "pos. Mode" 278 }
260 } 279
261 280 peakList <- getPeaklist(xa)
262 peakList <- getPeaklist(xa) 281 peakList <- cbind(groupnames(xa@xcmsSet), peakList)
263 peakList <- cbind(groupnames(xa@xcmsSet), peakList); colnames(peakList)[1] <- c("name"); 282 colnames(peakList)[1] <- c("name")
264 variableMetadata <- cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup", mode)]); 283 variableMetadata <- cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup", mode)])
265 variableMetadata <- variableMetadata[, !(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))] 284 variableMetadata <- variableMetadata[, !(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))]
266 285
267 #Test if there are more than two classes (conditions) 286 # Test if there are more than two classes (conditions)
268 if (nlevels(sampclass(xaP@xcmsSet)) == 2 & (!is.null(diffrepN)) & (!is.null(diffrepP))) { 287 if (nlevels(sampclass(xaP@xcmsSet)) == 2 && (!is.null(diffrepN)) && (!is.null(diffrepP))) {
269 diffrepP <- diffrepP[, c("name", "fold", "tstat", "pvalue")]; colnames(diffrepP) <- paste("P.", colnames(diffrepP), sep = "") 288 diffrepP <- diffrepP[, c("name", "fold", "tstat", "pvalue")]
270 diffrepN <- diffrepN[, c("name", "fold", "tstat", "pvalue")]; colnames(diffrepN) <- paste("N.", colnames(diffrepN), sep = "") 289 colnames(diffrepP) <- paste("P.", colnames(diffrepP), sep = "")
271 290 diffrepN <- diffrepN[, c("name", "fold", "tstat", "pvalue")]
272 variableMetadata <- merge(variableMetadata, diffrepP, by.x = "name", by.y = "P.name") 291 colnames(diffrepN) <- paste("N.", colnames(diffrepN), sep = "")
273 variableMetadata <- merge(variableMetadata, diffrepN, by.x = "name", by.y = "N.name") 292
274 } 293 variableMetadata <- merge(variableMetadata, diffrepP, by.x = "name", by.y = "P.name")
275 294 variableMetadata <- merge(variableMetadata, diffrepN, by.x = "name", by.y = "N.name")
276 rownames(variableMetadata) <- NULL 295 }
277 #TODO: checker colnames(variableMetadata)[1:2] = c("name", "mz/rt"); 296
278 297 rownames(variableMetadata) <- NULL
279 variableMetadata <- RTSecondToMinute(variableMetadata, convertRTMinute) 298 # TODO: checker colnames(variableMetadata)[1:2] = c("name", "mz/rt");
280 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = numDigitsRT, numDigitsMZ = numDigitsMZ) 299
281 300 variableMetadata <- RTSecondToMinute(variableMetadata, convertRTMinute)
282 #If the user want to keep only the metabolites which match a difference 301 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = numDigitsRT, numDigitsMZ = numDigitsMZ)
283 if (keep_meta) { 302
284 variableMetadata <- variableMetadata[variableMetadata[, c(mode)] != "", ] 303 # If the user want to keep only the metabolites which match a difference
285 } 304 if (keep_meta) {
286 305 variableMetadata <- variableMetadata[variableMetadata[, c(mode)] != "", ]
287 #Write the output into a tsv file 306 }
288 write.table(variableMetadata, sep = "\t", quote = FALSE, row.names = FALSE, file = variableMetadataOutput) 307
289 return(variableMetadata); 308 # Write the output into a tsv file
290 309 write.table(variableMetadata, sep = "\t", quote = FALSE, row.names = FALSE, file = variableMetadataOutput)
310 return(variableMetadata)
291 } 311 }
292 312
293 # This function get the raw file path from the arguments 313 # This function get the raw file path from the arguments
294 getRawfilePathFromArguments <- function(singlefile, zipfile, args) { 314 getRawfilePathFromArguments <- function(singlefile, zipfile, args) {
295 if (!is.null(args$zipfile)) zipfile <- args$zipfile 315 if (!is.null(args$zipfile)) zipfile <- args$zipfile
296 if (!is.null(args$zipfilePositive)) zipfile <- args$zipfilePositive 316 if (!is.null(args$zipfilePositive)) zipfile <- args$zipfilePositive
297 if (!is.null(args$zipfileNegative)) zipfile <- args$zipfileNegative 317 if (!is.null(args$zipfileNegative)) zipfile <- args$zipfileNegative
298 318
299 if (!is.null(args$singlefile_galaxyPath)) { 319 if (!is.null(args$singlefile_galaxyPath)) {
300 singlefile_galaxyPaths <- args$singlefile_galaxyPath; 320 singlefile_galaxyPaths <- args$singlefile_galaxyPath
301 singlefile_sampleNames <- args$singlefile_sampleName 321 singlefile_sampleNames <- args$singlefile_sampleName
302 } 322 }
303 if (!is.null(args$singlefile_galaxyPathPositive)) { 323 if (!is.null(args$singlefile_galaxyPathPositive)) {
304 singlefile_galaxyPaths <- args$singlefile_galaxyPathPositive; 324 singlefile_galaxyPaths <- args$singlefile_galaxyPathPositive
305 singlefile_sampleNames <- args$singlefile_sampleNamePositive 325 singlefile_sampleNames <- args$singlefile_sampleNamePositive
306 } 326 }
307 if (!is.null(args$singlefile_galaxyPathNegative)) { 327 if (!is.null(args$singlefile_galaxyPathNegative)) {
308 singlefile_galaxyPaths <- args$singlefile_galaxyPathNegative; 328 singlefile_galaxyPaths <- args$singlefile_galaxyPathNegative
309 singlefile_sampleNames <- args$singlefile_sampleNameNegative 329 singlefile_sampleNames <- args$singlefile_sampleNameNegative
310 } 330 }
311 if (exists("singlefile_galaxyPaths")) { 331 if (exists("singlefile_galaxyPaths")) {
312 singlefile_galaxyPaths <- unlist(strsplit(singlefile_galaxyPaths, ",")) 332 singlefile_galaxyPaths <- unlist(strsplit(singlefile_galaxyPaths, ","))
313 singlefile_sampleNames <- unlist(strsplit(singlefile_sampleNames, ",")) 333 singlefile_sampleNames <- unlist(strsplit(singlefile_sampleNames, ","))
314 334
315 singlefile <- NULL 335 singlefile <- NULL
316 for (singlefile_galaxyPath_i in seq_len(length(singlefile_galaxyPaths))) { 336 for (singlefile_galaxyPath_i in seq_len(length(singlefile_galaxyPaths))) {
317 singlefile_galaxyPath <- singlefile_galaxyPaths[singlefile_galaxyPath_i] 337 singlefile_galaxyPath <- singlefile_galaxyPaths[singlefile_galaxyPath_i]
318 singlefile_sampleName <- singlefile_sampleNames[singlefile_galaxyPath_i] 338 singlefile_sampleName <- singlefile_sampleNames[singlefile_galaxyPath_i]
319 singlefile[[singlefile_sampleName]] <- singlefile_galaxyPath 339 singlefile[[singlefile_sampleName]] <- singlefile_galaxyPath
320 } 340 }
321 } 341 }
322 for (argument in c("zipfile", "zipfilePositive", "zipfileNegative", 342 for (argument in c(
323 "singlefile_galaxyPath", "singlefile_sampleName", 343 "zipfile", "zipfilePositive", "zipfileNegative",
324 "singlefile_galaxyPathPositive", "singlefile_sampleNamePositive", 344 "singlefile_galaxyPath", "singlefile_sampleName",
325 "singlefile_galaxyPathNegative", "singlefile_sampleNameNegative")) { 345 "singlefile_galaxyPathPositive", "singlefile_sampleNamePositive",
326 args[[argument]] <- NULL 346 "singlefile_galaxyPathNegative", "singlefile_sampleNameNegative"
327 } 347 )) {
328 return(list(zipfile = zipfile, singlefile = singlefile, args = args)) 348 args[[argument]] <- NULL
349 }
350 return(list(zipfile = zipfile, singlefile = singlefile, args = args))
329 } 351 }
330 352
331 353
332 # This function retrieve the raw file in the working directory 354 # This function retrieve the raw file in the working directory
333 # - if zipfile: unzip the file with its directory tree 355 # - if zipfile: unzip the file with its directory tree
334 # - if singlefiles: set symlink with the good filename 356 # - if singlefiles: set symlink with the good filename
335 retrieveRawfileInTheWorkingDir <- function(singlefile, zipfile) { 357 retrieveRawfileInTheWorkingDir <- function(singlefile, zipfile) {
336 if (!is.null(singlefile) && (length("singlefile") > 0)) { 358 if (!is.null(singlefile) && (length("singlefile") > 0)) {
337 for (singlefile_sampleName in names(singlefile)) { 359 for (singlefile_sampleName in names(singlefile)) {
338 singlefile_galaxyPath <- singlefile[[singlefile_sampleName]] 360 singlefile_galaxyPath <- singlefile[[singlefile_sampleName]]
339 if (!file.exists(singlefile_galaxyPath)) { 361 if (!file.exists(singlefile_galaxyPath)) {
340 error_message <- paste("Cannot access the sample:", singlefile_sampleName, "located:", singlefile_galaxyPath, ". Please, contact your administrator ... if you have one!") 362 error_message <- paste("Cannot access the sample:", singlefile_sampleName, "located:", singlefile_galaxyPath, ". Please, contact your administrator ... if you have one!")
341 print(error_message); stop(error_message) 363 print(error_message)
342 } 364 stop(error_message)
343 365 }
344 file.symlink(singlefile_galaxyPath, singlefile_sampleName) 366
345 } 367 file.symlink(singlefile_galaxyPath, singlefile_sampleName)
346 directory <- "." 368 }
347 369 directory <- "."
348 } 370 }
349 if (!is.null(zipfile) && (zipfile != "")) { 371 if (!is.null(zipfile) && (zipfile != "")) {
350 if (!file.exists(zipfile)) { 372 if (!file.exists(zipfile)) {
351 error_message <- paste("Cannot access the Zip file:", zipfile, ". Please, contact your administrator ... if you have one!") 373 error_message <- paste("Cannot access the Zip file:", zipfile, ". Please, contact your administrator ... if you have one!")
352 print(error_message) 374 print(error_message)
353 stop(error_message) 375 stop(error_message)
354 } 376 }
355 377
356 #unzip 378 # unzip
357 suppressWarnings(unzip(zipfile, unzip = "unzip")) 379 suppressWarnings(unzip(zipfile, unzip = "unzip"))
358 380
359 #get the directory name 381 # get the directory name
360 filesInZip <- unzip(zipfile, list = T); 382 filesInZip <- unzip(zipfile, list = TRUE)
361 directories <- unique(unlist(lapply(strsplit(filesInZip$Name, "/"), function(x) x[1]))); 383 directories <- unique(unlist(lapply(strsplit(filesInZip$Name, "/"), function(x) x[1])))
362 directories <- directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir] 384 directories <- directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir]
363 directory <- "." 385 directory <- "."
364 if (length(directories) == 1) directory <- directories 386 if (length(directories) == 1) directory <- directories
365 387
366 cat("files_root_directory\t", directory, "\n") 388 cat("files_root_directory\t", directory, "\n")
367 389 }
368 } 390 return(directory)
369 return(directory) 391 }
370 }