Mercurial > repos > recetox > isolib

comparison isolib.R @ 5:964b4559eb1b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/RECETOX/galaxytools/tree/master/tools/isolib commit 0765a69d2180d0cfda663cc1b4585a8935142169
author recetox
date Thu, 15 May 2025 16:42:56 +0000
parents 2b1118bce0b1
children f0fe957df1cc
comparison
equal deleted inserted replaced
4:2b1118bce0b1 5:964b4559eb1b
4 library(MetaboCoreUtils) 4 library(MetaboCoreUtils)
5 library(readr) 5 library(readr)
6 library(tidyselect) 6 library(tidyselect)
7 7
8 8
9 isotopes <- data.frame(
10 element = character(),
11 abundance = numeric(),
12 isotope = character()
13 )
14 element <- character()
15 abundance <- numeric()
16 adducts <- data.frame(
17 Name = character(),
18 Mult = numeric(),
19 Formula_add = character(),
20 Formula_ded = character(),
21 Charge = numeric(),
22 Ion_mode = character(),
23 Mass = numeric()
24 )
25 mz <- numeric()
26 peaks <- numeric()
27 isos <- character()
28 retention_time <- numeric()
29 full_formula <- character()
30 name <- character()
31
32
9 parse_args <- function() { 33 parse_args <- function() {
10 args <- commandArgs(trailingOnly = TRUE) 34 args <- commandArgs(trailingOnly = TRUE)
11 35
12 compound_table <- read_tsv( 36 compound_table <- read_tsv(
13 file = args[1], 37 file = args[1],
14 col_types = "ccd", 38 col_types = "ccd",
15 col_select = all_of(c("name", "formula")) | any_of("rt") 39 col_select = all_of(c("name", "formula")) | any_of("rt")
16 ) 40 )
41 # Handle missing or empty rel_to argument
42 rel_to_value <- if (length(args) >= 8 && args[8] != "") {
43 if (args[8] == "none") 0 else as.numeric(args[8])
44 } else {
45 0 # Default value is 0
46 }
47
48 if (!rel_to_value %in% c(0, 1, 2, 3, 4)) {
49 stop(
50 "Invalid value for rel_to. Expected 'none' (0),",
51 " or a numeric value between 0 and 4."
52 )
53 }
17 54
18 parsed <- list( 55 parsed <- list(
19 compound_table = compound_table, 56 compound_table = compound_table,
20 adducts_to_use = c(unlist(strsplit(args[2], ",", fixed = TRUE))), 57 adducts_to_use = c(unlist(strsplit(args[2], ",", fixed = TRUE))),
21 threshold = as.numeric(args[3]), 58 threshold = as.numeric(args[3]),
22 append_adducts = args[4], 59 append_adducts = args[4],
23 append_isotopes = args[5], 60 append_isotopes = args[5],
24 out_format = args[6], 61 out_format = args[6],
25 outfile = args[7] 62 outfile = args[7],
26 ) 63 rel_to = rel_to_value
27 return(parsed) 64 )
65 parsed
28 } 66 }
29 67
30 generate_isotope_spectra <- function(compound_table, 68 generate_isotope_spectra <- function(compound_table,
31 adducts_to_use, 69 adducts_to_use,
32 append_adducts, 70 append_adducts,
33 threshold) { 71 threshold,
72 rel_to) {
34 data(isotopes) 73 data(isotopes)
35 data(adducts) 74 data(adducts)
36 75
37 monoisotopic <- isotopes |> 76 monoisotopic <- isotopes |>
38 dplyr::group_by(element) |> 77 dplyr::group_by(element) |>
39 dplyr::slice_max(abundance, n = 1) |> 78 dplyr::slice_max(abundance, n = 1) |>
40 dplyr::filter(!stringr::str_detect(element, "\\[|\\]")) 79 dplyr::filter(!stringr::str_detect(element, "\\[|\\]"))
41 80
42 chemforms <- check_chemform(isotopes, compound_table$formula)[, 2] 81 chemforms <- enviPat::check_chemform(isotopes, compound_table$formula)[, 2]
43 spectra <- data.frame() 82 spectra <- data.frame()
44 83
45 for (current in adducts_to_use) { 84 for (current in adducts_to_use) {
46 adduct <- adducts[adducts$Name == current, ] 85 adduct <- adducts[adducts$Name == current, ]
47 multiplied_chemforms <- multiform(chemforms, adduct$Mult) 86 multiplied_chemforms <- enviPat::multiform(chemforms, adduct$Mult)
48 87
49 if (adduct$Ion_mode == "negative") { 88 if (adduct$Ion_mode == "negative") {
50 merged_chemforms <- subform(multiplied_chemforms, adduct$Formula_ded) 89 merged_chemforms <- enviPat::subform(
90 multiplied_chemforms,
91 adduct$Formula_ded
92 )
51 } else { 93 } else {
52 merged_chemforms <- mergeform(multiplied_chemforms, adduct$Formula_add) 94 merged_chemforms <- enviPat::mergeform(
95 multiplied_chemforms,
96 adduct$Formula_add
97 )
53 } 98 }
54 99
55 charge_string <- paste0( 100 charge_string <- paste0(
56 if (adduct$Charge > 0) "+" else "-", 101 if (adduct$Charge > 0) "+" else "-",
57 if (abs(adduct$Charge) > 1) abs(adduct$Charge) else "" 102 if (abs(adduct$Charge) > 1) abs(adduct$Charge) else ""
58 ) 103 )
59 adduct_string <- paste0("[", adduct$Name, "]", charge_string) 104 adduct_string <- paste0("[", adduct$Name, "]", charge_string)
60 precursor_mz <- calculateMass(multiplied_chemforms) + adduct$Mass 105 precursor_mass <- MetaboCoreUtils::calculateMass(multiplied_chemforms)
106 precursor_mz <- precursor_mass + adduct$Mass
61 107
62 if (append_adducts == TRUE) { 108 if (append_adducts == TRUE) {
63 names <- paste( 109 names <- paste(
64 compound_table$name, 110 compound_table$name,
65 paste0("(", adduct$Name, ")"), 111 paste0("(", adduct$Name, ")"),
86 patterns <- enviPat::isopattern( 132 patterns <- enviPat::isopattern(
87 isotopes = isotopes, 133 isotopes = isotopes,
88 chemforms = merged_chemforms, 134 chemforms = merged_chemforms,
89 charge = adduct$Charge, 135 charge = adduct$Charge,
90 threshold = threshold, 136 threshold = threshold,
137 rel_to = rel_to,
91 ) 138 )
92 139
93 mzs <- list() 140 mzs <- list()
94 intensities <- list() 141 intensities <- list()
95 isos <- list() 142 isos <- list()
124 spectra_df$mz <- mzs 171 spectra_df$mz <- mzs
125 spectra_df$intensity <- intensities 172 spectra_df$intensity <- intensities
126 spectra_df$isotopes <- isos 173 spectra_df$isotopes <- isos
127 spectra <- rbind(spectra, spectra_df) 174 spectra <- rbind(spectra, spectra_df)
128 } 175 }
129 return(spectra) 176 spectra
130 } 177 }
131 178
132 write_to_msp <- function(spectra, file) { 179 write_to_msp <- function(spectra, file) {
133 sps <- Spectra(dplyr::select(spectra, -isotopes)) 180 sps <- Spectra::Spectra(dplyr::select(spectra, -isotopes))
134 export(sps, MsBackendMsp(), file = file) 181 Spectra::export(sps, MsBackendMsp::MsBackendMsp(), file = file)
135 } 182 }
136 183
137 write_to_table <- function(spectra, file, append_isotopes) { 184 write_to_table <- function(spectra, file, append_isotopes) {
138 entries <- spectra |> 185 entries <- spectra |>
139 dplyr::rowwise() |> 186 dplyr::rowwise() |>
164 args <- parse_args() 211 args <- parse_args()
165 spectra <- generate_isotope_spectra( 212 spectra <- generate_isotope_spectra(
166 args$compound_table, 213 args$compound_table,
167 args$adducts_to_use, 214 args$adducts_to_use,
168 args$append_adducts, 215 args$append_adducts,
169 args$threshold 216 args$threshold,
217 args$rel_to
170 ) 218 )
171 219
172 if (args$out_format == "msp") { 220 if (args$out_format == "msp") {
173 write_to_msp(spectra, args$outfile) 221 write_to_msp(spectra, args$outfile)
174 } else if (args$out_format == "tabular") { 222 } else if (args$out_format == "tabular") {