Mercurial > repos > workflow4metabolomics > camera_groupcorr

comparison CAMERA_groupCorr.R @ 0:cb57be5de070 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:29:48 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:cb57be5de070
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 "NULL" strings to actual NULL values
29 convertNullString <- function(x) {
30 if (x == "NULL") {
31 return(NULL)
32 }
33 return(x)
34 }
35
36 # Function to convert string to numeric lists
37 convert_psg_list <- function(x) {
38 # Check if x is NULL or has zero length before further processing
39 if (is.null(x) || length(x) == 0) {
40 return(NULL)
41 }
42
43 # Force conversion to character
44 x <- as.character(x)
45
46 if (grepl("^[0-9]+$", x)) {
47 # If the string represents a single numeric value
48 return(as.numeric(x))
49 } else {
50 # Convert string representation of a list to a numeric vector
51 # Use a regular expression to split by common separators
52 return(as.numeric(unlist(strsplit(x, "[,;\\s]+"))))
53 }
54 }
55
56 for (arg in names(args)) {
57 args[[arg]] <- convertNullString(args[[arg]])
58 }
59
60 # Convert only the 'psg_list' element in args
61 args$psg_list <- convert_psg_list(args$psg_list)
62
63 print("Argument types:")
64 print(sapply(args, class))
65
66 # Check if the image file exists
67 if (!file.exists(args$image)) {
68 stop("The RData file does not exist: ", args$image)
69 }
70
71 # ----- PROCESSING INFILE -----
72
73 # Load the RData file
74 load(args$image)
75 args$image <- NULL
76
77 # Save arguments to generate a report
78 if (!exists("listOFlistArguments")) listOFlistArguments <- list()
79 listOFlistArguments[[format(Sys.time(), "%y%m%d-%H:%M:%S_groupCorr")]] <- args
80
81 # We unzip automatically the chromatograms from the zip files.
82 if (!exists("zipfile")) zipfile <- NULL
83 if (!exists("singlefile")) singlefile <- NULL
84 rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args)
85 zipfile <- rawFilePath$zipfile
86 singlefile <- rawFilePath$singlefile
87 args <- rawFilePath$args
88
89 print(paste("singlefile :", singlefile))
90 if (!is.null(singlefile)) {
91 directory <- retrieveRawfileInTheWorkingDir(singlefile, zipfile)
92 }
93
94 # Ensure the xa object is loaded
95 if (!exists("xa")) {
96 stop("The xa object was not found in the RData file.")
97 }
98
99 print("xa object loaded:")
100 print(xa)
101
102 # Apply the groupCorr function to the xsAnnotate object
103 print("Calling groupCorr function:")
104 xa <- groupCorr(xa, cor_eic_th = args$cor_eic_th, pval = args$pval, graphMethod = args$graphMethod, calcIso = args$calcIso, calcCiS = args$calcCiS, calcCaS = args$calcCaS, psg_list = args$psg_list, cor_exp_th = args$cor_exp_th, intval = args$intval)
105
106 print("Result of groupCorr function:")
107 print(xa)
108
109 # Extract the list of annotated peaks
110 peakList <- getPeaklist(xa, intval = args$intval)
111
112 if (length(phenoData@data$sample_name) == 1) {
113 peakList$name <- make.unique(paste0("M", round(peakList[, "mz"], 0), "T", round(peakList[, "rt"], 0)), "_")
114 variableMetadata <- peakList[, c("name", setdiff(names(peakList), "name"))]
115 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
116 } else {
117 names_default <- groupnames(xa@xcmsSet, mzdec = 0, rtdec = 0) # Names without decimals
118 names_custom <- groupnames(xa@xcmsSet, mzdec = args$numDigitsMZ, rtdec = args$numDigitsRT) # Names with "x" decimals
119
120 variableMetadata <- data.frame(
121 name = names_default,
122 name_custom = names_custom,
123 stringsAsFactors = FALSE
124 )
125 variableMetadata <- cbind(variableMetadata, peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))])
126 }
127
128 if (!exists("RTinMinute")) RTinMinute <- FALSE
129
130 if (args$convertRTMinute == TRUE && RTinMinute == FALSE) {
131 RTinMinute <- TRUE
132 variableMetadata <- RTSecondToMinute(variableMetadata = variableMetadata, convertRTMinute = args$convertRTMinute)
133 }
134
135 # Save the extracted peak list as a TSV file named 'variableMetadata.tsv'
136 output_file_tsv <- "variableMetadata.tsv"
137 write.table(variableMetadata, file = output_file_tsv, sep = "\t", row.names = FALSE, quote = FALSE)
138
139 # Save the updated xsAnnotate object
140 output_file_RData <- "camera_groupCorr.RData"
141 objects2save <- c("xa", "variableMetadata", "listOFlistArguments", "zipfile", "singlefile", "RTinMinute", "phenoData")
142 save(list = objects2save[objects2save %in% ls()], file = output_file_RData)
143
144 cat("Output files generated:", output_file_tsv, "and", output_file_RData, "\n")