Mercurial > repos > lecorguille > xcms_merge

comparison lib.r @ 10:47e953d9da82 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/xcms commit 49203f8a5271fa5e6bb889e907df71ebf7757309
author lecorguille
date Thu, 08 Mar 2018 05:52:52 -0500
parents 6b5504f877ff
children 67ab853b89f3
comparison
equal deleted inserted replaced
9:35b9bb3205d8 10:47e953d9da82
72 } 72 }
73 73
74 #@author G. Le Corguille 74 #@author G. Le Corguille
75 # Draw the plotChromPeakDensity 3 per page in a pdf file 75 # Draw the plotChromPeakDensity 3 per page in a pdf file
76 getPlotAdjustedRtime <- function(xdata) { 76 getPlotAdjustedRtime <- function(xdata) {
77
77 pdf(file="raw_vs_adjusted_rt.pdf", width=16, height=12) 78 pdf(file="raw_vs_adjusted_rt.pdf", width=16, height=12)
79
78 # Color by group 80 # Color by group
79 group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))] 81 group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))]
80 names(group_colors) <- unique(xdata$sample_group) 82 names(group_colors) <- unique(xdata$sample_group)
81 plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group]) 83 plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group])
82 legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1) 84 legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1)
85
83 # Color by sample 86 # Color by sample
84 plotAdjustedRtime(xdata, col = rainbow(length(xdata@phenoData@data$sample_name))) 87 plotAdjustedRtime(xdata, col = rainbow(length(xdata@phenoData@data$sample_name)))
85 legend("topright", legend=xdata@phenoData@data$sample_name, col=rainbow(length(xdata@phenoData@data$sample_name)), cex=0.8, lty=1) 88 legend("topright", legend=xdata@phenoData@data$sample_name, col=rainbow(length(xdata@phenoData@data$sample_name)), cex=0.8, lty=1)
89
86 dev.off() 90 dev.off()
87 } 91 }
88 92
89 #@author G. Le Corguille 93 #@author G. Le Corguille
90 # value: intensity values to be used into, maxo or intb 94 # value: intensity values to be used into, maxo or intb
102 write.table(variableMetadata, file=variableMetadataOutput,sep="\t",quote=F,row.names=F) 106 write.table(variableMetadata, file=variableMetadataOutput,sep="\t",quote=F,row.names=F)
103 write.table(dataMatrix, file=dataMatrixOutput,sep="\t",quote=F,row.names=F) 107 write.table(dataMatrix, file=dataMatrixOutput,sep="\t",quote=F,row.names=F)
104 108
105 } 109 }
106 110
107 #@author Y. Guitton 111 #@author G. Le Corguille
108 getBPC <- function(file,rtcor=NULL, ...) { 112 getPlotChromatogram <- function(xdata, pdfname="Chromatogram.pdf", aggregationFun = "max") {
109 object <- xcmsRaw(file) 113
110 sel <- profRange(object, ...) 114 chrom <- chromatogram(xdata, aggregationFun = aggregationFun)
111 cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE])) 115 if (aggregationFun == "sum")
112 #plotChrom(xcmsRaw(file), base=T) 116 type="Total Ion Chromatograms"
113 } 117 else
114 118 type="Base Peak Intensity Chromatograms"
115 #@author Y. Guitton 119
116 getBPCs <- function (xcmsSet=NULL, pdfname="BPCs.pdf",rt=c("raw","corrected"), scanrange=NULL) { 120 adjusted="Raw"
117 cat("Creating BIC pdf...\n") 121 if (hasAdjustedRtime(xdata))
118 122 adjusted="Adjusted"
119 if (is.null(xcmsSet)) { 123
120 cat("Enter an xcmsSet \n") 124 main <- paste(type,":",adjusted,"data")
121 stop() 125
122 } else { 126 pdf(pdfname, width=16, height=10)
123 files <- filepaths(xcmsSet) 127
124 } 128 # Color by group
125 129 group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))]
126 phenoDataClass <- as.vector(levels(xcmsSet@phenoData[,"class"])) #sometime phenoData have more than 1 column use first as class 130 names(group_colors) <- unique(xdata$sample_group)
127 131 plot(chrom, col = group_colors[chrom$sample_group], main=main)
128 classnames <- vector("list",length(phenoDataClass)) 132 legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1)
129 for (i in 1:length(phenoDataClass)){ 133
130 classnames[[i]] <- which( xcmsSet@phenoData[,"class"]==phenoDataClass[i]) 134 # Color by sample
131 } 135 plot(chrom, col = rainbow(length(xdata@phenoData@data$sample_name)), main=main)
132 136 legend("topright", legend=xdata@phenoData@data$sample_name, col=rainbow(length(xdata@phenoData@data$sample_name)), cex=0.8, lty=1)
133 N <- dim(phenoData(xcmsSet))[1] 137
134 138 dev.off()
135 TIC <- vector("list",N) 139 }
136 140
137 141 #@author G. Le Corguille
138 for (j in 1:N) { 142 getPlotTICs <- function(xdata, pdfname="TICs.pdf") {
139 143 getPlotChromatogram(xdata, pdfname, aggregationFun = "sum")
140 TIC[[j]] <- getBPC(files[j]) 144 }
141 #good for raw 145
142 # seems strange for corrected 146 #@author G. Le Corguille
143 #errors if scanrange used in xcmsSetgeneration 147 getPlotBPIs <- function(xdata, pdfname="BPIs.pdf") {
144 if (!is.null(xcmsSet) && rt == "corrected") 148 getPlotChromatogram(xdata, pdfname, aggregationFun = "max")
145 rtcor <- xcmsSet@rt$corrected[[j]] 149 }
146 else
147 rtcor <- NULL
148
149 TIC[[j]] <- getBPC(files[j],rtcor=rtcor)
150 # TIC[[j]][,1]<-rtcor
151 }
152
153
154
155 pdf(pdfname,w=16,h=10)
156 cols <- rainbow(N)
157 lty <- 1:N
158 pch <- 1:N
159 #search for max x and max y in BPCs
160 xlim <- range(sapply(TIC, function(x) range(x[,1])))
161 ylim <- range(sapply(TIC, function(x) range(x[,2])))
162 ylim <- c(-ylim[2], ylim[2])
163
164
165 ##plot start
166
167 if (length(phenoDataClass)>2){
168 for (k in 1:(length(phenoDataClass)-1)){
169 for (l in (k+1):length(phenoDataClass)){
170 #print(paste(phenoDataClass[k],"vs",phenoDataClass[l],sep=" "))
171 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k]," vs ",phenoDataClass[l], sep=""), xlab="Retention Time (min)", ylab="BPC")
172 colvect <- NULL
173 for (j in 1:length(classnames[[k]])) {
174 tic <- TIC[[classnames[[k]][j]]]
175 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
176 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
177 colvect <- append(colvect,cols[classnames[[k]][j]])
178 }
179 for (j in 1:length(classnames[[l]])) {
180 # i <- class2names[j]
181 tic <- TIC[[classnames[[l]][j]]]
182 points(tic[,1]/60, -tic[,2], col=cols[classnames[[l]][j]], pch=pch[classnames[[l]][j]], type="l")
183 colvect <- append(colvect,cols[classnames[[l]][j]])
184 }
185 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col=colvect, lty=lty, pch=pch)
186 }
187 }
188 }#end if length >2
189
190 if (length(phenoDataClass)==2){
191 k <- 1
192 l <- 2
193 colvect <- NULL
194 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k],"vs",phenoDataClass[l], sep=""), xlab="Retention Time (min)", ylab="BPC")
195
196 for (j in 1:length(classnames[[k]])) {
197
198 tic <- TIC[[classnames[[k]][j]]]
199 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
200 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
201 colvect<-append(colvect,cols[classnames[[k]][j]])
202 }
203 for (j in 1:length(classnames[[l]])) {
204 # i <- class2names[j]
205 tic <- TIC[[classnames[[l]][j]]]
206 points(tic[,1]/60, -tic[,2], col=cols[classnames[[l]][j]], pch=pch[classnames[[l]][j]], type="l")
207 colvect <- append(colvect,cols[classnames[[l]][j]])
208 }
209 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col=colvect, lty=lty, pch=pch)
210
211 }#end length ==2
212
213 #case where only one class
214 if (length(phenoDataClass)==1){
215 k <- 1
216 ylim <- range(sapply(TIC, function(x) range(x[,2])))
217 colvect <- NULL
218 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k], sep=""), xlab="Retention Time (min)", ylab="BPC")
219
220 for (j in 1:length(classnames[[k]])) {
221 tic <- TIC[[classnames[[k]][j]]]
222 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
223 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
224 colvect <- append(colvect,cols[classnames[[k]][j]])
225 }
226
227 legend("topright",paste(basename(files[c(classnames[[k]])])), col=colvect, lty=lty, pch=pch)
228
229 }#end length ==1
230
231 dev.off() #pdf(pdfname,w=16,h=10)
232
233 invisible(TIC)
234 }
235
236
237
238 #@author Y. Guitton
239 getTIC <- function(file, rtcor=NULL) {
240 object <- xcmsRaw(file)
241 cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object, mzrange=range(object@env$mz))$intensity)
242 }
243
244 #overlay TIC from all files in current folder or from xcmsSet, create pdf
245 #@author Y. Guitton
246 getTICs <- function(xcmsSet=NULL,files=NULL, pdfname="TICs.pdf", rt=c("raw","corrected")) {
247 cat("Creating TIC pdf...\n")
248
249 if (is.null(xcmsSet)) {
250 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]")
251 filepattern <- paste(paste("\\.", filepattern, "$", sep=""), collapse="|")
252 if (is.null(files))
253 files <- getwd()
254 info <- file.info(files)
255 listed <- list.files(files[info$isdir], pattern=filepattern, recursive=TRUE, full.names=TRUE)
256 files <- c(files[!info$isdir], listed)
257 } else {
258 files <- filepaths(xcmsSet)
259 }
260
261 phenoDataClass <- as.vector(levels(xcmsSet@phenoData[,"class"])) #sometime phenoData have more than 1 column use first as class
262 classnames <- vector("list",length(phenoDataClass))
263 for (i in 1:length(phenoDataClass)){
264 classnames[[i]] <- which( xcmsSet@phenoData[,"class"]==phenoDataClass[i])
265 }
266
267 N <- length(files)
268 TIC <- vector("list",N)
269
270 for (i in 1:N) {
271 if (!is.null(xcmsSet) && rt == "corrected")
272 rtcor <- xcmsSet@rt$corrected[[i]] else
273 rtcor <- NULL
274 TIC[[i]] <- getTIC(files[i], rtcor=rtcor)
275 }
276
277 pdf(pdfname, w=16, h=10)
278 cols <- rainbow(N)
279 lty <- 1:N
280 pch <- 1:N
281 #search for max x and max y in TICs
282 xlim <- range(sapply(TIC, function(x) range(x[,1])))
283 ylim <- range(sapply(TIC, function(x) range(x[,2])))
284 ylim <- c(-ylim[2], ylim[2])
285
286
287 ##plot start
288 if (length(phenoDataClass)>2){
289 for (k in 1:(length(phenoDataClass)-1)){
290 for (l in (k+1):length(phenoDataClass)){
291 #print(paste(phenoDataClass[k],"vs",phenoDataClass[l],sep=" "))
292 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k]," vs ",phenoDataClass[l], sep=""), xlab="Retention Time (min)", ylab="TIC")
293 colvect <- NULL
294 for (j in 1:length(classnames[[k]])) {
295 tic <- TIC[[classnames[[k]][j]]]
296 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
297 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
298 colvect <- append(colvect,cols[classnames[[k]][j]])
299 }
300 for (j in 1:length(classnames[[l]])) {
301 # i=class2names[j]
302 tic <- TIC[[classnames[[l]][j]]]
303 points(tic[,1]/60, -tic[,2], col=cols[classnames[[l]][j]], pch=pch[classnames[[l]][j]], type="l")
304 colvect <- append(colvect,cols[classnames[[l]][j]])
305 }
306 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col=colvect, lty=lty, pch=pch)
307 }
308 }
309 }#end if length >2
310 if (length(phenoDataClass)==2){
311 k <- 1
312 l <- 2
313
314 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k],"vs",phenoDataClass[l], sep=""), xlab="Retention Time (min)", ylab="TIC")
315 colvect <- NULL
316 for (j in 1:length(classnames[[k]])) {
317 tic <- TIC[[classnames[[k]][j]]]
318 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
319 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
320 colvect <- append(colvect,cols[classnames[[k]][j]])
321 }
322 for (j in 1:length(classnames[[l]])) {
323 # i <- class2names[j]
324 tic <- TIC[[classnames[[l]][j]]]
325 points(tic[,1]/60, -tic[,2], col=cols[classnames[[l]][j]], pch=pch[classnames[[l]][j]], type="l")
326 colvect <- append(colvect,cols[classnames[[l]][j]])
327 }
328 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col=colvect, lty=lty, pch=pch)
329
330 }#end length ==2
331
332 #case where only one class
333 if (length(phenoDataClass)==1){
334 k <- 1
335 ylim <- range(sapply(TIC, function(x) range(x[,2])))
336
337 plot(0, 0, type="n", xlim=xlim/60, ylim=ylim, main=paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k], sep=""), xlab="Retention Time (min)", ylab="TIC")
338 colvect <- NULL
339 for (j in 1:length(classnames[[k]])) {
340 tic <- TIC[[classnames[[k]][j]]]
341 # points(tic[,1]/60, tic[,2], col=cols[i], pch=pch[i], type="l")
342 points(tic[,1]/60, tic[,2], col=cols[classnames[[k]][j]], pch=pch[classnames[[k]][j]], type="l")
343 colvect <- append(colvect,cols[classnames[[k]][j]])
344 }
345
346 legend("topright",paste(basename(files[c(classnames[[k]])])), col=colvect, lty=lty, pch=pch)
347
348 }#end length ==1
349
350 dev.off() #pdf(pdfname,w=16,h=10)
351
352 invisible(TIC)
353 }
354
355 150
356 151
357 # Get the polarities from all the samples of a condition 152 # Get the polarities from all the samples of a condition
358 #@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM 153 #@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM
359 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM 154 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM