Skip to main content

Topic: replacement for getEIC in XCMS3 (Read 908 times) previous topic - next topic

  • krista
  • [*]
replacement for getEIC in XCMS3
I am in the process of updating to use XCMS3. So far, I like XCMS3 and most of the transition to XCMS3 has been easy in terms of how to set peak picking, retention time correction, and defining correspondence across a set of files (the old 'grouping'). However, one thing I miss is the ability to plot all of our mzRT features into a PDF that we could use to assess how well our parameters were working. I find multiple examples to extract a pre-determined mz / RT combination in the vignettes provided with XCMS3. Therefore, if I know the mz and retention time for a single peak, I can easily plot it. However, what I cannot figure out is out to get all the peaks across my set of samples.

I used to use the code listed below. In the example, xgF is the data after peak picking, group, retention time correction, and fillPeaks.

Any thoughts?
Cheers,
Krista


Quote
xs.fill <- xgF
#
xeic.raw <- getEIC(xs.fill, rt = "raw", groupidx= 1:nrow(xs.fill@groups))
xeic.corrected <- getEIC(xs.fill, rt = "corrected", groupidx= 1:nrow(xs.fill@groups))

#colorList should be as long as the number of possible levels in the dataset...
colorList <- c("deepskyblue","magenta","forestgreen","darkorchid","firebrick", "gold")

#plot first feature from diffreport (not sorted after pval!!)
cairo_pdf(file = fName_peaksPicked,onefile=TRUE)
for (i in 1:nrow(xs.fill@groups)){
  #for (i in 1:5){
 
  par(mfrow= c(2,1))
  #rt row data
  plot(xeic.raw, xs.fill,groupidx=i,col=colorList)
  #rt corrected data
  plot(xeic.corrected, xs.fill,groupidx=i,col=colorList)
 
}
dev.off()

My sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base

other attached packages:
[1] CAMERA_1.36.0 xcms_3.2.0 MSnbase_2.6.4 ProtGenerics_1.12.0 mzR_2.14.0
[6] Rcpp_0.12.19 BiocParallel_1.14.2 Biobase_2.40.0 BiocGenerics_0.26.0

loaded via a namespace (and not attached):
 [1] lattice_0.20-35 snow_0.4-3 digest_0.6.18 foreach_1.4.4
 [5] plyr_1.8.4 backports_1.1.2 acepack_1.4.1 mzID_1.18.0
 [9] stats4_3.5.1 ggplot2_3.0.0 BiocInstaller_1.30.0 pillar_1.3.0
[13] zlibbioc_1.26.0 rlang_0.2.2 lazyeval_0.2.1 rstudioapi_0.8
[17] data.table_1.11.8 S4Vectors_0.18.3 rpart_4.1-13 Matrix_1.2-14
[21] checkmate_1.8.5 preprocessCore_1.42.0 splines_3.5.1 stringr_1.3.1
[25] foreign_0.8-71 htmlwidgets_1.3 igraph_1.2.2 munsell_0.5.0
[29] compiler_3.5.1 pkgconfig_2.0.2 base64enc_0.1-3 multtest_2.36.0
[33] pcaMethods_1.72.0 htmltools_0.3.6 nnet_7.3-12 tibble_1.4.2
[37] gridExtra_2.3 htmlTable_1.12 RANN_2.6 Hmisc_4.1-1
[41] IRanges_2.14.12 codetools_0.2-15 XML_3.98-1.16 crayon_1.3.4
[45] MASS_7.3-50 grid_3.5.1 MassSpecWavelet_1.46.0 RBGL_1.56.0
[49] gtable_0.2.0 affy_1.58.0 magrittr_1.5 scales_1.0.0
[53] graph_1.58.2 stringi_1.2.4 impute_1.54.0 affyio_1.50.0
[57] doParallel_1.0.14 limma_3.36.5 latticeExtra_0.6-28 Formula_1.2-3
[61] RColorBrewer_1.1-2 iterators_1.0.10 tools_3.5.1 survival_2.42-6
[65] yaml_2.2.0 colorspace_1.3-2 cluster_2.0.7-1 vsn_3.48.1
[69] MALDIquant_1.18 knitr_1.20



  • CoreyG
  • [*]
Re: replacement for getEIC in XCMS3
Reply #1
Hi Krista,

I'm getting used to XCMS3, but I hope I can help you anyway.

You might be getting unstuck due to differences in XCMSnExp and xcmsSet classes. Using the 'as' function, you can convert an XCMSnExp object to an xcmsSet object.

I think all of your code will work if you replace the first line with:
xs.fill <- as(xgF,"xcmsSet")

Please let me (and everyone else) know if this works for you :)
  • Baker Heart and Diabetes Institute

Re: replacement for getEIC in XCMS3
Reply #2
Dear Krista,

I realized that that's something apparently missing in xcms version 3 at present. I will think of a solution and post it. Thanks for the feedback.

cheers, jo
  • Eurac Research, Bolzano, Italy

Re: replacement for getEIC in XCMS3
Reply #3
For now you can use the function below (just copy-paste into R). I will add it also to xcms and update the vignette to describe your use case.

Code: [Select]
featureChromatograms <- function(x, expandRt = 0, aggregationFun = "max",
                                 ...) {
    if (!hasFeatures(x))
        stop("No feature definitions present. Please run first 'groupChromPeaks'")
    pks <- chromPeaks(x)
    mat <- do.call(rbind, lapply(featureDefinitions(x)$peakidx, function(z) {
        pks_current <- pks[z, , drop = FALSE]
        c(range(pks_current[, c("rtmin", "rtmax")]),
          range(pks_current[, c("mzmin", "mzmax")]))
    }))
    colnames(mat) <- c("rtmin", "rtmax", "mzmin", "mzmax")
    chromatogram(x, rt = mat[, 1:2], mz = mat[, 3:4],
                 aggregationFun = aggregationFun, ...)
}

This function takes an XCMSnExp object, extracts chromatograms for each feature and returns it as an Chromatograms object. The Chromatograms arranges the XIC in a two-dimensional matrix, columns being the individual samples. Each row contains the XIC for one feature. It is thus straight forward to plot the data for one specific feature. The example code below illustrates it on the famous faahko data set:

Code: [Select]
library(xcms)
library(faahKO)
faahko_3_files <- c(system.file('cdf/KO/ko15.CDF', package = "faahKO"),
                    system.file('cdf/KO/ko16.CDF', package = "faahKO"),
                    system.file('cdf/KO/ko18.CDF', package = "faahKO"))

## Do a quick and dirty preprocessing
od <- readMSData(faahko_3_files, mode = "onDisk")
od <- findChromPeaks(od, param = CentWaveParam(peakwidth = c(30, 80),
        noise = 1000))
od <- adjustRtime(od, param = ObiwarpParam(binSize = 0.6))
od <- groupChromPeaks(od,
        param = PeakDensityParam(minFraction = 0.8, sampleGroups = rep(1, 3)))

## Extract ion chromatograms for each feature (need to post the code above or use xcms version > 3.3.4
chrs <- featureChromatograms(od)

## Extract also XICs without adjusted retention times
chrs_raw <- featureChromatograms(od, adjustedRtime = FALSE)

## Plot the XIC for the first feature using different colors for each file
par(mfrow = c(1, 2))
plot(chrs[1, ], col = c("red", "green", "blue"))
plot(chrs_raw[1, ], col = c("red", "green", "blue"))


In your use case you could simply make then a for loop over nrow(chrs) to plot the data for each feature.

Note: you could also use the highlightChromPeaks function to indicate the identified peaks in the individual XICs:

Code: [Select]
plot(chrs[1, ], col = c("red", "green", "blue"))
highlightChromPeaks(od, rt = range(lapply(chrs[1, ], rtime)), mz = range(lapply(chrs[1, ], mz)),
    border = c("red", "green", "blue"))



Hope that helps.
cheers, jo
  • Eurac Research, Bolzano, Italy

  • krista
  • [*]
Re: replacement for getEIC in XCMS3
Reply #4
Hello,

CoreyG:  I have swapped between the XCMSnExp and xcmsSet using the as function. I used that to have CAMERA work for isotopes and adducts. It 'sort-of' works for the code that I posted, although I get an error message during the plotting that I cannot quite sort out:
Error in legend(rtrange[i, 2], max(maxint), legtext, col = unique(col[sampidx]), : 
 'legend' is of length 0

I am not sure where I would set the 'legend' parameters, so I am not sure how to get around this.

Jo: thanks for thinking about how to implement this in XCMS3. I used to use the PDF generated during the centWave peakPicking parameters, but at some point (even prior to XCMS3) that went away. Either is fine with me. I just like to see all/most of my features to see how I am doing beyond a pre-defined set of mz/RT values.

Thanks in advance. I look forward to moving forward with the transition to XCMS3.

Cheers,
Krista


  • krista
  • [*]
Re: replacement for getEIC in XCMS3
Reply #5
Hi Jo,

Thanks for the code. This really helps me out. I will take a look at it and see how things work. Either way I will make sure to post a reply later so that we can see how this works.

Cheers,
Krista

  • krista
  • [*]
Re: replacement for getEIC in XCMS3
Reply #6
Jo,
This is perfect. Exactly what I needed. The plots are helpful, and the code to loop over each row in chrs generates the PDF I want. For others, here's what I have:

Code: [Select]
library(xcms)
library(RColorBrewer)

#setup the name of the PDF file that I want:
fName_peaksPicked <- 'yourFilenamehere.pdf'

#preset the colors, use ncol to figure out how many samples I have (and hence # of colors needed)nColors <- ncol(chrs)
useColors <- brewer.pal(nColors,'YlGnBu')

cairo_pdf(file = fName_peaksPicked,onefile=TRUE)
for (i in 1:nrow(chrs)){
    par(mfrow = c(1, 2))
    plot(chrs[i, ], col = useColors)
    plot(chrs_raw[i, ], col = useColors)
   
    # You could also use the highlightChromPeaks function to indicate the identified
    #peaks in the individual XICs:
    if (TRUE) {
      plot(chrs[i, ], col = useColors)
      highlightChromPeaks(xgF, rt = range(lapply(chrs[i, ], rtime)),
                        mz = range(lapply(chrs[i, ], mz)),
                        border = useColors)
    }
 }
dev.off()



Thanks again!
Krista

  • krista
  • [*]
Re: replacement for getEIC in XCMS3
Reply #7
Ack..I just realized I only posted half of the code in my previous post. The first part uses Jo's function as follows:

Code: [Select]
featureChromatograms <- function(x, expandRt = 0, aggregationFun = "max",
                                 ...) {
    if (!hasFeatures(x))
        stop("No feature definitions present. Please run first 'groupChromPeaks'")
    pks <- chromPeaks(x)
    mat <- do.call(rbind, lapply(featureDefinitions(x)$peakidx, function(z) {
        pks_current <- pks[z, , drop = FALSE]
        c(range(pks_current[, c("rtmin", "rtmax")]),
          range(pks_current[, c("mzmin", "mzmax")]))
    }))
    colnames(mat) <- c("rtmin", "rtmax", "mzmin", "mzmax")
    chromatogram(x, rt = mat[, 1:2], mz = mat[, 3:4],
                 aggregationFun = aggregationFun, ...)
}


#now use the code from Jo:
#extract the EICs from the grouped XCMSnExp object:
chrs <- featureChromatograms(xgF)

## Extract also XICs without adjusted retention times
chrs_raw <- featureChromatograms(xgF, adjustedRtime = FALSE)

Then, I use the code I posted earlier today to do the actual plotting (see above). I hope others can take advantage of this.
Krista

Re: replacement for getEIC in XCMS3
Reply #8
Just to update: we've added the function to the xcms codebase. It will be officially available with the next Bioconductor release (3.8) which will be released end of October 2018.

cheers, jo
  • Eurac Research, Bolzano, Italy

  • krista
  • [*]
Re: replacement for getEIC in XCMS3
Reply #9
Excellent, I will keep an eye out for the new version of XCMS3.
Thanks for all your help.
Krista