Mercurial > repos > workflow4metabolomics > xcms_refine
changeset 0:eb115eb8f25c draft
"planemo upload for repository https://github.com/workflow4metabolomics/xcms commit d4bb3c31e2beca6b4059758c9e76d42c26f8aa55"
| author | workflow4metabolomics |
|---|---|
| date | Thu, 17 Jun 2021 07:42:55 +0000 |
| parents | |
| children | 20f8ebc3a391 |
| files | README.rst abims_xcms_refine.xml lib.r macros.xml macros_xcms.xml repository_dependencies.xml test-data/faahKO-single.xset.merged.RData test-data/faahKO_reduce.zip test-data/ko15.CDF test-data/ko16.CDF test-data/wt15.CDF test-data/wt16.CDF xcms_refine.r |
| diffstat | 13 files changed, 1180 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.rst Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,10 @@ + +Changelog/News +-------------- + +.. _News: https://bioconductor.org/packages/release/bioc/news/xcms/NEWS + +**Version 3.12.0.0 - 08/02/2021** + +- NEW: first version of tool wrapper +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/abims_xcms_refine.xml Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,264 @@ +<tool id="abims_xcms_refine" name="xcms refineChromPeaks (refine)" version="3.12.0+galaxy0"> + + <description>Remove or merge chromatographic peaks based on specific criteria.</description> + + <macros> + <import>macros.xml</import> + <import>macros_xcms.xml</import> + </macros> + + <expand macro="requirements"/> + + <expand macro="stdio"/> + + <command><![CDATA[ + @COMMAND_RSCRIPT@/xcms_refine.r + image '$image' + + method $methods.method + #if $methods.method == "CleanPeaks": + maxPeakwidth $methods.maxPeakwidth + #elif $methods.method == "FilterIntensity": + threshold $methods.threshold + nValues $methods.nValues + value $methods.value + #else: + expandRt $methods.expandRt + expandMz $methods.expandMz + ppm $methods.ppm + minProp $methods.minProp + #end if + + msLevel $msLevel + + @COMMAND_FILE_LOAD@ + + @COMMAND_LOG_EXIT@ + ]]></command> + + <inputs> + <param name="image" type="data" format="rdata.xcms.findchrompeaks" label="@INPUT_IMAGE_LABEL@" help="@INPUT_IMAGE_HELP@ from: findChromPeaks" /> + <conditional name="methods"> + <param name="method" type="select" label="Method for peak refinement" help="See the help section below"> + <option value="CleanPeaks" selected="true">CleanPeaks - remove peaks that are too wide</option> + <option value="FilterIntensity">FilterIntensity - remove peaks with too low intensity</option> + <option value="MergeNeighboringPeaks">MergeNeighboringPeaks - combine peaks that are too close together</option> + </param> + <when value="CleanPeaks"> + <param argument="maxPeakwidth" type="float" value="10" label="MaxPeakWidth" help="maximum peak width in seconds. See CleanPeaks manual (linked in help below)." /> + </when> + <when value="FilterIntensity"> + <param argument="threshold" type="float" value="0" label="minimal required intensity" help="minimal required intensity for a peak to be retained. Defaults to 0.0." /> + <param argument="nValues" type="integer" value="1" label="minimum number of data points per peak" help="number of data points (per chromatographic peak) that have to be >= threshold. Defaults to 1." /> + <param name="value" type="select" label="Definition of peak intensity" help="See the FilterIntensity manual (linked in help below). 'intb' requires centWave processed data"> + <option value="maxo" selected="true">maxo - use maximum peak intensities</option> + <option value="into">into - use integrated peak intensities</option> + <option value="intb">intb - use baseline-corrected integrated peak intensities</option> + </param> + </when> + <when value="MergeNeighboringPeaks"> + <param argument="expandRt" type="float" value="2" label="rt expansion (seconds)" help="seconds of expansion of retention time window (on both sides) to check for overlapping peaks. Defaults to 2." /> + <param argument="expandMz" type="float" value="0" label="m/z expansion (m/z units)" help="'m/z units' of expansion of m/z window (on both sides) to check for overlapping peaks. Defaults to 0." /> + <param argument="ppm" type="float" value="10" label="ppm expansion (ppm units)" help="'relative 'm/z units' (in parts per million) of expansion of m/z window (on both sides) to check for overlapping peaks. Defaults to 10." /> + <param argument="minProp" type="float" value="0.75" label="minimum proportion at midpoint" help="See MergeNeighboringPeaks manual (linked in help below). Minimum proportion midway between peaks relative to the smaller peak's 'maxo' (maximal intensity at peak apex). Defaults to 0.75." /> + </when> + </conditional> + <param argument="msLevel" type="integer" value="1" label="msLevel for refinement" help="the MS level on which refinement should be performed. Defaults to 1." /> + + <expand macro="input_file_load"/> + </inputs> + + <outputs> + <data name="xsetRData" format="rdata.xcms.findchrompeaks" label="${image.name.rsplit('.',1)[0]}.refine.RData" from_work_dir="xcmsSet.RData" /> + </outputs> + + <tests> + <!-- + test 1 - MergeNeighboringPeaks method + --> + <test> + <param name="image" value="faahKO-single.xset.merged.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="MergeNeighboringPeaks"/> + <param name="expandRt" value="2" /> + <param name="expandMz" value="0" /> + <param name="ppm" value="10" /> + <param name="minProp" value="0.75" /> + </conditional> + <expand macro="test_file_load_single"/> + <assert_stderr> + <has_text text="Evaluating 547 peaks in file wt15.CDF for merging ... OK" /> + <has_text text="Evaluating 547 peaks in file wt16.CDF for merging ... OK" /> + <has_text text="Evaluating 558 peaks in file ko15.CDF for merging ... OK"/> + <has_text text="Evaluating 589 peaks in file ko16.CDF for merging ... OK" /> + <has_text text="Merging reduced 9251 chromPeaks to 9093" /> + </assert_stderr> + <assert_stdout> + <has_text text="object with 4 samples" /> + <has_text text="Time range: 2506.1-4477.9 seconds (41.8-74.6 minutes)" /> + <has_text text="Mass range: 200.1-600 m/z" /> + <has_text text="Peaks: 9093 (about 2273 per sample)" /> + </assert_stdout> + </test> + <!-- + test 2 - FilterIntensity method + --> + <test> + <param name="image" value="faahKO-single.xset.merged.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="FilterIntensity"/> + <param name="threshold" value="1e5" /> + <param name="nValues" value="3" /> + <param name="value" value="into" /> + </conditional> + <expand macro="test_file_load_single"/> + <assert_stderr> + <has_text text="Removed 9018 chromatographic peaks" /> + </assert_stderr> + <assert_stdout> + <has_text text="object with 4 samples" /> + <has_text text="Time range: 2657.9-4124.2 seconds (44.3-68.7 minutes)" /> + <has_text text="Mass range: 256.1-599.4 m/z" /> + <has_text text="Peaks: 233 (about 58 per sample)" /> + </assert_stdout> + </test> + <!-- + test 3 - CleanPeaks method + --> + <test> + <param name="image" value="faahKO-single.xset.merged.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="CleanPeaks"/> + <param name="maxPeakwidth" value="8" /> + </conditional> + <expand macro="test_file_load_single"/> + <assert_stderr> + <has_text text="Removed 6037 of 9251 chromatographic peaks" /> + </assert_stderr> + <assert_stdout> + <has_text text="object with 4 samples" /> + <has_text text="Time range: 2524.9-4473.2 seconds (42.1-74.6 minutes)" /> + <has_text text="Mass range: 200.2-600 m/z" /> + <has_text text="Peaks: 3214 (about 804 per sample)" /> + </assert_stdout> + </test> + </tests> + + <help><![CDATA[ + +.. class:: infomark + +**Authors** Johannes Rainer Johannes.Rainer@eurac.edu and Mar Garcia-Aloy mar.garcia@fmach.it + +@HELP_AUTHORS_WRAPPERS@ + +--------------------------------------------------- + + +===================== +xcms refineChromPeaks +===================== + +----------- +Description +----------- + +After peak identification with xcms findChromPeaks (xcmsSet), this tool refines those peaks. +It either removes peaks that are too wide or removes peaks with too low intensity or combines peaks that are too close together. + +Note well that refineChromPeaks methods will always remove feature definitions, +because a call to this method can change or remove identified chromatographic peaks, which may be part of features. +Therefore it must only be run immediately after findChromPeaks (xcmsSet). + +----------------- +Workflow position +----------------- + +**Upstream tools** + +==================================== ======================== ============================== +Name Output file Format +==================================== ======================== ============================== +xcms.findChromPeaks Merger (single) xset.merged.RData rdata.xcms.findchrompeaks +------------------------------------ ------------------------ ------------------------------ +xcms.findChromPeaks (zip) ``*``.raw.xset.RData rdata.xcms.findchrompeaks +==================================== ======================== ============================== + + +**Downstream tools** + +==================================== ======================== ========================= +Name Output file Format +==================================== ======================== ========================= +xcms.findChromPeaks Merger (single) ``*``.raw.xset.RData rdata.xcms.findchrompeaks +------------------------------------ ------------------------ ------------------------- +xcms.groupChromPeaks (zip) ``*``.raw.xset.RData rdata.xcms.findchrompeaks +==================================== ======================== ========================= + +**General schema of the metabolomic workflow** + +This tool would appear between findChromPeaks and group, but only when the choice of peaks needs to be refined, e.g., to accommodate HILIC data. + +--------------------------------------------------- + +---------- +Parameters +---------- + +Method to use for refinement +---------------------------- + +**CleanPeaks** + + | Remove chromatographic peaks with too large rt width. + | See the CleanPeaks_manual_ + +**FilterIntensity** + + | Remove chromatographic peaks based on intensity. + | See the FilterIntensity_manual_ + +**MergeNeighboringPeaks** + + | Merge neighboring and overlapping chromatographic peaks. + | See the MergeNeighboringPeaks_manual_ + +.. _CleanPeaks_manual: https://rdrr.io/bioc/xcms/man/refineChromPeaks-clean.html +.. _FilterIntensity_manual: https://rdrr.io/bioc/xcms/man/refineChromPeaks-filter-intensity.html +.. _MergeNeighboringPeaks_manual: https://rdrr.io/bioc/xcms/man/refineChromPeaks-merge.html + +@HELP_XCMS_MANUAL@ + +------------ +Output files +------------ + +xset.RData: rdata.xcms.findchrompeaks format + + | (single) RData files that are necessary in the second step of the workflow "xcms.groupChromPeaks" - must be merged first using "xcms.findChromPeaks Merger" + | (zip) RData file that is necessary in the second step of the workflow "xcms.groupChromPeaks". + + +--------------------------------------------------- + +Changelog/News +-------------- + +.. _News: https://bioconductor.org/packages/release/bioc/news/xcms/NEWS + + +**Version 3.6.1.0 - 08/02/2021** + +- NEW: first version of tool wrapper + + + ]]></help> + + <citations> + <citation type="doi">10.5281/zenodo.3909299</citation> + <expand macro="citation_w4m"/> + <citation type="doi">10.1021/ac051437y</citation> + </citations> + + +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib.r Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,451 @@ +#@authors ABiMS TEAM, Y. Guitton +# lib.r for Galaxy Workflow4Metabolomics xcms tools + +#@author G. Le Corguille +# solve an issue with batch if arguments are logical TRUE/FALSE +parseCommandArgs <- function(...) { + args <- batch::parseCommandArgs(...) + for (key in names(args)) { + if (args[key] %in% c("TRUE", "FALSE")) + args[key] <- as.logical(args[key]) + } + return(args) +} + +#@author G. Le Corguille +# This function will +# - load the packages +# - display the sessionInfo +loadAndDisplayPackages <- function(pkgs) { + for (pkg in pkgs) suppressPackageStartupMessages(stopifnot(library(pkg, quietly = TRUE, logical.return = TRUE, character.only = TRUE))) + + sessioninfo <- sessionInfo() + cat(sessioninfo$R.version$version.string, "\n") + cat("Main packages:\n") + for (pkg in names(sessioninfo$otherPkgs)) { + cat(paste(pkg, packageVersion(pkg)), "\t") + } + cat("\n") + cat("Other loaded packages:\n") + for (pkg in names(sessioninfo$loadedOnly)) { + cat(paste(pkg, packageVersion(pkg)), "\t") + } + cat("\n") +} + +#@author G. Le Corguille +# This function merge several chromBPI or chromTIC into one. +mergeChrom <- function(chrom_merged, chrom) { + if (is.null(chrom_merged)) return(NULL) + chrom_merged@.Data <- cbind(chrom_merged@.Data, chrom@.Data) + return(chrom_merged) +} + +#@author G. Le Corguille +# This function merge several xdata into one. +mergeXData <- function(args) { + chromTIC <- NULL + chromBPI <- NULL + chromTIC_adjusted <- NULL + chromBPI_adjusted <- NULL + md5sumList <- NULL + for (image in args$images) { + + load(image) + # Handle infiles + if (!exists("singlefile")) singlefile <- NULL + if (!exists("zipfile")) zipfile <- NULL + rawFilePath <- retrieveRawfileInTheWorkingDir(singlefile, zipfile, args) + zipfile <- rawFilePath$zipfile + singlefile <- rawFilePath$singlefile + + if (exists("raw_data")) xdata <- raw_data + 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.*") + + cat(sampleNamesList$sampleNamesOrigin, "\n") + + if (!exists("xdata_merged")) { + xdata_merged <- xdata + singlefile_merged <- singlefile + md5sumList_merged <- md5sumList + sampleNamesList_merged <- sampleNamesList + chromTIC_merged <- chromTIC + chromBPI_merged <- chromBPI + chromTIC_adjusted_merged <- chromTIC_adjusted + chromBPI_adjusted_merged <- chromBPI_adjusted + } else { + if (is(xdata, "XCMSnExp")) xdata_merged <- c(xdata_merged, xdata) + else if (is(xdata, "OnDiskMSnExp")) xdata_merged <- xcms:::.concatenate_OnDiskMSnExp(xdata_merged, xdata) + else stop("\n\nERROR: The RData either a OnDiskMSnExp object called raw_data or a XCMSnExp object called xdata") + + singlefile_merged <- c(singlefile_merged, singlefile) + md5sumList_merged$origin <- rbind(md5sumList_merged$origin, md5sumList$origin) + sampleNamesList_merged$sampleNamesOrigin <- c(sampleNamesList_merged$sampleNamesOrigin, sampleNamesList$sampleNamesOrigin) + sampleNamesList_merged$sampleNamesMakeNames <- c(sampleNamesList_merged$sampleNamesMakeNames, sampleNamesList$sampleNamesMakeNames) + chromTIC_merged <- mergeChrom(chromTIC_merged, chromTIC) + chromBPI_merged <- mergeChrom(chromBPI_merged, chromBPI) + chromTIC_adjusted_merged <- mergeChrom(chromTIC_adjusted_merged, chromTIC_adjusted) + chromBPI_adjusted_merged <- mergeChrom(chromBPI_adjusted_merged, chromBPI_adjusted) + } + } + rm(image) + xdata <- xdata_merged; rm(xdata_merged) + singlefile <- singlefile_merged; rm(singlefile_merged) + md5sumList <- md5sumList_merged; rm(md5sumList_merged) + sampleNamesList <- sampleNamesList_merged; rm(sampleNamesList_merged) + + if (!is.null(args$sampleMetadata)) { + cat("\tXSET PHENODATA SETTING...\n") + sampleMetadataFile <- args$sampleMetadata + sampleMetadata <- getDataFrameFromFile(sampleMetadataFile, header = F) + xdata@phenoData@data$sample_group <- sampleMetadata$V2[match(xdata@phenoData@data$sample_name, sampleMetadata$V1)] + + if (any(is.na(pData(xdata)$sample_group))) { + sample_missing <- pData(xdata)$sample_name[is.na(pData(xdata)$sample_group)] + error_message <- paste("Those samples are missing in your sampleMetadata:", paste(sample_missing, collapse = " ")) + print(error_message) + stop(error_message) + } + } + + if (!is.null(chromTIC_merged)) { + chromTIC <- chromTIC_merged; chromTIC@phenoData <- xdata@phenoData + } + if (!is.null(chromBPI_merged)) { + chromBPI <- chromBPI_merged; chromBPI@phenoData <- xdata@phenoData + } + if (!is.null(chromTIC_adjusted_merged)) { + chromTIC_adjusted <- chromTIC_adjusted_merged; chromTIC_adjusted@phenoData <- xdata@phenoData + } + if (!is.null(chromBPI_adjusted_merged)) { + chromBPI_adjusted <- chromBPI_adjusted_merged; chromBPI_adjusted@phenoData <- xdata@phenoData + } + + return(list("xdata" = xdata, "singlefile" = singlefile, "md5sumList" = md5sumList, "sampleNamesList" = sampleNamesList, "chromTIC" = chromTIC, "chromBPI" = chromBPI, "chromTIC_adjusted" = chromTIC_adjusted, "chromBPI_adjusted" = chromBPI_adjusted)) +} + +#@author G. Le Corguille +# This function convert if it is required the Retention Time in minutes +RTSecondToMinute <- function(variableMetadata, convertRTMinute) { + if (convertRTMinute) { + #converting the retention times (seconds) into minutes + print("converting the retention times into minutes in the variableMetadata") + variableMetadata[, "rt"] <- variableMetadata[, "rt"] / 60 + variableMetadata[, "rtmin"] <- variableMetadata[, "rtmin"] / 60 + variableMetadata[, "rtmax"] <- variableMetadata[, "rtmax"] / 60 + } + return(variableMetadata) +} + +#@author G. Le Corguille +# This function format ions identifiers +formatIonIdentifiers <- function(variableMetadata, numDigitsRT = 0, numDigitsMZ = 0) { + splitDeco <- strsplit(as.character(variableMetadata$name), "_") + idsDeco <- sapply(splitDeco, + function(x) { + deco <- unlist(x)[2]; if (is.na(deco)) return("") else return(paste0("_", deco)) + } + ) + namecustom <- make.unique(paste0("M", round(variableMetadata[, "mz"], numDigitsMZ), "T", round(variableMetadata[, "rt"], numDigitsRT), idsDeco)) + variableMetadata <- cbind(name = variableMetadata$name, namecustom = namecustom, variableMetadata[, !(colnames(variableMetadata) %in% c("name"))]) + return(variableMetadata) +} + +#@author G. Le Corguille +# This function convert the remain NA to 0 in the dataMatrix +naTOzeroDataMatrix <- function(dataMatrix, naTOzero) { + if (naTOzero) { + dataMatrix[is.na(dataMatrix)] <- 0 + } + return(dataMatrix) +} + +#@author G. Le Corguille +# Draw the plotChromPeakDensity 3 per page in a pdf file +getPlotChromPeakDensity <- function(xdata, param = NULL, mzdigit = 4) { + pdf(file = "plotChromPeakDensity.pdf", width = 16, height = 12) + + par(mfrow = c(3, 1), mar = c(4, 4, 1, 0.5)) + + if (length(unique(xdata$sample_group)) < 10) { + group_colors <- brewer.pal(length(unique(xdata$sample_group)), "Set1") + }else{ + group_colors <- hcl.colors(length(unique(xdata$sample_group)), palette = "Dark 3") + } + names(group_colors) <- unique(xdata$sample_group) + col_per_samp <- as.character(xdata$sample_group) + for (i in seq_len(length(group_colors))) { + col_per_samp[col_per_samp == (names(group_colors)[i])] <- group_colors[i] + } + + xlim <- c(min(featureDefinitions(xdata)$rtmin), max(featureDefinitions(xdata)$rtmax)) + for (i in seq_len(nrow(featureDefinitions(xdata)))) { + mzmin <- featureDefinitions(xdata)[i, ]$mzmin + mzmax <- featureDefinitions(xdata)[i, ]$mzmax + plotChromPeakDensity(xdata, param = param, mz = c(mzmin, mzmax), col = col_per_samp, pch = 16, xlim = xlim, main = paste(round(mzmin, mzdigit), round(mzmax, mzdigit))) + legend("topright", legend = names(group_colors), col = group_colors, cex = 0.8, lty = 1) + } + + dev.off() +} + +#@author G. Le Corguille +# Draw the plotChromPeakDensity 3 per page in a pdf file +getPlotAdjustedRtime <- function(xdata) { + + pdf(file = "raw_vs_adjusted_rt.pdf", width = 16, height = 12) + + # Color by group + if (length(unique(xdata$sample_group)) < 10) { + group_colors <- brewer.pal(length(unique(xdata$sample_group)), "Set1") + } else { + group_colors <- hcl.colors(length(unique(xdata$sample_group)), palette = "Dark 3") + } + if (length(group_colors) > 1) { + names(group_colors) <- unique(xdata$sample_group) + plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group]) + legend("topright", legend = names(group_colors), col = group_colors, cex = 0.8, lty = 1) + } + + # Color by sample + plotAdjustedRtime(xdata, col = rainbow(length(xdata@phenoData@data$sample_name))) + legend("topright", legend = xdata@phenoData@data$sample_name, col = rainbow(length(xdata@phenoData@data$sample_name)), cex = 0.8, lty = 1) + + dev.off() +} + +#@author G. Le Corguille +# value: intensity values to be used into, maxo or intb +getPeaklistW4M <- function(xdata, intval = "into", convertRTMinute = F, numDigitsMZ = 4, numDigitsRT = 0, naTOzero = T, variableMetadataOutput, dataMatrixOutput, sampleNamesList) { + dataMatrix <- featureValues(xdata, method = "medret", value = intval) + colnames(dataMatrix) <- make.names(tools::file_path_sans_ext(colnames(dataMatrix))) + dataMatrix <- cbind(name = groupnames(xdata), dataMatrix) + variableMetadata <- featureDefinitions(xdata) + colnames(variableMetadata)[1] <- "mz"; colnames(variableMetadata)[4] <- "rt" + variableMetadata <- data.frame(name = groupnames(xdata), variableMetadata) + + variableMetadata <- RTSecondToMinute(variableMetadata, convertRTMinute) + variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = numDigitsRT, numDigitsMZ = numDigitsMZ) + dataMatrix <- naTOzeroDataMatrix(dataMatrix, naTOzero) + + # FIX: issue when the vector at peakidx is too long and is written in a new line during the export + variableMetadata[, "peakidx"] <- vapply(variableMetadata[, "peakidx"], FUN = paste, FUN.VALUE = character(1), collapse = ",") + + write.table(variableMetadata, file = variableMetadataOutput, sep = "\t", quote = F, row.names = F) + write.table(dataMatrix, file = dataMatrixOutput, sep = "\t", quote = F, row.names = F) + +} + +#@author G. Le Corguille +# It allow different of field separators +getDataFrameFromFile <- function(filename, header = T) { + myDataFrame <- read.table(filename, header = header, sep = ";", stringsAsFactors = F) + if (ncol(myDataFrame) < 2) myDataFrame <- read.table(filename, header = header, sep = "\t", stringsAsFactors = F) + if (ncol(myDataFrame) < 2) myDataFrame <- read.table(filename, header = header, sep = ",", stringsAsFactors = F) + if (ncol(myDataFrame) < 2) { + error_message <- "Your tabular file seems not well formatted. The column separators accepted are ; , and tabulation" + print(error_message) + stop(error_message) + } + return(myDataFrame) +} + +#@author G. Le Corguille +# Draw the BPI and TIC graphics +# colored by sample names or class names +getPlotChromatogram <- function(chrom, xdata, pdfname = "Chromatogram.pdf", aggregationFun = "max") { + + if (aggregationFun == "sum") + type <- "Total Ion Chromatograms" + else + type <- "Base Peak Intensity Chromatograms" + + adjusted <- "Raw" + if (hasAdjustedRtime(xdata)) + adjusted <- "Adjusted" + + main <- paste(type, ":", adjusted, "data") + + pdf(pdfname, width = 16, height = 10) + + # Color by group + if (length(unique(xdata$sample_group)) < 10) { + group_colors <- brewer.pal(length(unique(xdata$sample_group)), "Set1") + }else{ + group_colors <- hcl.colors(length(unique(xdata$sample_group)), palette = "Dark 3") + } + if (length(group_colors) > 1) { + names(group_colors) <- unique(xdata$sample_group) + plot(chrom, col = group_colors[chrom$sample_group], main = main, peakType = "none") + legend("topright", legend = names(group_colors), col = group_colors, cex = 0.8, lty = 1) + } + + # Color by sample + plot(chrom, col = rainbow(length(xdata@phenoData@data$sample_name)), main = main, peakType = "none") + legend("topright", legend = xdata@phenoData@data$sample_name, col = rainbow(length(xdata@phenoData@data$sample_name)), cex = 0.8, lty = 1) + + dev.off() +} + + +# Get the polarities from all the samples of a condition +#@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM +getSampleMetadata <- function(xdata = NULL, sampleMetadataOutput = "sampleMetadata.tsv") { + cat("Creating the sampleMetadata file...\n") + + #Create the sampleMetada dataframe + sampleMetadata <- xdata@phenoData@data + rownames(sampleMetadata) <- NULL + colnames(sampleMetadata) <- c("sample_name", "class") + + sampleNamesOrigin <- sampleMetadata$sample_name + sampleNamesMakeNames <- make.names(sampleNamesOrigin) + + if (any(duplicated(sampleNamesMakeNames))) { + write("\n\nERROR: Usually, R has trouble to deal with special characters in its column names, so it rename them using make.names().\nIn your case, at least two columns after the renaming obtain the same name, thus XCMS will collapse those columns per name.", stderr()) + for (sampleName in sampleNamesOrigin) { + write(paste(sampleName, "\t->\t", make.names(sampleName)), stderr()) + } + stop("\n\nERROR: One or more of your files will not be import by xcmsSet. It may due to bad characters in their filenames.") + } + + if (!all(sampleNamesOrigin == sampleNamesMakeNames)) { + cat("\n\nWARNING: Usually, R has trouble to deal with special characters in its column names, so it rename them using make.names()\nIn your case, one or more sample names will be renamed in the sampleMetadata and dataMatrix files:\n") + for (sampleName in sampleNamesOrigin) { + cat(paste(sampleName, "\t->\t", make.names(sampleName), "\n")) + } + } + + sampleMetadata$sample_name <- sampleNamesMakeNames + + + #For each sample file, the following actions are done + for (fileIdx in seq_len(length(fileNames(xdata)))) { + #Check if the file is in the CDF format + if (!mzR:::netCDFIsFile(fileNames(xdata))) { + + # If the column isn't exist, with add one filled with NA + if (is.null(sampleMetadata$polarity)) sampleMetadata$polarity <- NA + + #Extract the polarity (a list of polarities) + polarity <- fData(xdata)[fData(xdata)$fileIdx == fileIdx, "polarity"] + #Verify if all the scans have the same polarity + uniq_list <- unique(polarity) + if (length(uniq_list) > 1) { + polarity <- "mixed" + } else { + polarity <- as.character(uniq_list) + } + + #Set the polarity attribute + sampleMetadata$polarity[fileIdx] <- polarity + } + + } + + write.table(sampleMetadata, sep = "\t", quote = FALSE, row.names = FALSE, file = sampleMetadataOutput) + + return(list("sampleNamesOrigin" = sampleNamesOrigin, "sampleNamesMakeNames" = sampleNamesMakeNames)) + +} + + +# This function will compute MD5 checksum to check the data integrity +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +getMd5sum <- function(files) { + cat("Compute md5 checksum...\n") + library(tools) + return(as.matrix(md5sum(files))) +} + +# This function retrieve the raw file in the working directory +# - if zipfile: unzip the file with its directory tree +# - if singlefiles: set symlink with the good filename +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +retrieveRawfileInTheWorkingDir <- function(singlefile, zipfile, args, prefix = "") { + + if (!(prefix %in% c("", "Positive", "Negative", "MS1", "MS2"))) stop("prefix must be either '', 'Positive', 'Negative', 'MS1' or 'MS2'") + + # single - if the file are passed in the command arguments -> refresh singlefile + if (!is.null(args[[paste0("singlefile_galaxyPath", prefix)]])) { + singlefile_galaxyPaths <- unlist(strsplit(args[[paste0("singlefile_galaxyPath", prefix)]], "\\|")) + singlefile_sampleNames <- unlist(strsplit(args[[paste0("singlefile_sampleName", prefix)]], "\\|")) + + singlefile <- NULL + for (singlefile_galaxyPath_i in seq_len(length(singlefile_galaxyPaths))) { + singlefile_galaxyPath <- singlefile_galaxyPaths[singlefile_galaxyPath_i] + singlefile_sampleName <- singlefile_sampleNames[singlefile_galaxyPath_i] + # In case, an url is used to import data within Galaxy + singlefile_sampleName <- tail(unlist(strsplit(singlefile_sampleName, "/")), n = 1) + singlefile[[singlefile_sampleName]] <- singlefile_galaxyPath + } + } + # zipfile - if the file are passed in the command arguments -> refresh zipfile + if (!is.null(args[[paste0("zipfile", prefix)]])) + zipfile <- args[[paste0("zipfile", prefix)]] + + # single + if (!is.null(singlefile) && (length("singlefile") > 0)) { + files <- vector() + for (singlefile_sampleName in names(singlefile)) { + singlefile_galaxyPath <- singlefile[[singlefile_sampleName]] + if (!file.exists(singlefile_galaxyPath)) { + error_message <- paste("Cannot access the sample:", singlefile_sampleName, "located:", singlefile_galaxyPath, ". Please, contact your administrator ... if you have one!") + print(error_message); stop(error_message) + } + + if (!suppressWarnings(try(file.link(singlefile_galaxyPath, singlefile_sampleName), silent = T))) + file.copy(singlefile_galaxyPath, singlefile_sampleName) + files <- c(files, singlefile_sampleName) + } + } + # zipfile + if (!is.null(zipfile) && (zipfile != "")) { + if (!file.exists(zipfile)) { + error_message <- paste("Cannot access the Zip file:", zipfile, ". Please, contact your administrator ... if you have one!") + print(error_message) + stop(error_message) + } + suppressWarnings(unzip(zipfile, unzip = "unzip")) + + #get the directory name + suppressWarnings(filesInZip <- unzip(zipfile, list = T)) + directories <- unique(unlist(lapply(strsplit(filesInZip$Name, "/"), function(x) x[1]))) + directories <- directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir] + directory <- "." + if (length(directories) == 1) directory <- directories + + cat("files_root_directory\t", directory, "\n") + + filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") + filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|") + info <- file.info(directory) + listed <- list.files(directory[info$isdir], pattern = filepattern, recursive = TRUE, full.names = TRUE) + files <- c(directory[!info$isdir], listed) + exists <- file.exists(files) + files <- files[exists] + + } + return(list(zipfile = zipfile, singlefile = singlefile, files = files)) + +} + + +# This function retrieve a xset like object +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +getxcmsSetObject <- function(xobject) { + # XCMS 1.x + if (class(xobject) == "xcmsSet") + return(xobject) + # XCMS 3.x + if (class(xobject) == "XCMSnExp") { + # Get the legacy xcmsSet object + suppressWarnings(xset <- as(xobject, "xcmsSet")) + if (!is.null(xset@phenoData$sample_group)) + sampclass(xset) <- xset@phenoData$sample_group + else + sampclass(xset) <- "." + return(xset) + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros.xml Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,52 @@ +<?xml version="1.0"?> +<macros> + <xml name="stdio"> + <stdio> + <exit_code range="1" level="fatal" /> + </stdio> + </xml> + + <!-- COMMAND --> + <token name="@COMMAND_RSCRIPT@">LC_ALL=C Rscript $__tool_directory__/</token> + + <token name="@COMMAND_LOG_EXIT@"> + ; + return=\$?; + cat 'log.txt'; + sh -c "exit \$return" + </token> + + <!-- INPUT_VALIDATORS --> + <xml name="input_validator_range_integer"> + <validator type="regex" message="The format is 'min,max'" >[0-9]+ *, *[0-9]+</validator> + </xml> + + <xml name="input_validator_range_float"> + <validator type="regex" message="The format is 'min,max'" >[0-9]+\.?[0-9]* *, *[0-9]+\.?[0-9]*</validator> + </xml> + + <xml name="input_validator_list_integer"> + <validator type="regex" message="The format is '1,2,4,6'" >[0-9, ]+</validator> + </xml> + + + <token name="@INPUT_IMAGE_LABEL@">RData file</token> + <token name="@INPUT_IMAGE_HELP@">It contains a xcms3::XCMSnExp object (named xdata)</token> + + + <!-- MISC --> + <token name="@HELP_AUTHORS_WRAPPERS@"> + +.. class:: infomark + +**Galaxy integration** ABiMS TEAM - SU/CNRS - Station biologique de Roscoff and Yann Guitton - LABERCA +Part of Workflow4Metabolomics.org [W4M] + + | Contact support@workflow4metabolomics.org for any questions or concerns about the Galaxy implementation of this tool. + + </token> + + <xml name="citation_w4m"> + <citation type="doi">10.1093/bioinformatics/btu813</citation> + </xml> +</macros>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros_xcms.xml Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,272 @@ +<?xml version="1.0"?> +<macros> + + <token name="@TOOL_VERSION@">3.12.0</token> + <xml name="requirements"> + <requirements> + <requirement type="package" version="@TOOL_VERSION@">bioconductor-xcms</requirement> + <requirement type="package" version="1.1_5">r-batch</requirement> + <requirement type="package" version="1.1_2">r-rcolorbrewer</requirement> + <requirement type="package" version="6.0">unzip</requirement> + <yield /> + </requirements> + </xml> + + <!-- FILE_LOAD for planemo test --> + <token name="@COMMAND_FILE_LOAD@"> + #if $file_load_section.file_load_conditional.file_load_select == "yes": + #if $file_load_section.file_load_conditional.input[0].is_of_type("mzxml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzdata") or $file_load_section.file_load_conditional.input[0].is_of_type("netcdf"): + #set singlefile_galaxyPath = '|'.join( [ str( $single_file ) for $single_file in $file_load_section.file_load_conditional.input ] ) + #set singlefile_sampleName = '|'.join( [ str( $single_file.name ) for $single_file in $file_load_section.file_load_conditional.input ] ) + + singlefile_galaxyPath '$singlefile_galaxyPath' singlefile_sampleName '$singlefile_sampleName' + #else + zipfile '$file_load_section.file_load_conditional.input' + #end if + #end if + </token> + + <xml name="input_file_load"> + <section name="file_load_section" title="Resubmit your raw dataset or your zip file"> + <conditional name="file_load_conditional"> + <param name="file_load_select" type="select" label="Resubmit your dataset or your zip file" help="Use only if you get a message saying that your original dataset or zip file have been deleted on the server." > + <option value="no" >no need</option> + <option value="yes" >yes</option> + </param> + <when value="no"> + </when> + <when value="yes"> + <param name="input" type="data" format="mzxml,mzml,mzdata,netcdf,no_unzip.zip,zip" multiple="true" label="File(s) from your history containing your chromatograms" help="Single file mode for the format: mzxml, mzml, mzdata and netcdf. Zip file mode for the format: no_unzip.zip, zip. See the help section below." /> + </when> + </conditional> + </section> + </xml> + + <xml name="test_file_load_zip"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="faahKO_reduce.zip" ftype="zip" /> + </conditional> + </section> + </xml> + + <xml name="test_file_load_zip_sacuri"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="sacuri_dir_root.zip" ftype="zip" /> + </conditional> + </section> + </xml> + + <xml name="test_file_load_single"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="wt15.CDF,ko16.CDF,ko15.CDF,wt16.CDF" ftype="netcdf" /> + </conditional> + </section> + </xml> + + <xml name="test_file_load_single_file" token_filename=""> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="@FILENAME@.CDF" ftype="netcdf" /> + </conditional> + </section> + </xml> + + <!-- PEAKLIST --> + <token name="@COMMAND_PEAKLIST@"> + #if $peaklist.peaklistBool + convertRTMinute $peaklist.convertRTMinute + numDigitsMZ $peaklist.numDigitsMZ + numDigitsRT $peaklist.numDigitsRT + intval $peaklist.intval + naTOzero $peaklist.naTOzero + #end if + </token> + + <xml name="input_peaklist"> + <param name="convertRTMinute" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="Convert retention time (seconds) into minutes" help="Convert the columns rtmed, rtmin and rtmax into minutes"/> + <param name="numDigitsMZ" type="integer" value="4" label="Number of decimal places for mass values reported in ions' identifiers." help="A minimum of 4 decimal places is recommended. Useful to avoid duplicates within identifiers" /> + <param name="numDigitsRT" type="integer" value="0" label="Number of decimal places for retention time values reported in ions' identifiers." help="Useful to avoid duplicates within identifiers" /> + <param name="intval" type="select" label="Reported intensity values" help="[intval] See the help section below"> + <option value="into" selected="true">into</option> + <option value="maxo">maxo</option> + <option value="intb">intb</option> + </param> + <param name="naTOzero" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="If NA values remain, replace them by 0 in the dataMatrix" help="Mandatory for some of the downstream tools (data processing, statistics) that do not accept NA values"/> + </xml> + + <xml name="input_peaklist_section"> + <section name="peaklist" title="Peak List" expanded="True"> + <param name="peaklistBool" type="hidden" label="Get the Peak List" value="True" /> + <expand macro="input_peaklist"/> + </section> + </xml> + + <xml name="input_peaklist_conditional"> + <conditional name="peaklist"> + <param name="peaklistBool" type="boolean" label="Get the Peak List" /> + <when value="true"> + <expand macro="input_peaklist"/> + </when> + <when value="false" /> + </conditional> + </xml> + + <xml name="output_peaklist" token_function=""> + <data name="variableMetadata" format="tabular" label="${image.name[:-6]}.@FUNCTION@.variableMetadata.tsv" from_work_dir="variableMetadata.tsv" > + <filter>(peaklist['peaklistBool'])</filter> + </data> + <data name="dataMatrix" format="tabular" label="${image.name[:-6]}.@FUNCTION@.dataMatrix.tsv" from_work_dir="dataMatrix.tsv" > + <filter>(peaklist['peaklistBool'])</filter> + </data> + </xml> + + <token name="@HELP_PEAKLIST@"> + +Get a Peak List +--------------- + +If 'true', the module generates two additional files corresponding to the peak list: +- the variable metadata file (corresponding to information about extracted ions such as mass or retention time) +- the data matrix (corresponding to related intensities) + +**decimal places for [mass or retention time] values in identifiers** + + | Ions' identifiers are constructed as MxxxTyyy where 'xxx' is the ion median mass and 'yyy' the ion median retention time. + | Two parameters are used to adjust the number of decimal places wanted in identifiers for mass and retention time respectively. + | Theses parameters do not affect decimal places in columns other than the identifier one. + +**Reported intensity values** + + | This parameter determines which values should be reported as intensities in the dataMatrix table; it correspond to xcms 'intval' parameter: + | - into: integrated area of original (raw) peak + | - maxo: maximum intensity of original (raw) peak + | - intb: baseline corrected integrated peak area (only available if peak detection was done by ‘findPeaks.centWave’) + + </token> + + <token name="@HELP_PEAKLIST_OUTPUT@"> +xset.variableMetadata.tsv : tabular format + + | Table containing information about ions; can be used as one input of **Quality_Metrics** or **Generic_filter** modules. + +xset.dataMatrix.tsv : tabular format + + | Table containing ions' intensities; can be used as one input of **Quality_Metrics** or **Generic_filter** modules. + </token> + + <!-- CENTWAVE --> + <token name="@COMMAND_CENTWAVE@"> + ppm $methods.ppm + peakwidth "c($methods.peakwidth)" + + ## Advanced + snthresh $methods.CentWaveAdv.snthresh + prefilter "c($methods.CentWaveAdv.prefilter)" + mzCenterFun $methods.CentWaveAdv.mzCenterFun + integrate $methods.CentWaveAdv.integrate + mzdiff $methods.CentWaveAdv.mzdiff + fitgauss $methods.CentWaveAdv.fitgauss + noise $methods.CentWaveAdv.noise + verboseColumns $methods.CentWaveAdv.verboseColumns + </token> + + <xml name="input_centwave"> + <param argument="ppm" type="float" value="25" label="Max tolerated ppm m/z deviation in consecutive scans in ppm" help="for the initial ROI definition." /> + <param argument="peakwidth" type="text" value="20,50" label="Min,Max peak width in seconds" help="with the expected approximate peak width in chromatographic space."> + <expand macro="input_validator_range_float"/> + </param> + </xml> + + <xml name="input_centwaveAdv"> + <param argument="snthresh" type="integer" value="10" label="Signal to Noise ratio cutoff" /> + <param argument="prefilter" type="text" value="3,100" label="Prefilter step for for the first analysis step (ROI detection)" help="Separate by coma k, I. Mass traces are only retained if they contain at least ‘k‘ peaks with intensity ‘>= I‘."> + <expand macro="input_validator_range_integer"/> + </param> + <param argument="mzCenterFun" type="select" label="Name of the function to calculate the m/z center of the chromatographic peak" > + <option value="wMean">intensity weighted mean of the peak's m/z values</option> + <option value="mean">mean of the peak's m/z values</option> + <option value="apex">use the m/z value at the peak apex</option> + <option value="wMeanApex3">intensity weighted mean of the m/z value at the peak apex and the m/z values left and right of it</option> + <option value="meanApex3">mean of the m/z value of the peak apex and the m/z values left and right of it</option> + </param> + <param argument="integrate" type="select" label="Integration method" > + <option value="1">peak limits are found through descent on the mexican hat filtered data (more robust, but less exact)</option> + <option value="2">peak limits based on real data (more accurate but prone to noise)</option> + </param> + <param argument="mzdiff" type="float" value="-0.001" label="Minimum difference in m/z for peaks with overlapping retention times" help="can be negative to allow overlap" /> + <param argument="fitgauss" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="fitgauss" help="whether or not a Gaussian should be fitted to each peak" /> + <param argument="noise" type="integer" value="0" label="Noise filter" help="allowing to set a minimum intensity required for centroids to be considered in the first analysis step (centroids with intensity lower than ‘noise’ are omitted from ROI detection)." /> + <param argument="verboseColumns" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="verbose Columns" help="whether additional peak meta data columns should be returned" /> + </xml> + + <token name="@COMMAND_CENTWAVEADVROI@"> + #if $sectionROI.roiList: + roiList '$sectionROI.roiList' + firstBaselineCheck $sectionROI.firstBaselineCheck + #if $sectionROI.roiScales != "": + roiScales "c($sectionROI.roiScales)" + #end if + #end if + </token> + + <xml name="input_centwaveAdvROI" token_optional="true"> + <param argument="roiList" type="data" format="tabular" optional="@OPTIONAL@" label="List of regions-of-interest (ROI) representing detected mass traces" help="If ROIs are submitted the first analysis step is omitted and chromatographic peak detection is performed on the submitted ROIs. Each ROI is expected to have the following elements specified: ‘scmin’ (start scan index), ‘scmax’ (end scan index), ‘mzmin’ (minimum m/z), ‘mzmax’ (maximum m/z), ‘length’ (number of scans), ‘intensity’ (summed intensity)." /> + <param argument="firstBaselineCheck" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Is continuous data within regions of interest is checked to be above the first baseline." /> + <param argument="roiScales" type="text" value="" optional="true" label="Numeric vector defining the scale for each region of interest in ‘roiList’" help="Length equal to ‘roiList’ - Should be used for the centWave-wavelets (format 0.9,1,0.2)"> + <expand macro="input_validator_range_float"/> + </param> + </xml> + + <!-- MISC --> + <token name="@HELP_AUTHORS@"> +.. class:: infomark + +**Authors** Colin A. Smith csmith@scripps.edu, Ralf Tautenhahn rtautenh@gmail.com, Steffen Neumann sneumann@ipb-halle.de, Paul Benton hpaul.benton08@imperial.ac.uk and Christopher Conley cjconley@ucdavis.edu + +@HELP_AUTHORS_WRAPPERS@ + +--------------------------------------------------- + + </token> + + <token name="@HELP_XCMS_MANUAL@"> + +For details and explanations concerning all the parameters and workflow of xcms_ package, see its manual_ and this example_ + +.. _xcms: https://bioconductor.org/packages/release/bioc/html/xcms.html +.. _manual: http://www.bioconductor.org/packages/release/bioc/manuals/xcms/man/xcms.pdf +.. _example: https://bioconductor.org/packages/release/bioc/vignettes/xcms/inst/doc/xcms.html +.. _News: https://bioconductor.org/packages/release/bioc/news/xcms/NEWS + + </token> + + <token name="@HELP_XCMS_NEWVERSION_3440@"> +**Version 3.4.4.0 - 08/02/2019** + +- UPGRADE: upgrade the xcms version from 3.0.0 to 3.4.4 (see XCMS News_) + </token> + <token name="@HELP_XCMS_NEWVERSION_3610@"> +**Version 3.6.1+galaxy* - 03/09/2019** + +- UPGRADE: upgrade the xcms version from 3.4.4 to 3.6.1 (see XCMS News_) + </token> + <token name="@HELP_XCMS_NEWVERSION_31200@"> +**Version 3.12.0+galaxy* - 03/03/2020** + +- UPGRADE: upgrade the xcms version from 3.6.1 to 3.12.0 (see XCMS News_) + </token> + + <xml name="citation"> + <citations> + <citation type="doi">10.1021/ac051437y</citation> + <expand macro="citation_w4m"/> + </citations> + </xml> +</macros>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/repository_dependencies.xml Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,5 @@ +<?xml version="1.0" ?> +<repositories> + <repository name="no_unzip_datatype" owner="lecorguille" toolshed="https://toolshed.g2.bx.psu.edu" changeset_revision="7800ba9a4c1e"/> + <repository name="rdata_xcms_datatypes" owner="lecorguille" toolshed="https://toolshed.g2.bx.psu.edu" changeset_revision="544f6d2329ac"/> +</repositories> \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/xcms_refine.r Thu Jun 17 07:42:55 2021 +0000 @@ -0,0 +1,126 @@ +#!/usr/bin/env Rscript + +# ----- LOG FILE ----- +log_file <- file("log.txt", open = "wt") +sink(log_file) +sink(log_file, type = "output") + + +# ----- 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("xcms", "batch", "RColorBrewer") +loadAndDisplayPackages(pkgs) +cat("\n\n"); + +# ----- ARGUMENTS ----- +cat("\tARGUMENTS INFO\n") +# interpretation of arguments given in command line as an R list of objects +args <- parseCommandArgs(evaluate = FALSE) +write.table(as.matrix(args), col.names = F, quote = F, sep = "\t") + +cat("\n\n") + +# ----- PROCESSING INFILE ----- +cat("\tARGUMENTS PROCESSING INFO\n") + +#saving the specific parameters +args_method <- args$method +args_image <- args$image +args_msLevel <- args$msLevel +param_args <- list() + +if (args_method == "CleanPeaks") { + param_args$maxPeakwidth <- args$maxPeakwidth +} else if (args_method == "FilterIntensity") { + param_args$threshold <- args$threshold + param_args$value <- args$value + param_args$nValues <- args$nValues +} else if (args_method == "MergeNeighboringPeaks") { + param_args$expandRt <- args$expandRt + param_args$expandMz <- args$expandMz + param_args$ppm <- args$ppm + param_args$minProp <- args$minProp +} + +cat("\n\n") + + +# ----- ARGUMENTS PROCESSING ----- +cat("\tINFILE PROCESSING INFO\n") + +#image is an .RData file necessary to use xset variable given by previous tools +load(args_image) +if (!exists("xdata")) stop("\n\nERROR: The RData doesn't contain any object called 'xdata'. Such RData as this might have been created by an old version of XMCS 2.*") + +# Handle infiles +if (!exists("singlefile")) singlefile <- NULL +if (!exists("zipfile")) zipfile <- NULL +rawFilePath <- retrieveRawfileInTheWorkingDir(singlefile, zipfile, args) +zipfile <- rawFilePath$zipfile +singlefile <- rawFilePath$singlefile + +cat("\n\n") + + +# ----- MAIN PROCESSING INFO ----- +cat("\tMAIN PROCESSING INFO\n") + + +cat("\t\tPREPARE PARAMETERS\n\n") + +if (args_method == "CleanPeaks") { + refineChromPeaksParam <- CleanPeaksParam(maxPeakwidth = param_args$maxPeakwidth) +} else if (args_method == "FilterIntensity") { + refineChromPeaksParam <- FilterIntensityParam( + threshold = param_args$threshold, + nValues = param_args$nValues, + value = param_args$value + ) +} else if (args_method == "MergeNeighboringPeaks") { + refineChromPeaksParam <- MergeNeighboringPeaksParam( + expandRt = param_args$expandRt, + expandMz = param_args$expandMz, + ppm = param_args$ppm, + minProp = param_args$minProp + ) +} + +cat(str(refineChromPeaksParam)) + +cat("\n\n\t\tCOMPUTE\n") + +xdata <- updateObject(xdata) + +xdata <- refineChromPeaks(xdata, param = refineChromPeaksParam) + +cat("\n\n") + +# ----- EXPORT ----- + +cat("\tXCMSnExp OBJECT INFO\n") +print(xdata) +cat("\n\n") + +cat("\txcmsSet OBJECT INFO\n") +# Get the legacy xcmsSet object +xset <- getxcmsSetObject(xdata) +print(xset) +cat("\n\n") + +#saving R data in .Rdata file to save the variables used in the present tool +objects2save <- c("xdata", "zipfile", "singlefile", "md5sumList", "sampleNamesList") +save(list = objects2save[objects2save %in% ls()], file = "xcmsSet.RData") + +cat("\n\n") + + +cat("\tDONE\n")