Mercurial > repos > yguitton > metams_rungc

comparison metams.r @ 1:142fbe102a9d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded version 2.0
author yguitton
date Wed, 24 May 2017 07:25:50 -0400
parents 2066efbafd7c
children c75532b75ba1
comparison
equal deleted inserted replaced
0:2066efbafd7c 1:142fbe102a9d
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file 1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
2 # metams.r version="0.99.9" 2 # metams.r version="2.0"
3 #created by Yann GUITTON 3 #created by Yann GUITTON
4 4 #use RI options
5 5
6 #Redirect all stdout to the log file 6 #Redirect all stdout to the log file
7 log_file=file("metams.log", open = "wt") 7 log_file=file("metams.log", open = "wt")
8 sink(log_file) 8 sink(log_file)
9 sink(log_file, type = "out") 9 sink(log_file, type = "out")
18 } 18 }
19 print("step1") 19 print("step1")
20 20
21 21
22 listArguments = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects 22 listArguments = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects
23 print("new version 8") 23 print("new version 2.0")
24 24
25 25
26 print(listArguments) 26 print(listArguments)
27 27
28
28 if (listArguments[["ri"]]!="NULL"){ 29 if (listArguments[["ri"]]!="NULL"){
29 RIarg=read.table(listArguments[["ri"]], h=T) 30 RIarg=read.table(listArguments[["ri"]])
31 if (ncol(RIarg) < 2) RIarg=read.table(listArguments[["ri"]], h=T, sep=";")
32 if (ncol(RIarg) < 2) RIarg=read.table(listArguments[["ri"]], h=T, sep="\t")
33 if (ncol(RIarg) < 2) RIarg=read.table(listArguments[["ri"]], h=T, sep=",")
34 if (ncol(RIarg) < 2) {
35 error_message="Your RI file seems not well formatted. The column separators accepted are ; , and tabulation"
36 print(error_message)
37 stop(error_message)
38 }
39 #to do check real column names
30 colnames(RIarg)<-c("rt","RI") 40 colnames(RIarg)<-c("rt","RI")
31 # print(RIarg) 41 # print(RIarg)
32 } else { 42 } else {
33 RIarg = NULL 43 RIarg = NULL
34 # cat("Ri= ",RIarg) 44 # cat("Ri= ",RIarg)
35 } 45 }
36 46
47 if (listArguments[["rishift"]]!="none"){
48 RIshift=listArguments[["rishift"]]
49 cat("Rishift used= ",RIshift, "\n")
50 } else {
51 RIshift = "none"
52 cat("Rishift NONE= ",RIshift, "\n")
53 }
37 54
38 DBarg=listArguments[["db"]] 55 DBarg=listArguments[["db"]]
39 # if (listArguments[["use_db"]]!="NULL"){ 56 # if (listArguments[["use_db"]]!="NULL"){
40 if (DBarg!="NULL"){ 57 if (DBarg!="NULL"){
41 DBarg=listArguments[["db"]] 58 DBarg=listArguments[["db"]]
42 cat("Db= ",DBarg) 59 cat("Db= ",DBarg, "\n")
43 } else { 60 } else {
44 DBarg = NULL 61 DBarg = NULL
45 cat("NO Db : ",DBarg) 62 cat("NO Db : ",DBarg, "\n")
46 } 63 }
47
48 64
49 65
50 66
51 #for unknown EIC printing 67 #for unknown EIC printing
52 68
53 if (listArguments[["unkn"]]!="NULL") { 69 if (listArguments[["unkn"]][1]!="NULL") {
70 unknarg<-listArguments[["unkn"]]
71 } else {
54 unknarg<-"" 72 unknarg<-""
55 } else {
56 unknarg<-listArguments[["unkn"]]
57 } 73 }
58 74
59 print(paste("Unkn:",unknarg)) 75 print(paste("Unkn:",unknarg))
60 76
61 #remove unused arguments 77 #remove unused arguments
62 listArguments[["unkn"]]<-NULL 78 listArguments[["unkn"]]<-NULL
63 listArguments[["db"]] <- NULL 79 listArguments[["db"]] <- NULL
64 listArguments[["ri"]] <- NULL 80 listArguments[["ri"]] <- NULL
81 listArguments[["rishift"]] <- NULL
65 82
66 print(" step2") 83 print(" step2")
67 84
68 #runGC accept either a list of files a zip folder or an xset object from xcms.xcmsSet tool 85 #runGC accept either a list of files a zip folder or an xset object from xcms.xcmsSet tool
69 86
162 179
163 if (!is.null(DBarg)){ 180 if (!is.null(DBarg)){
164 manual <- read.msp(DBarg) 181 manual <- read.msp(DBarg)
165 DBarg <- createSTDdbGC(stdInfo = NULL, settings = TSQXLS.GC, manualDB = manual) 182 DBarg <- createSTDdbGC(stdInfo = NULL, settings = TSQXLS.GC, manualDB = manual)
166 } 183 }
167 nSlaves=listArguments[["nSlaves"]] 184
185 #use RI instead of rt for time comparison vs DB
186 if (RIshift!="none"){
187 TSQXLS.GC@match2DB.timeComparison<-"RI"
188 TSQXLS.GC@match2DB.RIdiff<-as.numeric(RIshift)
189 TSQXLS.GC@betweenSamples.timeComparison<-"RI"
190 TSQXLS.GC@betweenSamples.RIdiff<-as.numeric(RIshift)
191 }
192
193 nSlaves=listArguments[["nSlaves"]]
194
195
168 if(!metams_zip_file=="") { 196 if(!metams_zip_file=="") {
169 resGC<-runGC(files=samples,settings=TSQXLS.GC, rtrange=rtrange, DB= DBarg, removeArtefacts = TRUE, findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves) #default settings for GC from Wehrens et al 197 resGC<-runGC(files=samples,settings=TSQXLS.GC, rtrange=rtrange, DB= DBarg, removeArtefacts = TRUE, findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves) #default settings for GC from Wehrens et al
170 } 198 }
171 if(!is.null(xsetparam)){ 199 if(!is.null(xsetparam)){
172 settingslist=TSQXLS.GC 200 settingslist=TSQXLS.GC
187 215
188 if (listArguments[["settings"]]=="User_defined") { 216 if (listArguments[["settings"]]=="User_defined") {
189 listArguments[["settings"]]=NULL #delete from the list of arguments 217 listArguments[["settings"]]=NULL #delete from the list of arguments
190 fwhmparam=listArguments[["fwhm"]] 218 fwhmparam=listArguments[["fwhm"]]
191 rtdiffparam=listArguments[["rtdiff"]] 219 rtdiffparam=listArguments[["rtdiff"]]
192 #RIdiffparam=listArguments[["RIdiff"]] #only if timeComparison = "RI"
193 minfeatparam=listArguments[["minfeat"]] 220 minfeatparam=listArguments[["minfeat"]]
194 simthreshparam=listArguments[["simthreshold"]] 221 simthreshparam=listArguments[["simthreshold"]]
195 minclassfractionparam=listArguments[["minclassfraction"]] 222 minclassfractionparam=listArguments[["minclassfraction"]]
196 minclasssizeparam=listArguments[["minclasssize"]] 223 minclasssizeparam=listArguments[["minclasssize"]]
224
197 if (listArguments[["rtrange"]]!="NULL") { 225 if (listArguments[["rtrange"]]!="NULL") {
198 rtrange=listArguments[["rtrange"]] 226 rtrange=listArguments[["rtrange"]]
199 cat("rtrange= ",rtrange) 227 cat("rtrange= ",rtrange)
200 } else { 228 } else {
201 rtrange=NULL 229 rtrange=NULL
202 cat("rtrange= ",rtrange) 230 cat("rtrange= ",rtrange)
203 } 231 }
232
233
204 nSlaves=listArguments[["nSlaves"]] 234 nSlaves=listArguments[["nSlaves"]]
205 235
206 GALAXY.GC <- metaMSsettings("protocolName" = "GALAXY.GC", 236 GALAXY.GC <- metaMSsettings("protocolName" = "GALAXY.GC",
207 "chrom" = "GC", 237 "chrom" = "GC",
208 PeakPicking = list( 238 PeakPicking = list(
225 simthresh = simthreshparam, 255 simthresh = simthreshparam,
226 timeComparison = "rt", 256 timeComparison = "rt",
227 rtdiff = rtdiffparam, 257 rtdiff = rtdiffparam,
228 RIdiff = 5, 258 RIdiff = 5,
229 minfeat = minfeatparam) 259 minfeat = minfeatparam)
230 #to be used when DB option will be available 260
231 # metaSetting(GALAXY.GC, "matchIrrelevants") <- list( 261 #to used if contaminant filter
232 # irrelevantClasses = c("Bleeding", "Plasticizers"), 262
233 # timeComparison = "rt", 263 # metaSetting(GALAXY.GC, "matchIrrelevants") <- list(
234 # RIdiff = 2, 264 # irrelevantClasses = c("Bleeding", "Plasticizers"),
235 # rtdiff = rtdiffparam, 265 # timeComparison = "RI",
236 # simthresh = simthreshparam) 266 # RIdiff = RIdiffparam,
267 # rtdiff = rtdiffparam,
268 # simthresh = simthreshparam)
269
237 metaSetting(GALAXY.GC, "betweenSamples") <- list( 270 metaSetting(GALAXY.GC, "betweenSamples") <- list(
238 min.class.fraction = minclassfractionparam, 271 min.class.fraction = minclassfractionparam,
239 min.class.size = minclasssizeparam, 272 min.class.size = minclasssizeparam,
240 timeComparison = "rt", 273 timeComparison = "rt",
241 rtdiff = rtdiffparam, 274 rtdiff = rtdiffparam,
242 RIdiff = 2, 275 RIdiff = 2,
243 simthresh = simthreshparam) 276 simthresh = simthreshparam)
244 277
245 278 #ONLY use RI instead of rt for time comparison vs DB or samples
279 if (RIshift!="none"){
280 GALAXY.GC@match2DB.timeComparison<-"RI"
281 GALAXY.GC@match2DB.RIdiff<-as.numeric(RIshift)
282 GALAXY.GC@betweenSamples.timeComparison<-"RI"
283 GALAXY.GC@betweenSamples.RIdiff<-as.numeric(RIshift)
284 }
246 # files, xset, settings, rtrange = NULL, DB = NULL, 285 # files, xset, settings, rtrange = NULL, DB = NULL,
247 # removeArtefacts = TRUE, findUnknowns = nexp > 1, 286 # removeArtefacts = TRUE, findUnknowns = nexp > 1,
248 # returnXset = FALSE, RIstandards = NULL, nSlaves = 0 287 # returnXset = FALSE, RIstandards = NULL, nSlaves = 0
249 if (!is.null(DBarg)){ 288 if (!is.null(DBarg)){
250 manual <- read.msp(DBarg) 289 manual <- read.msp(DBarg)
300 print("BPCs") 339 print("BPCs")
301 c<-getBPC2s(files=samples, rt="raw", pdfname="BPCs_raw.pdf") 340 c<-getBPC2s(files=samples, rt="raw", pdfname="BPCs_raw.pdf")
302 341
303 print("Step QC plot") 342 print("Step QC plot")
304 343
344 #to do check if no peaks found
305 #Quality controls plots but only working in R (don't know why) 345 #Quality controls plots but only working in R (don't know why)
306 a<-plotUnknowns(resGC=resGC, unkn=unknarg) #use unknparam value 346 a<-plotUnknowns(resGC=resGC, unkn=unknarg) #use unknparam value
307 347
308 # create a mergpdf 348 # create a mergpdf
309 349