Mercurial > repos > computational-metabolomics > mspurity_frag4feature

comparison averageFragSpectra.R @ 0:85bc606fd219 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/computational-metabolomics/mspurity-galaxy commit 2579c8746819670348c378f86116f83703c493eb"
author computational-metabolomics
date Thu, 04 Mar 2021 12:22:55 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:85bc606fd219
1 library(optparse)
2 library(msPurity)
3 library(xcms)
4 print(sessionInfo())
5
6
7 get_av_spectra <- function(x) {
8
9 if (length(x$av_intra) > 0) {
10 av_intra_df <- plyr::ldply(x$av_intra)
11
12 if (nrow(av_intra_df) == 0) {
13 av_intra_df <- NULL
14 }else{
15 av_intra_df$method <- "intra"
16 }
17
18 }else{
19 av_intra_df <- NULL
20 }
21
22 if ((is.null(x$av_inter)) || (nrow(x$av_inter) == 0)) {
23 av_inter_df <- NULL
24 }else{
25 av_inter_df <- x$av_inter
26 av_inter_df$method <- "inter"
27 }
28
29 if ((is.null(x$av_all)) || (nrow(x$av_all) == 0)) {
30 av_all_df <- NULL
31 }else{
32 av_all_df <- x$av_all
33 av_all_df$method <- "all"
34 }
35
36 combined <- plyr::rbind.fill(av_intra_df, av_inter_df, av_all_df)
37
38 return(combined)
39 }
40
41
42 option_list <- list(
43 make_option("--out_rdata", type = "character"),
44 make_option("--out_peaklist", type = "character"),
45 make_option("--pa", type = "character"),
46 make_option("--av_level", type = "character"),
47 make_option("--minfrac", default = 0.5),
48 make_option("--minnum", default = 1),
49 make_option("--ppm", default = 5.0),
50 make_option("--snr", default = 0),
51 make_option("--ra", default = 0),
52 make_option("--av", default = "median", type = "character"),
53 make_option("--sumi", action = "store_true"),
54 make_option("--rmp", action = "store_true"),
55 make_option("--cores", default = 1)
56 )
57
58 opt <- parse_args(OptionParser(option_list = option_list))
59 print(opt)
60
61
62 load_r_data <- function(rdata_path, name) {
63 #loads an RData file, and returns the named xset object if it is there
64 load(rdata_path)
65 return(get(ls()[ls() %in% name]))
66 }
67
68 # Requires
69 pa <- load_r_data(opt$pa, "pa")
70
71 pa@cores <- opt$cores
72
73 if (is.null(opt$rmp)) {
74 rmp <- FALSE
75 }else{
76 rmp <- TRUE
77 }
78
79 if (is.null(opt$sumi)) {
80 sumi <- FALSE
81 }else{
82 sumi <- TRUE
83 }
84
85 if (opt$av_level == "intra") {
86 pa <- msPurity::averageIntraFragSpectra(pa,
87 minfrac = opt$minfrac,
88 minnum = opt$minnum,
89 ppm = opt$ppm,
90 snr = opt$snr,
91 ra = opt$ra,
92 av = opt$av,
93 sumi = sumi,
94 rmp = rmp,
95 cores = opt$cores)
96
97 } else if (opt$av_level == "inter") {
98
99 pa <- msPurity::averageInterFragSpectra(pa,
100 minfrac = opt$minfrac,
101 minnum = opt$minnum,
102 ppm = opt$ppm,
103 snr = opt$snr,
104 ra = opt$ra,
105 av = opt$av,
106 sumi = sumi,
107 rmp = rmp,
108 cores = opt$cores)
109 } else if (opt$av_level == "all") {
110
111 pa <- msPurity::averageAllFragSpectra(pa,
112 minfrac = opt$minfrac,
113 minnum = opt$minnum,
114 ppm = opt$ppm,
115 snr = opt$snr,
116 ra = opt$ra,
117 av = opt$av,
118 sumi = sumi,
119 rmp = rmp,
120 cores = opt$cores)
121 }
122
123 print(pa)
124 save(pa, file = opt$out_rdata)
125
126 if (length(pa) > 0) {
127
128 av_spectra <- plyr::ldply(pa@av_spectra, get_av_spectra)
129
130 if (nrow(av_spectra) == 0) {
131 message("No average spectra available")
132 } else {
133 colnames(av_spectra)[1] <- "grpid"
134 av_spectra$grpid <- names(pa@av_spectra)[av_spectra$grpid]
135
136 if ((length(pa@av_intra_params) > 0) || (length(pa@av_inter_params) > 0)) {
137 # Add some extra info (only required if av_intra or av_inter performed)
138 colnames(av_spectra)[2] <- "fileid"
139 av_spectra$avid <- seq_len(nrow(av_spectra))
140
141 filenames <- sapply(av_spectra$fileid,
142 function(x) names(pa@fileList)[as.integer(x)])
143 # filenames_galaxy <- sapply(
144 # av_spectra$fileid, function(x) basename(pa@fileList[as.integer(x)]))
145
146 av_spectra <- as.data.frame(
147 append(av_spectra, list(filename = filenames), after = 2))
148 }
149
150
151 print(head(av_spectra))
152 write.table(av_spectra, opt$out_peaklist, row.names = FALSE, sep = "\t")
153
154 }
155 }