Mercurial > repos > workflow4metabolomics > camera_groupfwhm

comparison CAMERA_groupFWHM.R @ 0:6aea0427511e draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 24d44ee26b7c23380c2b42fae2f7f6e58472100d
author workflow4metabolomics
date Sun, 24 Nov 2024 21:28:57 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6aea0427511e
1 #!/usr/bin/env Rscript
2
3 # ----- PACKAGE -----
4 cat("\tSESSION INFO\n")
5
6 # Import the different functions
7 source_local <- function(fname) {
8 argv <- commandArgs(trailingOnly = FALSE)
9 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
10 source(paste(base_dir, fname, sep = "/"))
11 }
12 source_local("lib.r")
13
14 pkgs <- c("CAMERA", "xcms", "multtest", "batch")
15 loadAndDisplayPackages(pkgs)
16 cat("\n\n")
17 # ----- ARGUMENTS -----
18 cat("\tARGUMENTS INFO\n")
19
20 args <- parseCommandArgs(evaluate = FALSE) # interpretation of arguments given in command line as an R list of objects
21 write.table(as.matrix(args), col.names = FALSE, quote = FALSE, sep = "\t")
22
23 cat("\n\n")
24
25 print("Arguments retrieved from the command line:")
26 print(args)
27
28 # Function to convert "NA" strings to actual NA values and string lists to numeric lists
29 convertStringToNumeric <- function(x) {
30 # Force conversion to character
31 x <- as.character(x)
32
33 if (x == "NA") {
34 return(NA)
35 } else if (grepl("^[0-9]+$", x)) {
36 # If the string represents a single numeric value
37 return(as.numeric(x))
38 } else {
39 # Convert string representation of a list to a numeric vector
40 # Use a regular expression to split by common separators
41 return(as.numeric(unlist(strsplit(x, "[,;\\s]+"))))
42 }
43 }
44
45 # Convert only the 'sample' element in args
46 args$sample <- convertStringToNumeric(args$sample)
47
48 print("Argument types:")
49 print(sapply(args, class))
50
51 # Check if the image file exists
52 if (!file.exists(args$image)) {
53 stop("The RData file does not exist: ", args$image)
54 }
55
56 # ----- PROCESSING INFILE -----
57
58 # Load the RData file (it should contain the xset object, typically an xcmsSet or XCMSnExp)
59 load(args$image)
60 args$image <- NULL
61
62 # Save arguments to generate a report
63 if (!exists("listOFlistArguments")) listOFlistArguments <- list()
64 listOFlistArguments[[format(Sys.time(), "%y%m%d-%H:%M:%S_groupFWHM")]] <- args
65
66 # We unzip automatically the chromatograms from the zip files.
67 if (!exists("zipfile")) zipfile <- NULL
68 if (!exists("singlefile")) singlefile <- NULL
69 rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args)
70 zipfile <- rawFilePath$zipfile
71 singlefile <- rawFilePath$singlefile
72 args <- rawFilePath$args
73
74 print(paste("singlefile :", singlefile))
75 if (!is.null(singlefile)) {
76 directory <- retrieveRawfileInTheWorkingDir(singlefile, zipfile)
77 }
78
79 # If the xdata object exists, convert it to xcmsSet
80 if (exists("xdata")) {
81 phenoData <- xdata@phenoData
82 xset <- getxcmsSetObject(xdata)
83 }
84
85 if (!exists("xdata")) stop("\n\nERROR: The RData doesn't contain any object called 'xdata'. This RData should have been created by an old version of XMCS 2.*")
86
87 # Verification of a group step before doing the fillpeaks job.
88 if (dim(xdata@phenoData@data)[1] > 1) {
89 if (!hasFeatures(xdata)) stop("You must always do a group step after a retcor. Otherwise it won't work for the groupFWHM step")
90 } else {
91 print("Only one file in the phenoData keep xset as is")
92 }
93
94 # Convert the xset object to xsAnnotate using CAMERA
95 cat("Converting xset object to xsAnnotate...\n")
96 xsa <- xsAnnotate(xset, sample = args$sample, nSlaves = as.numeric(args$nSlaves), polarity = args$polarity)
97
98 print(paste0("All samples in xset object: ", paste(seq_along(xset@filepaths), collapse = ", ")))
99 print(paste0("Selected samples: ", paste(xsa@sample, collapse = ", ")))
100 print(paste0("Run in parallel mode (0 = disabled): ", paste(xsa@runParallel)))
101 print(paste0("Polarity: ", xsa@polarity))
102
103 # Apply the groupFWHM function with the parameters
104 cat("Applying groupFWHM...\n")
105 xa <- groupFWHM(xsa, sigma = as.numeric(args$sigma), perfwhm = as.numeric(args$perfwhm), intval = args$intval)
106
107 # Extract the list of annotated peaks
108 peakList <- getPeaklist(xa, intval = args$intval)
109
110 if (length(phenoData@data$sample_name) == 1) {
111 peakList$name <- make.unique(paste0("M", round(peakList[, "mz"], 0), "T", round(peakList[, "rt"], 0)), "_")
112 variableMetadata <- peakList[, c("name", setdiff(names(peakList), "name"))]
113 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
114 } else {
115 names_default <- groupnames(xa@xcmsSet, mzdec = 0, rtdec = 0) # Names without decimals
116 names_custom <- groupnames(xa@xcmsSet, mzdec = args$numDigitsMZ, rtdec = args$numDigitsRT) # Names with "x" decimals
117
118 variableMetadata <- data.frame(
119 name = names_default,
120 name_custom = names_custom,
121 stringsAsFactors = FALSE
122 )
123 variableMetadata <- cbind(variableMetadata, peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))])
124 }
125
126 if (!exists("RTinMinute")) RTinMinute <- FALSE
127
128 if (args$convertRTMinute && RTinMinute == FALSE) {
129 RTinMinute <- TRUE
130 variableMetadata <- RTSecondToMinute(variableMetadata = variableMetadata, convertRTMinute = args$convertRTMinute)
131 }
132
133 # Save the peak list to a TSV file
134 output_file_tsv <- "variableMetadata.tsv"
135 write.table(variableMetadata, file = output_file_tsv, sep = "\t", row.names = FALSE, quote = FALSE)
136
137 # Save the xsAnnotate object
138 output_file_RData <- "camera_fwhm.RData"
139 objects2save <- c("xa", "variableMetadata", "listOFlistArguments", "zipfile", "singlefile", "RTinMinute", "phenoData")
140 save(list = objects2save[objects2save %in% ls()], file = output_file_RData)
141
142 cat("Output files generated:", output_file_tsv, "and", output_file_RData, "\n")