Mercurial > repos > computational-metabolomics > mspurity_filterfragspectra
diff flagRemove.R @ 6:cc0f8ddad4a8 draft
"planemo upload for repository https://github.com/computational-metabolomics/mspurity-galaxy commit 2579c8746819670348c378f86116f83703c493eb"
author | computational-metabolomics |
---|---|
date | Thu, 04 Mar 2021 12:28:18 +0000 |
parents | 1e6b7c4d99cd |
children | a46824d13914 |
line wrap: on
line diff
--- a/flagRemove.R Fri Nov 13 10:16:03 2020 +0000 +++ b/flagRemove.R Thu Mar 04 12:28:18 2021 +0000 @@ -2,105 +2,100 @@ library(optparse) print(sessionInfo()) option_list <- list( - make_option(c("-o", "--out_dir"), type="character", default=getwd(), - help="Output folder for resulting files [default = %default]" + make_option(c("-o", "--out_dir"), type = "character", default = getwd(), + help = "Output folder for resulting files [default = %default]" ), - make_option(c("-x", "--xset_path"), type="character", default=file.path(getwd(),"xset.rds"), - help="The path to the xcmsSet object [default = %default]" + make_option(c("-x", "--xset_path"), type = "character", default = file.path(getwd(), "xset.rds"), + help = "The path to the xcmsSet object [default = %default]" ), - make_option("--polarity", default=NA, - help="polarity (just used for naming purpose for files being saved) [positive, negative, NA] [default %default]" + make_option("--polarity", default = NA, + help = "polarity (just used for naming purpose for files being saved) [positive, negative, NA] [default %default]" ), - make_option("--rsd_i_blank", default=100, - help="RSD threshold for the blank [default = %default]" + make_option("--rsd_i_blank", default = 100, + help = "RSD threshold for the blank [default = %default]" ), - make_option("--minfrac_blank", default=0.5, - help="minimum fraction of files for features needed for the blank [default = %default]" + make_option("--minfrac_blank", default = 0.5, + help = "minimum fraction of files for features needed for the blank [default = %default]" ), - make_option("--rsd_rt_blank", default=100, - help="RSD threshold for the RT of the blank [default = %default]" + make_option("--rsd_rt_blank", default = 100, + help = "RSD threshold for the RT of the blank [default = %default]" ), - make_option("--ithres_blank", default=0, - help="Intensity threshold for the blank [default = %default]" + make_option("--ithres_blank", default = 0, + help = "Intensity threshold for the blank [default = %default]" ), - make_option("--s2b", default=10, - help="fold change (sample/blank) needed for sample peak to be allowed. e.g. + make_option("--s2b", default = 10, + help = "fold change (sample/blank) needed for sample peak to be allowed. e.g. if s2b set to 10 and the recorded sample 'intensity' value was 100 and blank was 10. 1000/10 = 100, so sample has fold change higher than the threshold and the peak is not considered a blank [default = %default]" ), - make_option("--blank_class", default='blank', type="character", - help="A string representing the class that will be used for the blank.[default = %default]" + make_option("--blank_class", default = "blank", type = "character", + help = "A string representing the class that will be used for the blank.[default = %default]" ), - make_option("--egauss_thr", default=NA, - help="Threshold for filtering out non gaussian shaped peaks. Note this only works + make_option("--egauss_thr", default = NA, + help = "Threshold for filtering out non gaussian shaped peaks. Note this only works if the 'verbose columns' and 'fit gauss' was used with xcms [default = %default]" ), - make_option("--rsd_i_sample", default=100, - help="RSD threshold for the samples [default = %default]" + make_option("--rsd_i_sample", default = 100, + help = "RSD threshold for the samples [default = %default]" ), - make_option("--minfrac_sample", default=0.8, - help="minimum fraction of files for features needed for the samples [default = %default]" + make_option("--minfrac_sample", default = 0.8, + help = "minimum fraction of files for features needed for the samples [default = %default]" ), - make_option("--rsd_rt_sample", default=100, - help="RSD threshold for the RT of the samples [default = %default]" + make_option("--rsd_rt_sample", default = 100, + help = "RSD threshold for the RT of the samples [default = %default]" ), - make_option("--ithres_sample", default=5000, - help="Intensity threshold for the sample [default = %default]" + make_option("--ithres_sample", default = 5000, + help = "Intensity threshold for the sample [default = %default]" ), - make_option("--grp_rm_ids", default=NA, - help="vector of grouped_xcms peaks to remove (corresponds to the row from xcms::group output) + make_option("--grp_rm_ids", default = NA, + help = "vector of grouped_xcms peaks to remove (corresponds to the row from xcms::group output) [default = %default]" ), - make_option("--remove_spectra", action="store_true", - help=" TRUE if flagged spectra is to be removed [default = %default]" + make_option("--remove_spectra", action = "store_true", + help = "TRUE if flagged spectra is to be removed [default = %default]" ), - make_option("--minfrac_xcms", default=0.5, - help="minfrac for xcms grouping [default = %default]" + make_option("--minfrac_xcms", default = 0.5, + help = "minfrac for xcms grouping [default = %default]" ), - make_option("--mzwid", default=0.001, - help="mzwid for xcms grouping [default = %default]" + make_option("--mzwid", default = 0.001, + help = "mzwid for xcms grouping [default = %default]" ), - make_option("--bw", default=5, - help="bw for xcms grouping [default = %default]" + make_option("--bw", default = 5, + help = "bw for xcms grouping [default = %default]" ), - make_option("--temp_save", action="store_true", - help="Assign True if files for each step saved (for testing purposes) [default = %default]" + make_option("--temp_save", action = "store_true", + help = "Assign True if files for each step saved (for testing purposes) [default = %default]" ), - make_option("--samplelist", type="character", help="Sample list to determine the blank class") - - - - + make_option("--samplelist", type = "character", help = "Sample list to determine the blank class") ) - #make_option("--multilist", action="store_true" - # help="NOT CURRENTLY IMPLEMENTED: If paired blank removal is to be performed a - multilist - sample list file has to be provided" - #), +# nolint start +# make_option("--multilist", action="store_true" +# help="NOT CURRENTLY IMPLEMENTED: If paired blank removal is to be performed a - multilist - sample list file has to be provided" +# ), +# nolint end # store options -opt<- parse_args(OptionParser(option_list=option_list)) +opt <- parse_args(OptionParser(option_list = option_list)) opt <- replace(opt, opt == "NA", NA) - - - -if (is.null(opt$temp_save)){ - temp_save<-FALSE +if (is.null(opt$temp_save)) { + temp_save <- FALSE }else{ - temp_save<-TRUE + temp_save <- TRUE } -if (is.null(opt$remove_spectra)){ - remove_spectra<-FALSE +if (is.null(opt$remove_spectra)) { + remove_spectra <- FALSE }else{ - remove_spectra<-TRUE + remove_spectra <- TRUE } @@ -109,35 +104,35 @@ getxcmsSetObject <- function(xobject) { # XCMS 1.x if (class(xobject) == "xcmsSet") - return (xobject) + return(xobject) # XCMS 3.x if (class(xobject) == "XCMSnExp") { # Get the legacy xcmsSet object - suppressWarnings(xset <- as(xobject, 'xcmsSet')) - sampclass(xset) <- xset@phenoData$sample_group - return (xset) + suppressWarnings(xset <- as(xobject, "xcmsSet")) + xcms::sampclass(xset) <- xset@phenoData$sample_group + return(xset) } } -loadRData <- function(rdata_path, name){ +loadRData <- function(rdata_path, name) { #loads an RData file, and returns the named xset object if it is there load(rdata_path) return(get(ls()[ls() %in% name])) } -xset <- getxcmsSetObject(loadRData(opt$xset_path, c('xset','xdata'))) +xset <- getxcmsSetObject(loadRData(opt$xset_path, c("xset", "xdata"))) print(xset) -if (is.null(opt$samplelist)){ +if (is.null(opt$samplelist)) { blank_class <- opt$blank_class }else{ - samplelist <- read.table(opt$samplelist, sep='\t', header=TRUE) - samplelist_blank <- unique(samplelist$sample_class[samplelist$blank=='yes']) + samplelist <- read.table(opt$samplelist, sep = "\t", header = TRUE) + samplelist_blank <- unique(samplelist$sample_class[samplelist$blank == "yes"]) chosen_blank <- samplelist_blank[samplelist_blank %in% xset@phenoData$class] - if (length(chosen_blank)>1){ - print('ERROR: only 1 blank is currently allowed to be used with this tool') + if (length(chosen_blank) > 1) { + print("ERROR: only 1 blank is currently allowed to be used with this tool") quit() } blank_class <- as.character(chosen_blank) @@ -145,47 +140,47 @@ } -if (is.null(opt$multilist)){ +if (is.null(opt$multilist)) { ffrm_out <- flag_remove(xset, - pol=opt$polarity, - rsd_i_blank=opt$rsd_i_blank, - minfrac_blank=opt$minfrac_blank, - rsd_rt_blank=opt$rsd_rt_blank, - ithres_blank=opt$ithres_blank, - s2b=opt$s2b, - ref.class=blank_class, - egauss_thr=opt$egauss_thr, - rsd_i_sample=opt$rsd_i_sample, - minfrac_sample=opt$minfrac_sample, - rsd_rt_sample=opt$rsd_rt_sample, - ithres_sample=opt$ithres_sample, - minfrac_xcms=opt$minfrac_xcms, - mzwid=opt$mzwid, - bw=opt$bw, - out_dir=opt$out_dir, - temp_save=temp_save, - remove_spectra=remove_spectra, - grp_rm_ids=unlist(strsplit(as.character(opt$grp_rm_ids), split=", "))[[1]]) - print('flag remove finished') + pol = opt$polarity, + rsd_i_blank = opt$rsd_i_blank, + minfrac_blank = opt$minfrac_blank, + rsd_rt_blank = opt$rsd_rt_blank, + ithres_blank = opt$ithres_blank, + s2b = opt$s2b, + ref.class = blank_class, + egauss_thr = opt$egauss_thr, + rsd_i_sample = opt$rsd_i_sample, + minfrac_sample = opt$minfrac_sample, + rsd_rt_sample = opt$rsd_rt_sample, + ithres_sample = opt$ithres_sample, + minfrac_xcms = opt$minfrac_xcms, + mzwid = opt$mzwid, + bw = opt$bw, + out_dir = opt$out_dir, + temp_save = temp_save, + remove_spectra = remove_spectra, + grp_rm_ids = unlist(strsplit(as.character(opt$grp_rm_ids), split = ", "))[[1]]) + print("flag remove finished") xset <- ffrm_out[[1]] grp_peaklist <- ffrm_out[[2]] removed_peaks <- ffrm_out[[3]] - save.image(file=file.path(opt$out_dir, 'xset_filtered.RData'), version=2) + save.image(file = file.path(opt$out_dir, "xset_filtered.RData"), version = 2) # grpid needed for mspurity ID needed for deconrank... (will clean up at some up) - peak_pth <- file.path(opt$out_dir, 'peaklist_filtered.tsv') + peak_pth <- file.path(opt$out_dir, "peaklist_filtered.tsv") print(peak_pth) - write.table(data.frame('grpid'=rownames(grp_peaklist), 'ID'=rownames(grp_peaklist), grp_peaklist), - peak_pth, row.names=FALSE, sep='\t') + write.table(data.frame("grpid" = rownames(grp_peaklist), "ID" = rownames(grp_peaklist), grp_peaklist), + peak_pth, row.names = FALSE, sep = "\t") removed_peaks <- data.frame(removed_peaks) - write.table(data.frame('ID'=rownames(removed_peaks),removed_peaks), - file.path(opt$out_dir, 'removed_peaks.tsv'), row.names=FALSE, sep='\t') + write.table(data.frame("ID" = rownames(removed_peaks), removed_peaks), + file.path(opt$out_dir, "removed_peaks.tsv"), row.names = FALSE, sep = "\t") }else{ - + # nolint start # TODO #xsets <- split(xset, multilist_df$multlist) # @@ -193,12 +188,10 @@ # #for (mgrp in mult_grps){ # xset_i <- xsets[mgrp] - # xcms::group(xset_i, + # xcms::group(xset_i, # # } - + # nolint end } - -