Mercurial > repos > tomnl > msp_split

comparison cameradims.R @ 0:aa55ff1d76d4 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/dma-tools-galaxy commit dcfc95273101a7ef0405c2efb8d83f5d456ccd15
author tomnl
date Thu, 10 May 2018 07:20:04 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:aa55ff1d76d4
1 library(optparse)
2 library(cameraDIMS)
3
4 option_list <- list(
5 make_option(c("-i", "--in_file"), type="character"),
6 make_option(c("-o", "--out_dir"), type="character"),
7 make_option("--ppm_iso", default=5),
8 make_option("--ppm_adduct", default=5),
9 make_option("--maxiso", default=4),
10 make_option("--mzabs_iso", default=0.015),
11 make_option("--mzabs_adduct", default=0.015),
12 make_option("--maxcharge", default=3),
13 make_option("--maxmol", default=3),
14 make_option("--polarity", default='pos'),
15 make_option("--rule_path"),
16 make_option("--rule_type"),
17 make_option("--export_ruleset", action="store_true"),
18 make_option("--intensity_filter", default=0)
19 )
20
21
22 # store options
23 opt<- parse_args(OptionParser(option_list=option_list))
24
25 print(sessionInfo())
26 print(opt)
27
28 df <- read.table(opt$in_file, header = TRUE, sep='\t', stringsAsFactors=FALSE)
29
30 print('IN DATA')
31 print(nrow(df))
32 print(head(df))
33
34 devppm_adduct <- opt$ppm_adduct / 1000000
35 devppm_iso <- opt$ppm_iso / 1000000
36
37 paramiso <- list("ppm"=opt$ppm_iso,
38 "filter"=TRUE,
39 "maxcharge"=opt$maxcharge,
40 "maxiso"=opt$maxiso,
41 "mzabs"=opt$mzabs_iso,
42 "intval"='maxo',
43 "minfrac"=0.5,
44 'IM'=NULL,
45 'devppm'=devppm_iso)
46
47 paramadduct <- list("maxCharge"= opt$maxcharge,
48 "maxMol"= opt$maxmol,
49 'devppm'=devppm_adduct,
50 "mzabs"=opt$mzabs_adduct,
51 'IM'=NULL,
52 "filter"=TRUE,
53 'ppm'=opt$ppm_adduct,
54 "quasimolion"= c(1, 6, 8),
55 'polarity'=opt$polarity)
56
57
58
59
60 if(is.null(opt$export_ruleset)){
61 rule_export <- FALSE
62 }else{
63 rule_export <- TRUE
64 }
65
66 df$mz <- as.numeric(df$mz)
67
68 print(head(df$mz))
69
70 if ('intensity' %in% colnames(df)){
71 colnames(df)[colnames(df)=='intensity'] = 'i'
72 }
73
74 if (!'peakID' %in% colnames(df)){
75 df <- cbind('peakID'=1:nrow(df), df)
76 }
77
78 df$i <- as.numeric(df$i)
79
80 df <- df[df$i>opt$intensity_filter,]
81
82 if (nrow(df)==0){
83 print('No peaks left after filtering')
84 quit()
85 }
86
87
88 if (opt$rule_type=='extended_large'){
89
90 rule_sep=','
91 opt$rule_type <- 'user'
92 if(opt$polarity=='pos'){
93 opt$rule_path <- system.file("rules", "CAMERA_rules_PosFinal_PlusLi.csv", package = "cameraDIMS")
94 }else{
95 opt$rule_path <- system.file("rules", "CAMERA_rules_NegFinal_PlusLi.csv", package = "cameraDIMS")
96 }
97 }else{
98 rule_sep='\t'
99 }
100
101 cameraOut <- cameraDIMS(data=df,
102 params_iso=paramiso,
103 params_adduct=paramadduct,
104 rule_type=opt$rule_type,
105 rule_pth=opt$rule_path,
106 rule_sep=rule_sep,
107 rule_export=rule_export)
108
109
110 print(head(cameraOut[[1]]))
111 print(head(cameraOut[[2]]))
112
113
114 out_file1 <- file.path(opt$out_dir, 'camera_annotated_peaklist.txt')
115 out_file2 <- file.path(opt$out_dir, 'camera_annotated_map.txt')
116 out_file3 <- file.path(opt$out_dir, 'ruleset.txt')
117
118 write.table(cameraOut[[1]], out_file1, row.names=FALSE, sep='\t')
119 write.table(cameraOut[[2]], out_file2, row.names=FALSE, sep='\t')
120
121 if (rule_export){
122 write.table(cameraOut[[3]], out_file3, row.names=FALSE, sep='\t')
123 }