Mercurial > repos > computational-metabolomics > mspurity_createdatabase

comparison averageFragSpectra.R @ 0:f52287a06c02 draft

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 cb903cd93f9378cfb5eeb68512a54178dcea7bbc-dirty"
author computational-metabolomics
date Wed, 27 Nov 2019 13:44:58 -0500
parents
children 2f71b3495221
comparison
equal deleted inserted replaced
-1:000000000000 0:f52287a06c02
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
47 make_option("--av_level", type="character"),
48
49 make_option("--minfrac", default=0.5),
50 make_option("--minnum", default=1),
51 make_option("--ppm", default=5.0),
52
53 make_option("--snr", default=0),
54
55 make_option("--ra", default=0),
56
57 make_option("--av", default="median", type="character"),
58 make_option("--sumi", action="store_true"),
59
60 make_option("--rmp", action="store_true"),
61 make_option("--cores", default=1)
62 )
63
64 opt <- parse_args(OptionParser(option_list=option_list))
65 print(opt)
66
67
68 loadRData <- function(rdata_path, name){
69 #loads an RData file, and returns the named xset object if it is there
70 load(rdata_path)
71 return(get(ls()[ls() %in% name]))
72 }
73
74 # Requires
75 pa <- loadRData(opt$pa, 'pa')
76
77 pa@cores <- opt$cores
78
79 if(is.null(opt$rmp)){
80 rmp = FALSE
81 }else{
82 rmp = TRUE
83 }
84
85 if(is.null(opt$sumi)){
86
87 sumi = FALSE
88 }else{
89 sumi = TRUE
90
91 }
92
93
94 if(opt$av_level=="intra"){
95
96 pa <- msPurity::averageIntraFragSpectra(pa,
97 minfrac=opt$minfrac,
98 minnum=opt$minnum,
99 ppm=opt$ppm,
100 snr=opt$snr,
101 ra=opt$ra,
102 av=opt$av,
103 sumi=sumi,
104 rmp=rmp,
105 cores=opt$cores)
106
107 } else if(opt$av_level=="inter"){
108
109 pa <- msPurity::averageInterFragSpectra(pa,
110 minfrac=opt$minfrac,
111 minnum=opt$minnum,
112 ppm=opt$ppm,
113 snr=opt$snr,
114 ra=opt$ra,
115 av=opt$av,
116 sumi=sumi,
117 rmp=rmp,
118 cores=opt$cores)
119 } else if(opt$av_level=="all"){
120
121 pa <- msPurity::averageAllFragSpectra(pa,
122 minfrac=opt$minfrac,
123 minnum=opt$minnum,
124 ppm=opt$ppm,
125 snr=opt$snr,
126 ra=opt$ra,
127 av=opt$av,
128 sumi=sumi,
129 rmp=rmp,
130 cores=opt$cores)
131
132 }
133
134 print(pa)
135 save(pa, file=opt$out_rdata)
136
137
138 if (length(pa)>0){
139
140 av_spectra <- plyr::ldply(pa@av_spectra, get_av_spectra)
141
142 if (nrow(av_spectra)==0){
143 message('No average spectra available')
144 } else{
145 colnames(av_spectra)[1] <- 'grpid'
146 av_spectra$grpid <- names(pa@av_spectra)[av_spectra$grpid]
147
148 if((length(pa@av_intra_params)>0) || (length(pa@av_inter_params)>0) ){
149 # Add some extra info (only required if av_intra or av_inter performed)
150 colnames(av_spectra)[2] <- 'fileid'
151 av_spectra$avid <- 1:nrow(av_spectra)
152
153 filenames <- sapply(av_spectra$fileid, function(x) names(pa@fileList)[as.integer(x)])
154 # filenames_galaxy <- sapply(av_spectra$fileid, function(x) basename(pa@fileList[as.integer(x)]))
155
156 av_spectra = as.data.frame(append(av_spectra, list(filename = filenames), after=2))
157 }
158
159
160 print(head(av_spectra))
161 write.table(av_spectra, opt$out_peaklist, row.names=FALSE, sep='\t')
162
163 }
164 }
165