rformassspectrometry / Spectra

Low level infrastructure to handle MS spectra
https://rformassspectrometry.github.io/Spectra/
34 stars 24 forks source link

How to plot an EIC with the Spectra package? #241

Closed wkumler closed 2 years ago

wkumler commented 2 years ago

Hi all, neat work on this package so far and congrats on the publication in Metabolites! I have a dumb question about the package - is it possible to plot an EIC for a given mass range? It seems like many of the use cases in the tutorial are MSMS focused and I'm hoping to see how the new backends compare with existing packages. I'm basically trying to replicate the code below but with Spectra instead of RaMS or MSnbase:

fls <- msdata::proteomics(full.names=TRUE)[4]

library(RaMS)
msdata <- RaMS::grabMSdata(fls)
plot(msdata$MS1[mz%between%pmppm(504.3279, 5), c("rt", "int")], type="l")

library(MSnbase)
ms <- MSnbase::readMSData(fls, mode = "onDisk", msLevel. = 1)
plot(chromatogram(ms, mz=pmppm(504.3279, 5)))

So far I've tried a couple different methods but none seem to extract an EIC as I'm expecting:

library(Spectra)
sps_mzr <- Spectra(fls, source=MsBackendMzR())
sps_mzr |> filterMzRange(c(504.32537, 504.33042)) |> plotSpectra()
# Makes a plot for every MS1 scan

sps_df <- setBackend(v, MsBackendDataFrame())
sps_df |> filterMzRange(c(504.32537, 504.33042)) |> plotSpectra()
# Same as above, one plot per MS1 scan

chromatogram(sps_mzr, mz=c(504.32537, 504.33042))
chromatogram(sps_df, mz=c(504.32537, 504.33042))
# Unable to convert Spectra to chromatogram

as.data.frame(sps_mzr)
as.data.frame(sps_df)
# Unable to convert Spectra to data frame

I'm able to hack together a chromatogram with the below code but this doesn't feel like it obeys the spirit of Spectra at all:

sfs_filtered <- Spectra(fls, source=MsBackendMzR()) |> 
  filterMzRange(c(504.32537, 504.33042)) |>
  filterMsLevel(1)
mzs <- mz(sfs_filtered)
ints <- intensity(sfs_filtered)
rts <- rtime(sfs_filtered)
eic <- lapply(seq_along(rts), function(i){
  if(length(mzs[[i]])){
    cbind(rt=rts[[i]], mz=mzs[[i]], int=ints[[i]])
  }
}) |> do.call(what = "rbind")
plot(eic[,"rt"], eic[,"int"], type="l")

Is there a better way to assemble a chromatogram for a given mass?

jorainer commented 2 years ago

sorry - pushed enter too early - this is the edited/completed example:

Unfortunately there is no such functionality yet in Spectra. We've started working on the Chromatograms package, but unfortunately we didn't find the time recently to work on it (but it's high on the TODO list...).

Maybe a slightly improved version over your code above:

## need to aggregate intensities of all peaks in a spectrum
## Simple function that takes a peaks matrix as input `x` and aggregates 
## (sums) the intensities and returns a single row peaks matrix. For spectra
## with no peaks available, NA will be returned as intensity
.sum_intensities <- function(x, ...) {
    if (nrow(x)) {
        cbind(mz = NA_real_, intensity = sum(x[, "intensity"], na.rm = TRUE))
    } else cbind(mz = NA_real_, intensity = NA_real_)
}
sfs_agg <- addProcessing(sfs_filtered, .sum_intensities)
eic <- cbind(rtime(sfs_agg), unlist(intensity(sfs_agg), use.names = FALSE))

That version will do automatically per-file parallel processing and loading data chunk-wise into memory hence keeping a low memory footprint.

wkumler commented 2 years ago

I see, thank you! The new version you suggest is indeed more effective. Looking forward to the eventual release of Chromatograms!