Mercurial > repos > workflow4metabolomics > ms2snoop

comparison MS2snoop.R @ 4:856001213966 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 53543b5d911fc1f2d204f314a4d2aaf93a8c7715
author workflow4metabolomics
date Wed, 06 Jul 2022 10:38:39 +0000
parents c68c94865667
children 78d5a12406c2
comparison
equal deleted inserted replaced
3:c68c94865667 4:856001213966
13 #' 13 #'
14 #' @import optparse 14 #' @import optparse
15 #' 15 #'
16 16
17 17
18 assign("MS2SNOOP_VERSION", "1.1.0") 18 assign("MS2SNOOP_VERSION", "2.0.0")
19 lockBinding("MS2SNOOP_VERSION", globalenv()) 19 lockBinding("MS2SNOOP_VERSION", globalenv())
20 20
21 assign("MISSING_PARAMETER_ERROR", 1) 21 assign("MISSING_PARAMETER_ERROR", 1)
22 lockBinding("MISSING_PARAMETER_ERROR", globalenv()) 22 lockBinding("MISSING_PARAMETER_ERROR", globalenv())
23 23
85 fid, 85 fid,
86 sum_int, 86 sum_int,
87 vmz, 87 vmz,
88 cor_abs_int, 88 cor_abs_int,
89 refcol, 89 refcol,
90 c_name 90 c_name,
91 inchikey,
92 elemcomposition
91 ) { 93 ) {
92 ## du fait de la difference de nombre de colonne entre la dataframe qui 94 ## du fait de la difference de nombre de colonne entre la dataframe qui
93 ## inclue les scans en 1ere col, mzRef se decale de 1 95 ## inclue les scans en 1ere col, mzRef se decale de 1
94 refcol <- refcol - 1 96 refcol <- refcol - 1
95 ## compute relative intensities max=100% 97 ## compute relative intensities max=100%
154 156
155 ## prepare result file 157 ## prepare result file
156 corValid <- (round(cor_abs_int, 2) >= r_threshold) ##nolint object_name_linter 158 corValid <- (round(cor_abs_int, 2) >= r_threshold) ##nolint object_name_linter
157 cp_res <- data.frame( 159 cp_res <- data.frame(
158 rep(c_name, length(vmz)), 160 rep(c_name, length(vmz)),
161 rep(inchikey, length(vmz)),
162 rep(elemcomposition, length(vmz)),
159 rep(fid, length(vmz)), 163 rep(fid, length(vmz)),
160 vmz, 164 vmz,
161 cor_abs_int, 165 cor_abs_int,
162 sum_int[-1], 166 sum_int[-1],
163 rel_int, 167 rel_int,
164 corValid 168 corValid
165 ) 169 )
166 170
167 colnames(cp_res) <- c( 171 colnames(cp_res) <- c(
168 "compoundName", 172 "compoundName",
173 "inchikey",
174 "elemcomposition",
169 "fileid", 175 "fileid",
170 "fragments_mz", 176 "fragments_mz",
171 "CorWithPrecursor", 177 "CorWithPrecursor",
172 "AbsoluteIntensity", 178 "AbsoluteIntensity",
173 "relativeIntensity", 179 "relativeIntensity",
198 precursors, 204 precursors,
199 fragments, 205 fragments,
200 mzref, 206 mzref,
201 rtref, 207 rtref,
202 c_name, 208 c_name,
209 inchikey,
210 elemcomposition,
203 min_number_scan, 211 min_number_scan,
204 mzdecimal, 212 mzdecimal,
205 r_threshold = DEFAULT_EXTRACT_FRAGMENTS_R_THRESHOLD, 213 r_threshold = DEFAULT_EXTRACT_FRAGMENTS_R_THRESHOLD,
206 seuil_ra = DEFAULT_EXTRACT_FRAGMENTS_SEUIL_RA, 214 seuil_ra = DEFAULT_EXTRACT_FRAGMENTS_SEUIL_RA,
207 tolmz = DEFAULT_EXTRACT_FRAGMENTS_TOLMZ, 215 tolmz = DEFAULT_EXTRACT_FRAGMENTS_TOLMZ,
359 fid = mf[f], 367 fid = mf[f],
360 sum_int = sum_int, 368 sum_int = sum_int,
361 vmz = vmz, 369 vmz = vmz,
362 cor_abs_int = cor_abs_int, 370 cor_abs_int = cor_abs_int,
363 refcol = refcol, 371 refcol = refcol,
364 c_name = c_name 372 c_name = c_name,
373 inchikey = inchikey,
374 elemcomposition = elemcomposition
365 ) 375 )
366 if (f == 1) { 376 if (f == 1) {
367 res_comp <- res_comp_by_file 377 res_comp <- res_comp_by_file
368 } 378 }
369 } else { 379 } else {
539 ) 549 )
540 return(parser) 550 return(parser)
541 } 551 }
542 552
543 stop_with_status <- function(msg, status) { 553 stop_with_status <- function(msg, status) {
554 sink(stderr())
544 message(sprintf("Error: %s", msg)) 555 message(sprintf("Error: %s", msg))
545 message(sprintf("Error code: %s", status)) 556 message(sprintf("Error code: %s", status))
557 sink(NULL)
546 base::quit(status = status) 558 base::quit(status = status)
547 } 559 }
548 560
549 check_args_validity <- function(args) { ## nolint cyclocomp_linter 561 check_args_validity <- function(args) { ## nolint cyclocomp_linter
550 sysvars <- Sys.getenv() 562 sysvars <- Sys.getenv()
658 } 670 }
659 # } 671 # }
660 return(result) 672 return(result)
661 } 673 }
662 674
675 uniformize_columns <- function(df) {
676 cols <- colnames(df)
677 for (func in c(tolower)) {
678 cols <- func(cols)
679 }
680 colnames(df) <- cols
681 return(df)
682 }
683
663 main <- function(args) { 684 main <- function(args) {
664 if (args$version) { 685 if (args$version) {
665 cat(sprintf("%s\n", MS2SNOOP_VERSION)) 686 cat(sprintf("%s\n", MS2SNOOP_VERSION))
666 base::quit(status = 0) 687 base::quit(status = 0)
667 } 688 }
677 precursors <- get_csv_or_tsv(args$precursors) 698 precursors <- get_csv_or_tsv(args$precursors)
678 ## MSpurity fragments file 699 ## MSpurity fragments file
679 fragments <- get_csv_or_tsv(args$fragments) 700 fragments <- get_csv_or_tsv(args$fragments)
680 ## list of compounds : col1=Name of molecule, col2=m/z, col3=retention time 701 ## list of compounds : col1=Name of molecule, col2=m/z, col3=retention time
681 compounds <- get_csv_or_tsv(args$compounds) 702 compounds <- get_csv_or_tsv(args$compounds)
703
704 compounds <- uniformize_columns(compounds)
705 mandatory_columns <- c(
706 "compound_name",
707 "mz",
708 "rtsec",
709 "inchikey"
710 )
711 presents <- mandatory_columns %in% colnames(compounds)
712 if (!all(presents)) {
713 stop_with_status(
714 sprintf(
715 "Some columns are missing: %s",
716 paste(mandatory_columns[which(!presents)], collapse = ", ")
717 ),
718 BAD_PARAMETER_VALUE_ERROR
719 )
720 }
682 721
683 res_all <- NULL 722 res_all <- NULL
684 for (i in seq_len(nrow(compounds))) { 723 for (i in seq_len(nrow(compounds))) {
685 ## loop execution for all compounds in the compounds file 724 ## loop execution for all compounds in the compounds file
686 res_cor <- NULL 725 res_cor <- NULL
687 res_cor <- extract_fragments( 726 res_cor <- extract_fragments(
688 precursors = precursors, 727 precursors = precursors,
689 fragments = fragments, 728 fragments = fragments,
690 mzref = compounds[[2]][i], 729 mzref = compounds[["mz"]][i],
691 rtref = compounds[[3]][i], 730 rtref = compounds[["rtsec"]][i],
692 c_name = compounds[[1]][i], 731 c_name = compounds[["compound_name"]][i],
732 inchikey = compounds[["inchikey"]][i],
733 elemcomposition = compounds[["elemcomposition"]][i],
693 min_number_scan = args$min_number_scan, 734 min_number_scan = args$min_number_scan,
694 mzdecimal = args$mzdecimal, 735 mzdecimal = args$mzdecimal,
695 r_threshold = args$r_threshold, 736 r_threshold = args$r_threshold,
696 seuil_ra = args$seuil_ra, 737 seuil_ra = args$seuil_ra,
697 tolmz = args$tolmz, 738 tolmz = args$tolmz,
715 sep = "\t", 756 sep = "\t",
716 row.names = FALSE 757 row.names = FALSE
717 ) 758 )
718 } 759 }
719 760
720 unset_debug() 761 global_debug <- FALSE
721 unset_verbose() 762 global_verbose <- FALSE
722 args <- optparse::parse_args(create_parser()) 763 args <- optparse::parse_args(create_parser())
723 main(args) 764 main(args)
724 765
725 warnings() 766 warnings()