Mercurial > repos > workflow4metabolomics > camera_groupfwhm
diff CAMERA_groupFWHM.R @ 0:6aea0427511e draft default tip
planemo upload commit 24d44ee26b7c23380c2b42fae2f7f6e58472100d
author | workflow4metabolomics |
---|---|
date | Sun, 24 Nov 2024 21:28:57 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CAMERA_groupFWHM.R Sun Nov 24 21:28:57 2024 +0000 @@ -0,0 +1,142 @@ +#!/usr/bin/env Rscript + +# ----- PACKAGE ----- +cat("\tSESSION INFO\n") + +# Import the different functions +source_local <- function(fname) { + argv <- commandArgs(trailingOnly = FALSE) + base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) + source(paste(base_dir, fname, sep = "/")) +} +source_local("lib.r") + +pkgs <- c("CAMERA", "xcms", "multtest", "batch") +loadAndDisplayPackages(pkgs) +cat("\n\n") +# ----- ARGUMENTS ----- +cat("\tARGUMENTS INFO\n") + +args <- parseCommandArgs(evaluate = FALSE) # interpretation of arguments given in command line as an R list of objects +write.table(as.matrix(args), col.names = FALSE, quote = FALSE, sep = "\t") + +cat("\n\n") + +print("Arguments retrieved from the command line:") +print(args) + +# Function to convert "NA" strings to actual NA values and string lists to numeric lists +convertStringToNumeric <- function(x) { + # Force conversion to character + x <- as.character(x) + + if (x == "NA") { + return(NA) + } else if (grepl("^[0-9]+$", x)) { + # If the string represents a single numeric value + return(as.numeric(x)) + } else { + # Convert string representation of a list to a numeric vector + # Use a regular expression to split by common separators + return(as.numeric(unlist(strsplit(x, "[,;\\s]+")))) + } +} + +# Convert only the 'sample' element in args +args$sample <- convertStringToNumeric(args$sample) + +print("Argument types:") +print(sapply(args, class)) + +# Check if the image file exists +if (!file.exists(args$image)) { + stop("The RData file does not exist: ", args$image) +} + +# ----- PROCESSING INFILE ----- + +# Load the RData file (it should contain the xset object, typically an xcmsSet or XCMSnExp) +load(args$image) +args$image <- NULL + +# Save arguments to generate a report +if (!exists("listOFlistArguments")) listOFlistArguments <- list() +listOFlistArguments[[format(Sys.time(), "%y%m%d-%H:%M:%S_groupFWHM")]] <- args + +# We unzip automatically the chromatograms from the zip files. +if (!exists("zipfile")) zipfile <- NULL +if (!exists("singlefile")) singlefile <- NULL +rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args) +zipfile <- rawFilePath$zipfile +singlefile <- rawFilePath$singlefile +args <- rawFilePath$args + +print(paste("singlefile :", singlefile)) +if (!is.null(singlefile)) { + directory <- retrieveRawfileInTheWorkingDir(singlefile, zipfile) +} + +# If the xdata object exists, convert it to xcmsSet +if (exists("xdata")) { + phenoData <- xdata@phenoData + xset <- getxcmsSetObject(xdata) +} + +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.*") + +# Verification of a group step before doing the fillpeaks job. +if (dim(xdata@phenoData@data)[1] > 1) { + if (!hasFeatures(xdata)) stop("You must always do a group step after a retcor. Otherwise it won't work for the groupFWHM step") +} else { + print("Only one file in the phenoData keep xset as is") +} + +# Convert the xset object to xsAnnotate using CAMERA +cat("Converting xset object to xsAnnotate...\n") +xsa <- xsAnnotate(xset, sample = args$sample, nSlaves = as.numeric(args$nSlaves), polarity = args$polarity) + +print(paste0("All samples in xset object: ", paste(seq_along(xset@filepaths), collapse = ", "))) +print(paste0("Selected samples: ", paste(xsa@sample, collapse = ", "))) +print(paste0("Run in parallel mode (0 = disabled): ", paste(xsa@runParallel))) +print(paste0("Polarity: ", xsa@polarity)) + +# Apply the groupFWHM function with the parameters +cat("Applying groupFWHM...\n") +xa <- groupFWHM(xsa, sigma = as.numeric(args$sigma), perfwhm = as.numeric(args$perfwhm), intval = args$intval) + +# Extract the list of annotated peaks +peakList <- getPeaklist(xa, intval = args$intval) + +if (length(phenoData@data$sample_name) == 1) { + peakList$name <- make.unique(paste0("M", round(peakList[, "mz"], 0), "T", round(peakList[, "rt"], 0)), "_") + variableMetadata <- peakList[, c("name", setdiff(names(peakList), "name"))] + variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ) +} else { + names_default <- groupnames(xa@xcmsSet, mzdec = 0, rtdec = 0) # Names without decimals + names_custom <- groupnames(xa@xcmsSet, mzdec = args$numDigitsMZ, rtdec = args$numDigitsRT) # Names with "x" decimals + + variableMetadata <- data.frame( + name = names_default, + name_custom = names_custom, + stringsAsFactors = FALSE + ) + variableMetadata <- cbind(variableMetadata, peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]) +} + +if (!exists("RTinMinute")) RTinMinute <- FALSE + +if (args$convertRTMinute && RTinMinute == FALSE) { + RTinMinute <- TRUE + variableMetadata <- RTSecondToMinute(variableMetadata = variableMetadata, convertRTMinute = args$convertRTMinute) +} + +# Save the peak list to a TSV file +output_file_tsv <- "variableMetadata.tsv" +write.table(variableMetadata, file = output_file_tsv, sep = "\t", row.names = FALSE, quote = FALSE) + +# Save the xsAnnotate object +output_file_RData <- "camera_fwhm.RData" +objects2save <- c("xa", "variableMetadata", "listOFlistArguments", "zipfile", "singlefile", "RTinMinute", "phenoData") +save(list = objects2save[objects2save %in% ls()], file = output_file_RData) + +cat("Output files generated:", output_file_tsv, "and", output_file_RData, "\n")