Mercurial > repos > lecorguille > camera_annotate

comparison lib.r @ 18:cb923396e70f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 459ef7f63e313493aca32441bd821f09e36de48c
author lecorguille
date Thu, 29 Aug 2019 11:38:21 -0400
parents 73d82de36369
children 01459b73daf9
comparison
equal deleted inserted replaced
17:73d82de36369 18:cb923396e70f
1 # lib.r 1 # lib.r
2
3 #@author G. Le Corguille
4 # solve an issue with batch if arguments are logical TRUE/FALSE
5 parseCommandArgs <- function(...) {
6 args <- batch::parseCommandArgs(...)
7 for (key in names(args)) {
8 if (args[key] %in% c("TRUE","FALSE"))
9 args[key] = as.logical(args[key])
10 }
11 return(args)
12 }
13
14 #@author G. Le Corguille
15 # This function will
16 # - load the packages
17 # - display the sessionInfo
18 loadAndDisplayPackages <- function(pkgs) {
19 for(pkg in pkgs) suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
20
21 sessioninfo = sessionInfo()
22 cat(sessioninfo$R.version$version.string,"\n")
23 cat("Main packages:\n")
24 for (pkg in names(sessioninfo$otherPkgs)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
25 cat("Other loaded packages:\n")
26 for (pkg in names(sessioninfo$loadedOnly)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
27 }
2 28
3 # This function retrieve a xset like object 29 # This function retrieve a xset like object
4 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr 30 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr
5 getxcmsSetObject <- function(xobject) { 31 getxcmsSetObject <- function(xobject) {
6 # XCMS 1.x 32 # XCMS 1.x
32 system(paste0("gm convert ",filebase,"_box/*.png ",pdfBoxOutput)) 58 system(paste0("gm convert ",filebase,"_box/*.png ",pdfBoxOutput))
33 59
34 } 60 }
35 61
36 #@author G. Le Corguille 62 #@author G. Le Corguille
63 #The function create a zip archive from the different png generated by diffreport
64 diffreport_png2zip <- function() {
65 zip("eic.zip", dir(pattern="_eic"), zip=Sys.which("zip"))
66 zip("box.zip", dir(pattern="_box"), zip=Sys.which("zip"))
67 }
68
69 #The function create a zip archive from the different tabular generated by diffreport
70 diffreport_tabular2zip <- function() {
71 zip("tabular.zip", dir(pattern="tabular/*"), zip=Sys.which("zip"))
72 }
73
74 #@author G. Le Corguille
37 #This function convert if it is required the Retention Time in minutes 75 #This function convert if it is required the Retention Time in minutes
38 RTSecondToMinute <- function(variableMetadata, convertRTMinute) { 76 RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
39 if (convertRTMinute){ 77 if (convertRTMinute){
40 #converting the retention times (seconds) into minutes 78 #converting the retention times (seconds) into minutes
41 print("converting the retention times into minutes in the variableMetadata") 79 print("converting the retention times into minutes in the variableMetadata")
55 variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))]) 93 variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))])
56 return(variableMetadata) 94 return(variableMetadata)
57 } 95 }
58 96
59 #The function annotateDiffreport without the corr function which bugs 97 #The function annotateDiffreport without the corr function which bugs
60 annotatediff <- function(xset=xset, listArguments=listArguments, variableMetadataOutput="variableMetadata.tsv") { 98 annotatediff <- function(xset=xset, args=args, variableMetadataOutput="variableMetadata.tsv") {
61 # Resolve the bug with x11, with the function png 99 # Resolve the bug with x11, with the function png
62 options(bitmapType='cairo') 100 options(bitmapType='cairo')
63 101
64 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped. 102 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
65 res=try(is.null(xset@filled)) 103 res=try(is.null(xset@filled))
66 104
67 # ------ annot ------- 105 # ------ annot -------
68 listArguments[["calcCiS"]]=as.logical(listArguments[["calcCiS"]]) 106 args$calcCiS=as.logical(args$calcCiS)
69 listArguments[["calcIso"]]=as.logical(listArguments[["calcIso"]]) 107 args$calcIso=as.logical(args$calcIso)
70 listArguments[["calcCaS"]]=as.logical(listArguments[["calcCaS"]]) 108 args$calcCaS=as.logical(args$calcCaS)
71 109
72 # common parameters 110 # common parameters
73 listArguments4annotate = list(object=xset, 111 args4annotate = list(object=xset,
74 nSlaves=listArguments[["nSlaves"]],sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]], 112 nSlaves=args$nSlaves,sigma=args$sigma,perfwhm=args$perfwhm,
75 maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]], 113 maxcharge=args$maxcharge,maxiso=args$maxiso,minfrac=args$minfrac,
76 ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]], 114 ppm=args$ppm,mzabs=args$mzabs,quick=args$quick,
77 polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]]) 115 polarity=args$polarity,max_peaks=args$max_peaks,intval=args$intval)
78 116
79 # quick == FALSE 117 # quick == FALSE
80 if(listArguments[["quick"]]==FALSE) { 118 if(args$quick==FALSE) {
81 listArguments4annotate = append(listArguments4annotate, 119 args4annotate = append(args4annotate,
82 list(graphMethod=listArguments[["graphMethod"]],cor_eic_th=listArguments[["cor_eic_th"]],pval=listArguments[["pval"]], 120 list(graphMethod=args$graphMethod,cor_eic_th=args$cor_eic_th,pval=args$pval,
83 calcCiS=listArguments[["calcCiS"]],calcIso=listArguments[["calcIso"]],calcCaS=listArguments[["calcCaS"]])) 121 calcCiS=args$calcCiS,calcIso=args$calcIso,calcCaS=args$calcCaS))
84 # no ruleset 122 # no ruleset
85 if (!is.null(listArguments[["multiplier"]])) { 123 if (!is.null(args$multiplier)) {
86 listArguments4annotate = append(listArguments4annotate, 124 args4annotate = append(args4annotate,
87 list(multiplier=listArguments[["multiplier"]])) 125 list(multiplier=args$multiplier))
88 } 126 }
89 # ruleset 127 # ruleset
90 else { 128 else {
91 rulset=read.table(listArguments[["rules"]], h=T, sep=";") 129 rulset=read.table(args$rules, h=T, sep=";")
92 if (ncol(rulset) < 4) rulset=read.table(listArguments[["rules"]], h=T, sep="\t") 130 if (ncol(rulset) < 4) rulset=read.table(args$rules, h=T, sep="\t")
93 if (ncol(rulset) < 4) rulset=read.table(listArguments[["rules"]], h=T, sep=",") 131 if (ncol(rulset) < 4) rulset=read.table(args$rules, h=T, sep=",")
94 if (ncol(rulset) < 4) { 132 if (ncol(rulset) < 4) {
95 error_message="Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation" 133 error_message="Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation"
96 print(error_message) 134 print(error_message)
97 stop(error_message) 135 stop(error_message)
98 } 136 }
99 137
100 listArguments4annotate = append(listArguments4annotate, 138 args4annotate = append(args4annotate,
101 list(rules=rulset)) 139 list(rules=rulset))
102 } 140 }
103 } 141 }
104 142
105 143
106 # launch annotate 144 # launch annotate
107 xa = do.call("annotate", listArguments4annotate) 145 xa = do.call("annotate", args4annotate)
108 peakList=getPeaklist(xa,intval=listArguments[["intval"]]) 146 peakList=getPeaklist(xa,intval=args$intval)
109 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name"); 147 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
110 148
111 # --- Multi condition : diffreport --- 149 # --- Multi condition : diffreport ---
112 diffrepOri=NULL 150 diffrepOri=NULL
113 if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) { 151 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset))>=2) {
114 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped. 152 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
115 res=try(is.null(xset@filled)) 153 res=try(is.null(xset@filled))
116 classes=levels(sampclass(xset)) 154 classes=levels(sampclass(xset))
117 x=1:(length(classes)-1) 155 x=1:(length(classes)-1)
118 for (i in seq(along=x) ) { 156 for (i in seq(along=x) ) {
119 y=1:(length(classes)) 157 y=1:(length(classes))
120 for (n in seq(along=y)){ 158 for (n in seq(along=y)){
121 if(i+n <= length(classes)){ 159 if(i+n <= length(classes)){
122 filebase=paste(classes[i],class2=classes[i+n],sep="-vs-") 160 filebase=paste(classes[i],class2=classes[i+n],sep="-vs-")
123 161
124 diffrep=diffreport(object=xset,class1=classes[i],class2=classes[i+n],filebase=filebase,eicmax=listArguments[["eicmax"]],eicwidth=listArguments[["eicwidth"]],sortpval=TRUE,value=listArguments[["value"]],h=listArguments[["h"]],w=listArguments[["w"]],mzdec=listArguments[["mzdec"]],missing=0) 162 diffrep=diffreport(
163 object=xset,class1=classes[i],class2=classes[i+n],
164 filebase=filebase,eicmax=args$eicmax,eicwidth=args$eicwidth,
165 sortpval=TRUE,value=args$value,h=args$h,w=args$w,mzdec=args$mzdec,missing=0)
125 166
126 diffrepOri = diffrep 167 diffrepOri = diffrep
127 168
128 # renamming of the column rtmed to rt to fit with camera peaklist function output 169 # renamming of the column rtmed to rt to fit with camera peaklist function output
129 colnames(diffrep)[colnames(diffrep)=="rtmed"] <- "rt" 170 colnames(diffrep)[colnames(diffrep)=="rtmed"] <- "rt"
131 172
132 # combines results and reorder columns 173 # combines results and reorder columns
133 diffrep = merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F) 174 diffrep = merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F)
134 diffrep = cbind(diffrep[,!(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))],diffrep[,(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))]) 175 diffrep = cbind(diffrep[,!(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))],diffrep[,(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))])
135 176
136 diffrep = RTSecondToMinute(diffrep, listArguments[["convertRTMinute"]]) 177 diffrep = RTSecondToMinute(diffrep, args$convertRTMinute)
137 diffrep = formatIonIdentifiers(diffrep, numDigitsRT=listArguments[["numDigitsRT"]], numDigitsMZ=listArguments[["numDigitsMZ"]]) 178 diffrep = formatIonIdentifiers(diffrep, numDigitsRT=args$numDigitsRT, numDigitsMZ=args$numDigitsMZ)
138 179
139 if(listArguments[["sortpval"]]){ 180 if(args$sortpval){
140 diffrep=diffrep[order(diffrep$pvalue), ] 181 diffrep=diffrep[order(diffrep$pvalue), ]
141 } 182 }
142 183
143 dir.create("tabular") 184 dir.create("tabular", showWarnings = FALSE)
144 write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste("tabular/",filebase,"_tsv.tabular",sep="")) 185 write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste("tabular/",filebase,"_tsv.tabular",sep=""))
145 186
146 if (listArguments[["eicmax"]] != 0) { 187 if (args$eicmax != 0) {
147 diffreport_png2pdf(filebase) 188 if (args$png2 == "pdf")
189 diffreport_png2pdf(filebase)
148 } 190 }
149 } 191 }
150 } 192 }
151 } 193 }
194 if (args$png2 == "zip")
195 diffreport_png2zip()
196 if (args$tabular2 == "zip")
197 diffreport_tabular2zip()
152 } 198 }
153 199
154 # --- variableMetadata --- 200 # --- variableMetadata ---
155 variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))] 201 variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
156 variableMetadata = RTSecondToMinute(variableMetadata, listArguments[["convertRTMinute"]]) 202 variableMetadata = RTSecondToMinute(variableMetadata, args$convertRTMinute)
157 variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=listArguments[["numDigitsRT"]], numDigitsMZ=listArguments[["numDigitsMZ"]]) 203 variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=args$numDigitsRT, numDigitsMZ=args$numDigitsMZ)
158 # if we have 2 conditions, we keep stat of diffrep 204 # if we have 2 conditions, we keep stat of diffrep
159 if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))==2) { 205 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset))==2) {
160 variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F) 206 variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F)
161 if(exists("listArguments[[\"sortpval\"]]")){ 207 if(exists("args[[\"sortpval\"]]")){
162 variableMetadata=variableMetadata[order(variableMetadata$pvalue), ] 208 variableMetadata=variableMetadata[order(variableMetadata$pvalue), ]
163 } 209 }
164 } 210 }
165 211
166 variableMetadataOri=variableMetadata 212 variableMetadataOri=variableMetadata
169 return(list("xa"=xa,"diffrep"=diffrepOri,"variableMetadata"=variableMetadataOri)); 215 return(list("xa"=xa,"diffrep"=diffrepOri,"variableMetadata"=variableMetadataOri));
170 216
171 } 217 }
172 218
173 219
174 combinexsAnnos_function <- function(xaP, xaN, listOFlistArgumentsP,listOFlistArgumentsN, diffrepP=NULL,diffrepN=NULL,pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, convertRTMinute=F, numDigitsMZ=0, numDigitsRT=0, variableMetadataOutput="variableMetadata.tsv"){ 220 combinexsAnnos_function <- function(xaP, xaN, diffrepP=NULL,diffrepN=NULL,
221 pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, convertRTMinute=F, numDigitsMZ=0,
222 numDigitsRT=0, variableMetadataOutput="variableMetadata.tsv"){
175 223
176 #Load the two Rdata to extract the xset objects from positive and negative mode 224 #Load the two Rdata to extract the xset objects from positive and negative mode
177 cat("\tObject xset from positive mode\n") 225 cat("\tObject xset from positive mode\n")
178 print(xaP) 226 print(xaP)
179 cat("\n") 227 cat("\n")
196 stop("You must relauch the CAMERA.annotate step with the lastest version.") 244 stop("You must relauch the CAMERA.annotate step with the lastest version.")
197 } 245 }
198 246
199 if(pos){ 247 if(pos){
200 xa=xaP 248 xa=xaP
201 listOFlistArgumentsP=listOFlistArguments
202 mode="neg. Mode" 249 mode="neg. Mode"
203 } else { 250 } else {
204 xa=xaN 251 xa=xaN
205 listOFlistArgumentsN=listOFlistArguments
206 mode="pos. Mode" 252 mode="pos. Mode"
207 } 253 }
208 254
209 peakList=getPeaklist(xa) 255 peakList=getPeaklist(xa)
210 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name"); 256 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
237 return(variableMetadata); 283 return(variableMetadata);
238 284
239 } 285 }
240 286
241 # This function get the raw file path from the arguments 287 # This function get the raw file path from the arguments
242 getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) { 288 getRawfilePathFromArguments <- function(singlefile, zipfile, args) {
243 if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]] 289 if (!is.null(args$zipfile)) zipfile = args$zipfile
244 if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]] 290 if (!is.null(args$zipfilePositive)) zipfile = args$zipfilePositive
245 if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]] 291 if (!is.null(args$zipfileNegative)) zipfile = args$zipfileNegative
246 292
247 if (!is.null(listArguments[["singlefile_galaxyPath"]])) { 293 if (!is.null(args$singlefile_galaxyPath)) {
248 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]]; 294 singlefile_galaxyPaths = args$singlefile_galaxyPath;
249 singlefile_sampleNames = listArguments[["singlefile_sampleName"]] 295 singlefile_sampleNames = args$singlefile_sampleName
250 } 296 }
251 if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) { 297 if (!is.null(args$singlefile_galaxyPathPositive)) {
252 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]]; 298 singlefile_galaxyPaths = args$singlefile_galaxyPathPositive;
253 singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]] 299 singlefile_sampleNames = args$singlefile_sampleNamePositive
254 } 300 }
255 if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) { 301 if (!is.null(args$singlefile_galaxyPathNegative)) {
256 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]]; 302 singlefile_galaxyPaths = args$singlefile_galaxyPathNegative;
257 singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]] 303 singlefile_sampleNames = args$singlefile_sampleNameNegative
258 } 304 }
259 if (exists("singlefile_galaxyPaths")){ 305 if (exists("singlefile_galaxyPaths")){
260 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,",")) 306 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,","))
261 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,",")) 307 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,","))
262 308
265 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i] 311 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i]
266 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i] 312 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i]
267 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath 313 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath
268 } 314 }
269 } 315 }
270 for (argument in c("zipfile","zipfilePositive","zipfileNegative","singlefile_galaxyPath","singlefile_sampleName","singlefile_galaxyPathPositive","singlefile_sampleNamePositive","singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) { 316 for (argument in c("zipfile", "zipfilePositive", "zipfileNegative",
271 listArguments[[argument]]=NULL 317 "singlefile_galaxyPath", "singlefile_sampleName",
272 } 318 "singlefile_galaxyPathPositive", "singlefile_sampleNamePositive",
273 return(list(zipfile=zipfile, singlefile=singlefile, listArguments=listArguments)) 319 "singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) {
320 args[[argument]]=NULL
321 }
322 return(list(zipfile=zipfile, singlefile=singlefile, args=args))
274 } 323 }
275 324
276 325
277 # This function retrieve the raw file in the working directory 326 # This function retrieve the raw file in the working directory
278 # - if zipfile: unzip the file with its directory tree 327 # - if zipfile: unzip the file with its directory tree