Mercurial > repos > workflow4metabolomics > camera_findisotopes

comparison CAMERA_findIsotopes.R @ 0:a8aa53d54c88 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:30:26 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a8aa53d54c88
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 print("Argument types:")
29 print(sapply(args, class))
30
31 # Check if the image file exists
32 if (!file.exists(args$image)) {
33 stop("The RData file does not exist: ", args$image)
34 }
35
36 # ----- PROCESSING INFILE -----
37
38 # Load the RData file
39 load(args$image)
40 args$image <- NULL
41
42 # Save arguments to generate a report
43 if (!exists("listOFlistArguments")) listOFlistArguments <- list()
44 listOFlistArguments[[format(Sys.time(), "%y%m%d-%H:%M:%S_findIsotopes")]] <- args
45
46 # We unzip automatically the chromatograms from the zip files.
47 if (!exists("zipfile")) zipfile <- NULL
48 if (!exists("singlefile")) singlefile <- NULL
49 rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args)
50 zipfile <- rawFilePath$zipfile
51 singlefile <- rawFilePath$singlefile
52 args <- rawFilePath$args
53
54 print(paste("singlefile :", singlefile))
55 if (!is.null(singlefile)) {
56 directory <- retrieveRawfileInTheWorkingDir(singlefile, zipfile)
57 }
58
59 # Verify if the xa object is loaded
60 if (!exists("xa")) {
61 stop("The xa object was not found in the RData file.")
62 }
63
64 print("xa object loaded:")
65 print(xa)
66
67 print("calcIsotopeMatrix")
68 calcIsotopeMatrix <- function(maxiso = 4) {
69 if (!is.numeric(maxiso)) {
70 stop("Parameter maxiso is not numeric!\n")
71 } else if (maxiso < 1 || maxiso > 8) {
72 stop(paste(
73 "Parameter maxiso must be between 1 and 8. ",
74 "Otherwise, use your own IsotopeMatrix.\n"
75 ), sep = "")
76 }
77
78 isotopeMatrix <- matrix(NA, 8, 4)
79 colnames(isotopeMatrix) <- c("mzmin", "mzmax", "intmin", "intmax")
80
81 isotopeMatrix[1, ] <- c(1.000, 1.0040, 1.0, 150)
82 isotopeMatrix[2, ] <- c(0.997, 1.0040, 0.01, 200)
83 isotopeMatrix[3, ] <- c(1.000, 1.0040, 0.001, 200)
84 isotopeMatrix[4, ] <- c(1.000, 1.0040, 0.0001, 200)
85 isotopeMatrix[5, ] <- c(1.000, 1.0040, 0.00001, 200)
86 isotopeMatrix[6, ] <- c(1.000, 1.0040, 0.000001, 200)
87 isotopeMatrix[7, ] <- c(1.000, 1.0040, 0.0000001, 200)
88 isotopeMatrix[8, ] <- c(1.000, 1.0040, 0.00000001, 200)
89
90 return(isotopeMatrix[1:maxiso, , drop = FALSE])
91 }
92
93 isotopeMatrix <- calcIsotopeMatrix(args$maxiso)
94
95 # Apply the findIsotopes function on the xsAnnotate object
96 print("Calling findIsotopes function:")
97 xa <- findIsotopes(xa, maxcharge = args$maxcharge, maxiso = args$maxiso, ppm = args$ppm, mzabs = args$mzabs, intval = args$intval, minfrac = args$minfrac, isotopeMatrix = isotopeMatrix, filter = args$filter)
98
99 print("Result of the findIsotopes function:")
100 print(xa)
101
102 # Extract the list of annotated peaks
103 peakList <- getPeaklist(xa, intval = args$intval)
104
105 if (length(phenoData@data$sample_name) == 1) {
106 peakList$name <- make.unique(paste0("M", round(peakList[, "mz"], 0), "T", round(peakList[, "rt"], 0)), "_")
107 variableMetadata <- peakList[, c("name", setdiff(names(peakList), "name"))]
108 variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT = args$numDigitsRT, numDigitsMZ = args$numDigitsMZ)
109 } else {
110 names_default <- groupnames(xa@xcmsSet, mzdec = 0, rtdec = 0) # Names without decimals
111 names_custom <- groupnames(xa@xcmsSet, mzdec = args$numDigitsMZ, rtdec = args$numDigitsRT) # Names with "x" decimals
112
113 variableMetadata <- data.frame(
114 name = names_default,
115 name_custom = names_custom,
116 stringsAsFactors = FALSE
117 )
118 variableMetadata <- cbind(variableMetadata, peakList[, !(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))])
119 }
120
121 if (!exists("RTinMinute")) RTinMinute <- FALSE
122
123 if (args$convertRTMinute && RTinMinute == FALSE) {
124 RTinMinute <- TRUE
125 variableMetadata <- RTSecondToMinute(variableMetadata = variableMetadata, convertRTMinute = args$convertRTMinute)
126 }
127
128 # Saves the extracted peak list as a TSV file named 'variableMetadata.tsv'
129 output_file_tsv <- "variableMetadata.tsv"
130 write.table(variableMetadata, file = output_file_tsv, sep = "\t", row.names = FALSE, quote = FALSE)
131
132 # Save the updated xsAnnotate object
133 output_file_RData <- "camera_findIsotopes.RData"
134
135 objects2save <- c("xa", "variableMetadata", "listOFlistArguments", "zipfile", "singlefile", "RTinMinute", "phenoData")
136 save(list = objects2save[objects2save %in% ls()], file = output_file_RData)
137
138 cat("Output files generated:", output_file_tsv, "and", output_file_RData, "\n")