Skip to content

Commit

Permalink
Added my_plotEICs.R function. It implements a bug fix that addresses …
Browse files Browse the repository at this point in the history
…an error found in the color labels of CAMERA's plotEICs plots.
  • Loading branch information
wilmelz committed Apr 20, 2020
1 parent 2f98ec5 commit f64b5f9
Show file tree
Hide file tree
Showing 3 changed files with 184 additions and 9 deletions.
29 changes: 29 additions & 0 deletions R/file-location-workaround.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# We are currently receiving an error when running ParseCAMERA with gen.plots = TRUE
#Munge Data Modules
#if positive, run:
# ParseCAMERA(from.table = "Annotated", to.table = "output_parsed", CAMERA.obj = "anposGa")
#
# Error in .local(object, ...) :
# xcmsSource: file not found: E:\ORD Backup\Users\jmosley\Desktop\WLSSD 2014 Livers\mzML/Pooled QCs/Pooled_QC_Pos_10.mzML

#To get around this, set the variables in the script below and run *AFTER* running the InitWorkflow() module. Then execute the remaining script code to bypass the bug.

##Set these variables to be equivalent on your local system
data_dir <- "O:/Public/jmosl01/LUMA-Data"


##Execute this code only after setting the above variables for your local system
new_filepaths <- list.files(path = data_dir, recursive = TRUE, full.names = TRUE)

if(BLANK) new_filepaths <- new_filepaths[grepl("Blanks", new_filepaths)]
if(!BLANK) new_filepaths <- new_filepaths[!grepl("Blanks", new_filepaths)]

#If Positive mode, run this:
new_filepaths <- new_filepaths[grepl(ion.id[1],new_filepaths)]

anposGa@xcmsSet@filepaths <- new_filepaths

#If Negative mode, run this:
new_filepaths <- new_filepaths[grepl(ion.id[2],new_filepaths)]

annegGa@xcmsSet@filepaths <- new_filepaths
136 changes: 136 additions & 0 deletions R/my_plotEICs.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
#
# This function is a modified version of CAMERA's plotEICs. It implements a bug fix that
# addresses an error with the color labels in the plots. Bug fix was first proposed by
# Jonathan Mosley.
#
my_plotEICs <- function (object, pspec = 1:length(object@pspectra), maxlabel = 0,
sleep = 0, method = "bin")
{
smpls <- unique(object@psSamples[pspec])
xeic <- new("xcmsEIC")
xeic@rtrange <- matrix(nrow = length(pspec), ncol = 2)
xeic@mzrange <- matrix(nrow = length(pspec), ncol = 2)
pcpos <- 1
rtmargin <- 1
for (a in seq(along = smpls)) {
xraw <- xcmsRaw(object@xcmsSet@filepaths[smpls[a]],
profmethod = method)
pspecS <- pspec[which(object@psSamples[pspec] == smpls[a])]
peaks <- CAMERA:::getPeaks(object@xcmsSet, smpls[a])
eic <- lapply(pspecS, function(pc) {
pidx <- object@pspectra[[pc]]
pks <- peaks[pidx, , drop = FALSE]
gks <- object@groupInfo[pidx, , drop = FALSE]
nap <- which(is.na(pks[, 1]))
pks[nap, ] <- cbind(gks[nap, c(1:6), drop = FALSE],
matrix(nrow = length(nap), ncol = 5, 0))
bbox <- c(rtmin = min(pks[, "rtmin"]) - rtmargin,
rtmax = max(pks[, "rtmax"]) + rtmargin,
mzmin = min(pks[,"mzmin"]),
mzmax = max(pks[, "mzmax"]))
eic <- xcms:::getEIC(xraw, rtrange = pks[, c("rtmin","rtmax"), drop = FALSE],
mzrange = pks[, c("mzmin","mzmax"), drop = FALSE])
xeic@rtrange[pcpos, ] <<- bbox[c("rtmin", "rtmax")]
xeic@mzrange[pcpos, ] <<- bbox[c("mzmin", "mzmax")]
cat("-->", pcpos, "\n")
pcpos <<- pcpos + 1
eic@eic[[1]]
})
xeic@eic <- c(xeic@eic, eic)
}
names(xeic@eic) <- paste("Pseudospectrum ", pspec, sep = "")
if (maxlabel > 0) {
col <- rainbow(maxlabel)
}
else {
col <- c()
}
for (ps in seq(along = pspec)) {
EIC <- xeic@eic[[ps]]
pidx <- object@pspectra[[pspec[ps]]]
peaks <- CAMERA:::getPeaks(object@xcmsSet, object@psSamples[pspec[ps]])[pidx,, drop = FALSE]
grps <- object@groupInfo[pidx, ]
nap <- which(is.na(peaks[, 1]))
naps <- rep(FALSE, nrow(peaks))
if (length(nap) > 0) {
naps[nap] <- TRUE
peaks[nap, ] <- cbind(grps[nap, c(1:6), drop = FALSE],
matrix(nrow = length(nap), ncol = 5, 0))
}
main <- paste("Pseudospectrum ", pspec[ps], sep = "")
neics <- length(pidx)
lmaxlabel <- min(maxlabel, neics)
eicidx <- 1:neics
maxint <- numeric(length(eicidx))
for (j in eicidx) {
maxint[j] <- max(EIC[[j]][, "intensity"])
}
o <- order(maxint, decreasing = TRUE)
rt <- xeic@rtrange[ps, ]
rt.min <- round(mean(peaks[, "rtmin"]), digits = 3)
rt.med <- round(peaks[o[1], "rt"], digits = 3)
rt.max <- round(mean(peaks[, "rtmax"]), digits = 3)
plot(0, 0, type = "n", xlim = rt, ylim = c(0, max(maxint)),
xaxs = "i", xlab = "Retention Time (seconds)",
ylab = "Intensity", main = paste("Extracted Ion Chromatograms for ",
main, "\nTime: From", rt.min, "to", rt.max,
", mean", rt.med))
lcol <- rgb(0.6, 0.6, 0.6)
lcol <- c(col, rep(lcol, max(nrow(peaks) - maxlabel,0)))
cnt <- 1
for (j in eicidx[o]) {
pts <- xeic@eic[[ps]][[j]]
points(pts, type = "l", col = lcol[cnt])
peakrange <- peaks[, c("rtmin", "rtmax"), drop = FALSE]
ptsidx <- pts[, "rt"] >= peakrange[j, 1] & pts[,"rt"] <= peakrange[j, 2]
if (naps[j]) {
points(pts[ptsidx, ], type = "l", col = col[cnt],lwd = 1.3, lty = 3)
}
else {
points(pts[ptsidx, ], type = "l", col = col[cnt],lwd = 1.3)
}
cnt <- cnt + 1
}
pspectrum <- getpspectra(object, grp = pspec[ps])
mz <- pspectrum[o, "mz"]
if (lmaxlabel > 0 & "adduct" %in% colnames(pspectrum)) {
adduct <- sub("^ ", "", pspectrum[o, "adduct"])
mass <- sapply(strsplit(adduct, " "), function(x) {
x[2]
})
adduct <- sapply(strsplit(adduct, " "), function(x) {
x[1]
})
umass <- unique(na.omit(mass[1:maxlabel]))
adduct[is.na(adduct)] <- ""
test <- vector("list", length = length(mz))
mz <- format(pspectrum[o[1:lmaxlabel], "mz"],
digits = 5)
if (length(umass) > 0) {
for (i in 1:length(umass)) {
ini <- which(mass == umass[i])
for (ii in 1:length(ini)) {
firstpart <- strsplit(adduct[ini[ii]],"M")[[1]][1]
secondpart <- strsplit(adduct[ini[ii]],"M")[[1]][2]
masspart <- mz[ini[ii]]
test[[ini[ii]]] <- substitute(paste(masspart," ", firstpart, M[i], secondpart),
list(firstpart = firstpart,
i = i,
secondpart = secondpart,
masspart = masspart))
}
}
}
for (i in seq(along = lmaxlabel)) {
if (is.null(test[[i]])) {
test[[i]] <- mz[i]
}
}
leg <- as.expression(test[1:lmaxlabel])
legend("topright", legend = leg, col = lcol,lty = 1)
}
if (sleep > 0) {
Sys.sleep(sleep)
}
}
}
28 changes: 19 additions & 9 deletions R/plotting-functions.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
source("my_plotEICs.R")
#' @title Loop-based EIC Plotter for metabolite groups
#'
#' @export
Expand Down Expand Up @@ -33,7 +34,6 @@ plot_metgroup = function(CAMERA.obj, Sample.df, Peak.list, center, BLANK, gen.pl
if (missing(maxlabel))
maxlabel = 10


# Change the psspectra list in the xsAnnotate object to be the unique list of metabolite_groups
X <- split(Peak.list$EIC_ID, as.numeric(Peak.list$metabolite_group))
names(X) <- sort(unique(Peak.list$metabolite_group))
Expand Down Expand Up @@ -93,21 +93,26 @@ plot_metgroup = function(CAMERA.obj, Sample.df, Peak.list, center, BLANK, gen.pl
test.mat <- as.matrix(t(my.mat))
dimnames(test.mat) <- list(colnames(my.df[, sample.cols]), my.df$EIC_ID)
plot(c(0, 1), c(0, 1), ann = F, bty = "n", type = "n", xaxt = "n", yaxt = "n")
text(x = 0.5, y = 0.5, paste("Dendrograms can't be plotted for two features.\n", "Use the correlation plot to the left to \n",
"Tell if two features are from the same metabolite"), cex = 1.6, col = "black")
text(x = 0.5, y = 0.5, paste("Dendrograms can't be plotted for two features.\n", "Use the correlation plot to the left to tell \n",
"if two features are from the same metabolite"), cex = 1.6, col = "black")
res2 <- rcorr(test.mat)
M <- res2$r
P <- res2$P
# Insignificant correlations are leaved blank
order.hc2 <- corrplot::corrMatOrder(M, order = "hclust", hclust.method = "complete")
col.new = rainbow(maxlabel)[order.hc2]
corrplot(M, type = "upper", order = "hclust", hclust.method = "complete", p.mat = P, sig.level = 0.01, insig = "blank",
tl.col = col.new)
tl.col = col.new, cl.ratio = 0.15, cl.align.text = "l", cl.cex = 0.6, cl.offset = 0.2)

# Plots the EICs and Pseudo-spectra for each metabolite group containing more than one feature in a for loop
EIC.plots <- plotEICs(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) ## Plot the EICs
EIC.plots <- my_plotEICs(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) ## Plot the EICs

a1 <- new_CAMERA.obj@pspectra[[get.mg[i]]]
b1 <- new_CAMERA.obj@groupInfo[a1,"mz"]
mz1 <- min(b1) * 0.80
mz2 <- max(b1) * 1.20

plotPsSpectrum(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) #Plot the Mass Spectra
plotPsSpectrum(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0, cex.main = 0.8, cexMulti = 1.0, mzrange = c(mz1,mz2)) #Plot the Mass Spectra

} else {
my.mat <- as.matrix(my.df[, grep(paste(QC.id, paste(strsplit(sexes, "(?<=.[_])", perl = TRUE),
Expand All @@ -129,11 +134,16 @@ plot_metgroup = function(CAMERA.obj, Sample.df, Peak.list, center, BLANK, gen.pl
order.hc2 <- corrplot::corrMatOrder(M, order = "hclust", hclust.method = "complete")
col.new = rainbow(maxlabel)[order.hc2]
corrplot(M, type = "upper", order = "hclust", hclust.method = "complete", p.mat = P, sig.level = 0.01, insig = "blank",
tl.col = col.new)
tl.col = col.new, cl.ratio = 0.15, cl.align.text = "l", cl.cex = 0.6, cl.offset = 0.2)
# Plots the EICs and Pseudo-spectra for each metabolite group containing more than one feature in a for loop
EIC.plots <- plotEICs(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) ## Plot the EICs
EIC.plots <- my_plotEICs(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) ## Plot the EICs

a1 <- new_CAMERA.obj@pspectra[[get.mg[i]]]
b1 <- new_CAMERA.obj@groupInfo[a1,"mz"]
mz1 <- min(b1) * 0.80
mz2 <- max(b1) * 1.20

plotPsSpectrum(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0) #Plot the Mass Spectra
plotPsSpectrum(new_CAMERA.obj, pspec = get.mg[i], maxlabel = maxlabel, sleep = 0, cex.main = 0.8, cexMulti = 1.0, mzrange = c(mz1, mz2)) #Plot the Mass Spectra
}
}
}
Expand Down

4 comments on commit f64b5f9

@jmosl01
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Wilson,

I am going to rollback any code changes related to the my_plotEICs.R function, effectively doing away with the addition of this custom CAMERA function to LUMA. Instead, I will submit a bug fix to CAMERA for the known issue in CAMERA "Bug in labelling and/or legend of plotEICs" (sneumann/CAMERA#2)

@wilmelz
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Jonathan,

Sounds good. Thanks for the update.

@jmosl01
Copy link
Owner

@jmosl01 jmosl01 commented on f64b5f9 May 4, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Wilson,

I finished the bug fix and submitted a pull request to @sneumann: sneumann/CAMERA#55 to fix sneumann/CAMERA#2.

@wilmelz
Copy link
Collaborator Author

@wilmelz wilmelz commented on f64b5f9 May 4, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Jonathan,

That's great to know. Let's hope @sneumann is more responsive this time around. Thanks for the update.

Please sign in to comment.