Adafede
commented
1 year ago As a side note, here is some benchmark I did regarding the different backends:
#################################### Paths ####################################
file <- "CompDb.Hsapiens.HMDB.4.0"
file_db <- file |>
paste0(".sqlite")
file_mgf <- file |>
paste0(".mgf")
################ Taking small HMDB from CompoundDb as example #################
hmdb_file <-
system.file("sdf/HMDB_sub.sdf.gz", package = "CompoundDb")
cmps <- CompoundDb::compound_tbl_sdf(hmdb_file)
xml_path <- system.file("xml", package = "CompoundDb")
spctra <- CompoundDb::msms_spectra_hmdb(xml_path)
metad <-
CompoundDb::make_metadata(
source = "HMDB",
url = "http://www.hmdb.ca",
source_version = "4.0",
source_date = "2017-09",
organism = "Hsapiens"
)
################################## Functions ##################################
export_spectra <- function(file,
spectra,
dir = ".",
cmps = NULL,
metad = NULL) {
switch(
EXPR = gsub(
pattern = ".*\\.",
replacement = "",
x = file
),
"mgf" = {
Spectra::Spectra(object = spectra) |>
Spectra::export(
backend = MsBackendMgf::MsBackendMgf(),
file = file
)
},
"sqlite" = {
if (is.null(cmps)) {
cmps <-
data.frame(
compound_id = if (!is.null(spectra$compound_id)) {
spectra$compound_id
} else {
no <- rep(
NA,
nrow(spectra)
)
},
name = if (!is.null(spectra$name)) {
spectra$name
} else {
no <- rep(
NA,
nrow(spectra)
)
},
inchi = if (!is.null(spectra$inchi)) {
spectra$inchi
} else {
no <- rep(
NA,
nrow(spectra)
)
},
inchikey = if (!is.null(spectra$inchikey)) {
spectra$inchikey
} else {
no <- rep(
NA,
nrow(spectra)
)
},
exactmass = if (!is.null(spectra$exactmass)) {
spectra$exactmass
} else {
no <- rep(
NA_integer_,
nrow(spectra)
)
},
formula = if (!is.null(spectra$formula)) {
spectra$formula
} else {
no <- rep(
NA,
nrow(spectra)
)
},
synonyms = if (!is.null(spectra$synonyms)) {
spectra$synonyms
} else {
no <- rep(
NA,
nrow(spectra)
)
}
)
}
if (is.null(metad)) {
metad <-
data.frame(
name = c(
"source",
"url",
"source_version",
"source_date",
"organism"
),
value = c(
NA,
NA,
0,
Sys.time() |>
as.character(),
NA
)
)
}
CompoundDb::createCompDb(
x = cmps,
metadata = metad,
msms_spectra = spectra,
path = dir,
dbFile = file
)
}
)
}
extract_spectra <- function(object) {
## issues
incoherent_colnames <- c(
ms_level = "msLevel",
precursor_intensity = "precursorIntensity"
)
incoherent_logical <- c("predicted")
incoherent_integer <- c("spectrum_id")
##
peaks <- object |>
Spectra::peaksData() |>
data.frame() |>
dplyr::group_by(group) |>
dplyr::summarize_all(list)
spectra <- object |>
Spectra::spectraData() |>
data.frame()
spectra$mz <- peaks$mz
spectra$intensity <- peaks$intensity
## synonyms issue
spectra <- spectra |>
dplyr::group_by(dplyr::across(c(-dplyr::any_of("synonym")))) |>
dplyr::summarize(dplyr::across(.cols = where(is.list), .fns = as.character)) |>
dplyr::ungroup()
## columns types issue
spectra <- spectra |>
dplyr::mutate(dplyr::across(
.cols = dplyr::any_of(incoherent_logical),
.fns = as.logical
)) |>
dplyr::mutate(dplyr::across(
.cols = dplyr::any_of(incoherent_integer),
.fns = as.integer
))
## columns names issue
spectra <- spectra |>
dplyr::select(-c(dplyr::any_of(incoherent_colnames))) |>
dplyr::rename(dplyr::any_of(incoherent_colnames))
return(spectra)
}
import_spectra <- function(file) {
switch(
EXPR = gsub(
pattern = ".*\\.",
replacement = "",
x = file
),
"mgf" = {
file |>
MsBackendMgf::readMgf() |>
Spectra::Spectra()
},
"sqlite" = {
file |>
CompoundDb::CompDb() |>
CompoundDb::Spectra()
}
)
}
################################### Script ###################################
## Exporting
### To SQLite
export_spectra(
file = file_db,
spectra = spctra,
cmps = cmps,
metad = metad
)
### To MGF
export_spectra(file = file_mgf, spectra = spctra)
## Importing
### From SQLite
sps_from_sqlite <- file_db |>
import_spectra()
### From MGF
sps_from_mgf <- file_mgf |>
import_spectra()
################################# Parameters #################################
params_sim <- MetaboAnnotation::CompareSpectraParam(
ppm = 10,
tolerance = 0.02,
FUN = MsCoreUtils::gnps,
requirePrecursor = TRUE,
THRESHFUN = function(x) {
which(x >= 0.6)
}
)
################################## Benchmark ##################################
mbm_1 <- microbenchmark::microbenchmark(
matches_mgf_mgf = MetaboAnnotation::matchSpectra(
query = sps_from_mgf,
target = sps_from_mgf,
param = params_sim
),
matches_mgf_sqlite = MetaboAnnotation::matchSpectra(
query = sps_from_mgf,
target = sps_from_sqlite,
param = params_sim
),
matches_sqlite_mgf = MetaboAnnotation::matchSpectra(
query = sps_from_sqlite,
target = sps_from_mgf,
param = params_sim
),
matches_sqlite_sqlite = MetaboAnnotation::matchSpectra(
query = sps_from_sqlite,
target = sps_from_sqlite,
param = params_sim
),
times = 200
)
mbm_2 <- microbenchmark::microbenchmark(
matches_mgf_mgf_1 = MetaboAnnotation::matchSpectra(
query = file_mgf |>
import_spectra(),
target = sps_from_mgf,
param = params_sim
),
matches_mgf_mgf_2 = MetaboAnnotation::matchSpectra(
query = sps_from_mgf,
target = file_mgf |>
import_spectra(),
param = params_sim
),
matches_mgf_sqlite_1 = MetaboAnnotation::matchSpectra(
query = file_mgf |>
import_spectra(),
target = sps_from_sqlite,
param = params_sim
),
matches_mgf_sqlite_2 = MetaboAnnotation::matchSpectra(
query = sps_from_mgf,
target = file_db |>
import_spectra(),
param = params_sim
),
matches_sqlite_mgf_1 = MetaboAnnotation::matchSpectra(
query = file_db |>
import_spectra(),
target = sps_from_mgf,
param = params_sim
),
matches_sqlite_mgf_2 = MetaboAnnotation::matchSpectra(
query = sps_from_sqlite,
target = file_mgf |>
import_spectra(),
param = params_sim
),
matches_sqlite_sqlite_1 = MetaboAnnotation::matchSpectra(
query = file_db |>
import_spectra(),
target = sps_from_sqlite,
param = params_sim
),
matches_sqlite_sqlite_2 = MetaboAnnotation::matchSpectra(
query = sps_from_sqlite,
target = file_db |>
import_spectra(),
param = params_sim
),
times = 200
)
mbm_3 <- microbenchmark::microbenchmark(
matches_mgf_mgf = MetaboAnnotation::matchSpectra(
query = file_mgf |>
import_spectra(),
target = file_mgf |>
import_spectra(),
param = params_sim
),
matches_mgf_sqlite = MetaboAnnotation::matchSpectra(
query = file_mgf |>
import_spectra(),
target = file_db |>
import_spectra(),
param = params_sim
),
matches_sqlite_mgf = MetaboAnnotation::matchSpectra(
query = file_db |>
import_spectra(),
target = file_mgf |>
import_spectra(),
param = params_sim
),
matches_sqlite_sqlite = MetaboAnnotation::matchSpectra(
query = file_db |>
import_spectra(),
target = file_db |>
import_spectra(),
param = params_sim
),
times = 200
)
#################################### Plots ####################################
mbm_1 |>
ggplot2::autoplot()
mbm_2 |>
ggplot2::autoplot()
mbm_3 |>
ggplot2::autoplot()
################################### Cleanup ###################################
file.remove(file_mgf)
file.remove(file_db)
################################# Session Info #################################
sessionInfo()
michaelwitting
jorainer
Hi,
I am trying to regularly match 1000+ query spectra to 100,000+ target spectra. Currently, this is easily doable in minutes using
matchmsin python. I never obtained such results usingMetaboAnnotationin R.Here below a tentative script trying to describe what I would like to achieve. It currently uses
MsBackendMgf, but this is not the issue. I had to subset the data to remain in reasonable times (approx. 5 minutes).Happy to further develop and exchange if needed.
CODE: