MS data transformation timings
Below, I’ve pulled out the “transformation” step from each of the eight methods above that gets the MS data out of the mzML documents and into an R-readable structure. For the first four (MSnExp, onDiskMSnExp, Spectra, and RaMS) this creates an object in the R environment that is wiped each time R is restarted, so they must be recreated every time. For the last four, files are created on the hard drive that are robust across R sessions and only need to be created once.
Code
msnexp_make_fun <- function(){
readMSData(ms_files, mode="inMemory", msLevel. = 1)
}
ondisk_make_fun <- function(){
readMSData(ms_files, mode="onDisk", msLevel. = 1)
}
spectra_make_fun <- function(){
Spectra(ms_files, source=MsBackendMzR()) %>% filterMsLevel(1)
}
rams_make_fun <- function(){
grabMSdata(ms_files, grab_what="MS1")
}
tmzml_make_fun <- function(){
tmzml_names <- paste0(dirname(dirname(ms_files)), "/tmzMLs/",
gsub("mzXML", "tmzML", basename(ms_files)))
dir.create("vignettes/figures/ssc_vignette_renders/tmzMLs")
mapply(tmzmlMaker, ms_files, tmzml_names)
unlink("vignettes/figures/ssc_vignette_renders/tmzMLs", recursive = TRUE)
}
arrow_make_fun <- function(){
msdata <- grabMSdata(ms_files, grab_what="MS1")
write_dataset(msdata$MS1, path = "vignettes/figures/ssc_vignette_renders/pqds")
unlink("vignettes/figures/ssc_vignette_renders/pqds", recursive = TRUE)
}
sql_make_fun <- function(){
msdata <- grabMSdata(ms_files, grab_what="MS1")
MSdb <- dbConnect(RSQLite::SQLite(), "vignettes/figures/ssc_vignette_renders/MSdata.sqlite")
dbWriteTable(MSdb, "MS1", msdata$MS1, overwrite=TRUE)
dbDisconnect(MSdb)
unlink("vignettes/figures/ssc_vignette_renders/MSdata.sqlite")
}
sqlidx_make_fun <- function(){
msdata <- grabMSdata(ms_files, grab_what="MS1")
MSdb_idx <- dbConnect(RSQLite::SQLite(), "vignettes/figures/ssc_vignette_renders/MSdata_idx.sqlite")
dbWriteTable(MSdb_idx, "MS1", msdata$MS1, overwrite=TRUE)
rs <- dbSendQuery(MSdb_idx, "CREATE INDEX mz ON MS1 (mz)")
dbClearResult(rs)
dbDisconnect(MSdb_idx)
unlink("vignettes/figures/ssc_vignette_renders/MSdata_idx.sqlite")
}
make_timings <- microbenchmark(
msnexp_make_fun(), ondisk_make_fun(), spectra_make_fun(), rams_make_fun(),
tmzml_make_fun(), arrow_make_fun(), sql_make_fun(), sqlidx_make_fun(),
times = 10)
saveRDS(make_timings, "vignettes/figures/ssc_vignette_renders/make_timings.rds")Code
make_timings <- readRDS("vignettes/figures/ssc_vignette_renders/make_timings.rds")
make_timings %>%
as.data.frame() %>%
arrange(expr) %>%
mutate(expr=str_remove(expr, "_make_fun\\(\\)")) %>%
mutate(rep_type=case_when(
expr%in%c("msnexp", "ondisk", "spectra", "rams")~"Every R session",
TRUE~"Single-time only"
)) %>%
mutate(expr=factor(expr, levels=c("msnexp", "ondisk", "spectra", "rams",
"tmzml", "arrow", "sql", "sqlidx"),
labels=c("MSnExp", "OnDiskMSnExp", "Spectra", "RaMS",
"tmzMLs", "Arrow", "SQL", "SQL (indexed)"))) %>%
ggplot() +
geom_boxplot(aes(x=expr, y=time/1e9)) +
geom_hline(yintercept = 0) +
facet_wrap(~rep_type, nrow = 1, scales="free_x") +
labs(y="Seconds", x=NULL) +
theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5))
ggsave("vignettes/figures/ssc_vignette_renders/make_time_gp.png", width = 6.5, height = 4, units = "in", device = "png", dpi = 144)
Nothing super surprising here. It should be noted that I ran all of
these without using parallel processing to ensure a fair comparison, but
RaMS and the tmzMLs are almost perfectly parallelizable and can be
read/written in series for a huge time saving. Testing this with 5 cores
across the 10 files resulted in a 3-fold speed improvement (from ~250
seconds to ~85, placing it about on par with the SQL database) for
tmzMLs but not for default RaMS. This makes sense if I/O is the slow
step but reading can be done separately from writing, thus preventing
multiple simultaneous reads but allowing one tmzML to be written while
another is read. Also worth noting is that because
tmzmlMaker calls grabMSdata internally, the
RaMS values are a lower bound on the maximum tmzML speed. The MSnExp,
OnDiskMSnExp, and Spectra packages focus instead on parallelizing their
query step and I can’t seem to run those in parallel. The SQL and Arrow
options similarly write to a single file that’s presumably going to be
“locked” during writing and don’t have parallelization options available
either. Also worth noting is that performing this test on 10 files (~1GB
MS data) doesn’t really challenge memory requirements, so this could be
very different if a larger number of files was used.
Code
unpar_rams <- function(){
print("Unpar RaMS")
lapply(ms_files, grabMSdata)
}
unpar_tmzml <- function(){
print("Unpar tmzML")
tmzml_names <- paste0(dirname(dirname(ms_files)), "/tmzMLs/",
gsub("mzXML", "tmzML", basename(ms_files)))
dir.create("vignettes/figures/ssc_vignette_renders/tmzMLs")
pbapply::pbmapply(tmzmlMaker, ms_files, tmzml_names)
unlink("vignettes/figures/ssc_vignette_renders/tmzMLs", recursive = TRUE)
}
library(BiocParallel)
par_param <- SnowParam(workers = 5, progressbar = TRUE)
par_rams <- function(){
print("Par RaMS")
bplapply(ms_files, grabMSdata, BPPARAM = par_param)
}
par_tmzml <- function(){
print("Par tmzML")
tmzml_names <- paste0(dirname(dirname(ms_files)), "/tmzMLs/",
gsub("mzXML", "tmzML", basename(ms_files)))
dir.create("vignettes/figures/ssc_vignette_renders/tmzMLs")
bpmapply(tmzmlMaker, ms_files, tmzml_names, BPPARAM = par_param)
unlink("vignettes/figures/ssc_vignette_renders/tmzMLs", recursive = TRUE)
}
par_timings <- microbenchmark(par_rams(), unpar_rams(), par_tmzml(), unpar_tmzml(), times = 5)
saveRDS(par_timings, "vignettes/figures/ssc_vignette_renders/par_timings.rds")Code
par_timings <- readRDS("vignettes/figures/ssc_vignette_renders/par_timings.rds")
par_timings %>%
as.data.frame() %>%
separate(expr, into = c("sub_type", "par_type"), sep = "_") %>%
mutate(par_type=str_remove(par_type, "\\(\\)")) %>%
mutate(sub_type=factor(sub_type, levels=c("unpar", "par"),
labels=c("Sequential", "Parallel"))) %>%
mutate(par_type=factor(par_type, levels=c("rams", "tmzml"),
labels=c("RaMS", "tmzMLs"))) %>%
ggplot() +
geom_boxplot(aes(x=par_type, y=time/1e9)) +
geom_hline(yintercept = 0) +
facet_wrap(~sub_type, nrow = 1) +
labs(y="Seconds", x=NULL) +
theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5))
ggsave("vignettes/figures/ssc_vignette_renders/par_time_gp.png", width = 6.5, height = 3, units = "in", device = "png", dpi = 144)
