Mercurial > repos > yguitton > metams_rungc

comparison lib_metams.r @ 5:b8d4129dd2a6 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/metaMS commit c7a518686137f6d62b7415715152e8d5a9953ed7
author yguitton
date Fri, 06 Sep 2019 06:09:10 -0400
parents c10824185547
children d1ce2634135f
comparison
equal deleted inserted replaced
4:c10824185547 5:b8d4129dd2a6
160 } 160 }
161 161
162 ##ADDITIONS FROM Y. Guitton 162 ##ADDITIONS FROM Y. Guitton
163 getBPC <- function(file,rtcor=NULL, ...) { 163 getBPC <- function(file,rtcor=NULL, ...) {
164 object <- xcmsRaw(file) 164 object <- xcmsRaw(file)
165 sel <- profRange(object, ...) 165 sel <- profRange(object, ...)
166 cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE])) 166 cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE]))
167 } 167 }
168 168
169 getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) { 169 getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) {
170 require(xcms) 170 require(xcms)
171 171
172 #create sampleMetadata, get sampleMetadata and class 172 #create sampleMetadata, get sampleMetadata and class
173 if(!is.null(xset)) { 173 if(!is.null(xset)) {
174 #When files come from XCMS3 directly before metaMS 174 #When files come from XCMS3 directly before metaMS
175 sampleMetadata <- xset@phenoData 175 sampleMetadata <- xset@phenoData
176 } else { 176 } else {
182 for (i in 1:length(class)){ 182 for (i in 1:length(class)){
183 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i]) 183 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i])
184 } 184 }
185 185
186 N <- dim(sampleMetadata)[1] 186 N <- dim(sampleMetadata)[1]
187 TIC <- vector("list",N) 187 BPC <- vector("list",N)
188 188
189 for (j in 1:N) { 189 for (j in 1:N) {
190 TIC[[j]] <- getBPC(files[j]) 190 BPC[[j]] <- getBPC(files[j])
191 #good for raw 191 #good for raw
192 # seems strange for corrected 192 # seems strange for corrected
193 #errors if scanrange used in xcmsSetgeneration 193 #errors if scanrange used in xcmsSetgeneration
194 if (!is.null(xcmsSet) && rt == "corrected"){ 194 if (!is.null(xcmsSet) && rt == "corrected"){
195 rtcor <- xcmsSet@rt$corrected[[j]] 195 rtcor <- xcmsSet@rt$corrected[[j]]
196 }else{ 196 }else{
197 rtcor <- NULL 197 rtcor <- NULL
198 } 198 }
199 TIC[[j]] <- getBPC(files[j],rtcor=rtcor) 199 BPC[[j]] <- getBPC(files[j],rtcor=rtcor)
200 } 200 }
201 201
202 pdf(pdfname,w=16,h=10) 202 pdf(pdfname,w=16,h=10)
203 cols <- rainbow(N) 203 cols <- rainbow(N)
204 lty = 1:N 204 lty = 1:N
205 pch = 1:N 205 pch = 1:N
206 #search for max x and max y in BPCs 206 #search for max x and max y in BPCs
207 xlim = range(sapply(TIC, function(x) range(x[,1]))) 207
208 ylim = range(sapply(TIC, function(x) range(x[,2]))) 208 xlim = range(sapply(BPC, function(x) range(x[,1])))
209 ylim = range(sapply(BPC, function(x) range(x[,2])))
210
209 ylim = c(-ylim[2], ylim[2]) 211 ylim = c(-ylim[2], ylim[2])
210 212
211 ##plot start 213 ##plot start
212 if (length(class)>2){ 214 if (length(class)>2){
213 for (k in 1:(length(class)-1)){ 215 for (k in 1:(length(class)-1)){
214 for (l in (k+1):length(class)){ 216 for (l in (k+1):length(class)){
215 cat(paste(class[k],"vs",class[l],sep=" ","\n")) 217 cat(paste(class[k],"vs",class[l],sep=" ","\n"))
216 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") 218 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
217 colvect<-NULL 219 colvect<-NULL
218 for (j in 1:length(classnames[[k]])) { 220 for (j in 1:length(classnames[[k]])) {
219 tic <- TIC[[classnames[[k]][j]]] 221 bpc <- BPC[[classnames[[k]][j]]]
220 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 222 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
221 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 223 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
222 colvect<-append(colvect,cols[classnames[[k]][j]]) 224 colvect<-append(colvect,cols[classnames[[k]][j]])
223 } 225 }
224 for (j in 1:length(classnames[[l]])) { 226 for (j in 1:length(classnames[[l]])) {
225 # i=class2names[j] 227 # i=class2names[j]
226 tic <- TIC[[classnames[[l]][j]]] 228 bpc <- BPC[[classnames[[l]][j]]]
227 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 229 points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
228 colvect<-append(colvect,cols[classnames[[l]][j]]) 230 colvect<-append(colvect,cols[classnames[[l]][j]])
229 } 231 }
230 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) 232 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
231 } 233 }
232 } 234 }
237 l=2 239 l=2
238 colvect<-NULL 240 colvect<-NULL
239 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") 241 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
240 242
241 for (j in 1:length(classnames[[k]])) { 243 for (j in 1:length(classnames[[k]])) {
242 tic <- TIC[[classnames[[k]][j]]] 244 bpc <- BPC[[classnames[[k]][j]]]
243 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 245 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
244 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 246 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
245 colvect<-append(colvect,cols[classnames[[k]][j]]) 247 colvect<-append(colvect,cols[classnames[[k]][j]])
246 } 248 }
247 for (j in 1:length(classnames[[l]])) { 249 for (j in 1:length(classnames[[l]])) {
248 # i=class2names[j] 250 # i=class2names[j]
249 tic <- TIC[[classnames[[l]][j]]] 251 bpc <- BPC[[classnames[[l]][j]]]
250 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 252 points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
251 colvect<-append(colvect,cols[classnames[[l]][j]]) 253 colvect<-append(colvect,cols[classnames[[l]][j]])
252 } 254 }
253 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) 255 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
254 }#end length ==2 256 }#end length ==2
255 257
256 if (length(class)==1){ 258 if (length(class)==1){
257 k=1 259 k=1
258 ylim = range(sapply(TIC, function(x) range(x[,2]))) 260
261 ylim = range(sapply(BPC, function(x) range(x[,2])))
262
259 colvect<-NULL 263 colvect<-NULL
260 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC") 264 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
261 265
262 for (j in 1:length(classnames[[k]])) { 266 for (j in 1:length(classnames[[k]])) {
263 tic <- TIC[[classnames[[k]][j]]] 267 bpc <- BPC[[classnames[[k]][j]]]
264 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 268 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
265 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 269 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
266 colvect<-append(colvect,cols[classnames[[k]][j]]) 270 colvect<-append(colvect,cols[classnames[[k]][j]])
267 } 271 }
268 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) 272 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch)
269 }#end length ==1 273 }#end length ==1
270 dev.off() 274 dev.off()
295 299
296 N <- dim(sampleMetadata)[1] 300 N <- dim(sampleMetadata)[1]
297 TIC <- vector("list",N) 301 TIC <- vector("list",N)
298 302
299 for (i in 1:N) { 303 for (i in 1:N) {
300 cat(files[i],"\n")
301 if (!is.null(xcmsSet) && rt == "corrected") 304 if (!is.null(xcmsSet) && rt == "corrected")
302 rtcor <- xcmsSet@rt$corrected[[i]] 305 rtcor <- xcmsSet@rt$corrected[[i]]
303 else 306 else
304 rtcor <- NULL 307 rtcor <- NULL
305 TIC[[i]] <- getTIC(files[i],rtcor=rtcor) 308 TIC[[i]] <- getTIC(files[i],rtcor=rtcor)
316 319
317 ##plot start 320 ##plot start
318 if (length(class)>2){ 321 if (length(class)>2){
319 for (k in 1:(length(class)-1)){ 322 for (k in 1:(length(class)-1)){
320 for (l in (k+1):length(class)){ 323 for (l in (k+1):length(class)){
321 print(paste(class[k],"vs",class[l],sep=" ")) 324 cat(paste(class[k],"vs",class[l],"\n",sep=" "))
322 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") 325 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
323 colvect<-NULL 326 colvect<-NULL
324 for (j in 1:length(classnames[[k]])) { 327 for (j in 1:length(classnames[[k]])) {
325 tic <- TIC[[classnames[[k]][j]]] 328 tic <- TIC[[classnames[[k]][j]]]
326 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 329 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")