Note: This vignette is built without network access. Online databases are
read from stored copies included with MetMashR, and REST API lookups use
stored response caches with cache_mode = "offline", so a value that is not in
a cache is left as NA instead of being queried. To query the online services
yourself, use your own cache and the default cache_mode = "update". The stored
copies were made with inst/scripts/make_vignette_snapshots.R.
# Getting Started
The latest versions of struct and MetMashR that are compatible
with your
current R version can be installed using BiocManager.
# install BiocManager if not present
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
# install MetMashR and dependencies
BiocManager::install("MetMashR")
Once installed you can activate the packages in the usual way:
# load the packages
library(MetMashR)
library(ggplot2)
library(structToolbox)
library(dplyr)
library(DT)
library(ggplotify)
MTox700+ is a list of toxicologically relevant metabolites derived from publications, public databases and relevant toxicological assays.
In this vignette you will learn how to combine MTox700+ with the PathBank pathway database to see which chemical classes MTox700+ contains and how many of its metabolites have pathway information in PathBank.
The MTox700+ database can be imported using the MTox700plus_database object,
which downloads a release of the database from GitHub. In this vignette we
read a stored copy of the database with rds_database(). A table describing
the stored copies included with MetMashR records which release was stored
and when it was downloaded.
# where and when the stored copies were retrieved
snapshots <- readRDS(system.file(
"extdata/vignette_snapshots/snapshot_info.rds",
package = "MetMashR"
))
snapshots[
snapshots$resource %in% c("MTox700+", "PathBank"),
c("resource", "version", "retrieved")
]
#> resource version retrieved
#> MTox700plus.rds MTox700+ v1.0 2024-01-25
#> PathBank_primary_metabolites.rds PathBank primary 2023-11-23
# prep object (stored copy of MTox700+)
MT <- rds_database(
source = system.file(
"extdata/vignette_snapshots/MTox700plus.rds",
package = "MetMashR"
),
tag = "MTox700+",
.writable = FALSE
)
# import
df <- read_database(MT)
# show contents
.DT(df)
This example uses MTox700+ release v1.0, downloaded on 2024-01-25.
# prepare workflow that uses MTox700+ as a source
M <-
import_source() +
trim_whitespace(
column_name = ".all",
which = "both",
whitespace = "[\\h\\v]"
)
# apply
M <- model_apply(M, MT)
The chemical (or “metabolite”) space covered by the MTox700+ database can be explored in several ways using the data included in the database, which contains information about the structural classification of the metabolites based on ChemOnt (a chemical taxonomy) and ClassyFire (software to compute the taxonomy of a structure) [10.1186/s13321-016-0174-y].
In this plot we show the number of metabolites in the MTox700+ database that are assigned to a “superclass” of molecules.
# initialise chart object
C <- annotation_bar_chart(
factor_name = "superclass",
label_rotation = TRUE,
label_location = "outside",
label_type = "percent",
legend = TRUE
)
# plot
g <-
chart_plot(C, predicted(M)) +
ylim(c(0, 600)) +
guides(
fill = guide_legend(ncol = 1, title = NULL),byrow=TRUE) +
theme(
legend.position = "right",
legend.margin = margin(),
legend.key.size = unit(0.5, "cm"),
legend.text = element_text(size = 10)
)
# layout
leg <- cowplot::get_legend(g)
cowplot::plot_grid(g + theme(legend.position = "none"), leg,
nrow = 1,
rel_widths = c(50, 50)
)
To explore the biological space covered by the metabolites in MTox700+ we need mash the database with additional information about the biological pathways that the metabolites are part of.
We use PathBank for this purpose. MetMashR
already includes a database source for PathBank.
MetMashR provides the PathBank_metabolite_database() object to import the
PathBank database. You can choose to import:
The “complete” database is a >50mb download, and unzipped is >1Gb. Unzipping and caching of the database is handled by [BiocFileCache].
For the vignette we use a stored copy of the “primary” PathBank database, keeping only the columns used below.
We can use the database in two ways:
To explore the biological space covered by MetMashR we will do both.
It is useful to visualise the overlap between PathBank and MTox700+. MTox700+ is a much smaller database due to it being a curated list of metabolites with toxicologial relevance, and PathBank is more general.
In th example below we import PathBank as a source, and use a venn diagram to compare the overlap between inchikey identifiers in PathBank and MTox700+.
# object M already contains the MTox700+ database as a source
# prepare PathBank as a source (stored copy of the primary database)
P <- rds_database(
source = system.file(
"extdata/vignette_snapshots/PathBank_primary_metabolites.rds",
package = "MetMashR"
),
tag = "PathBank",
.writable = FALSE
)
# import
P <- read_source(P)
# prepare chart
C <- annotation_venn_chart(
factor_name = c("inchikey", "InChI.Key"), legend = FALSE,
fill_colour = ".group",
line_colour = "white"
)
# plot
g1 <- chart_plot(C, predicted(M), P)
C <- annotation_upset_chart(
factor_name = c("inchikey", "InChI.Key"),
)
g2 <- as.ggplot(chart_plot(C, predicted(M), P)) +
theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))
cowplot::plot_grid(g1, g2, nrow = 1, labels = c("Venn diagram", "UpSet plot"))
PathBank does not report a version number, so we report the date the stored copy was downloaded: 2023-11-23.
These charts compare full InChIKeys in MTox700+ with those in the selected PathBank primary-pathway dataset. To report the overlap, we count the distinct, non-missing keys shared by both resources and divide by the number of distinct, non-missing MTox700+ keys.
# distinct, non-missing InChIKeys
keys <- function(x) {
x <- unique(x)
x[!is.na(x) & x != ""]
}
mtox_keys <- keys(predicted(M)$data$inchikey)
pathbank_keys <- keys(P$data$InChI.Key)
shared_keys <- intersect(mtox_keys, pathbank_keys)
# MTox700+ records with no InChIKey
n_missing <- sum(is.na(predicted(M)$data$inchikey) |
predicted(M)$data$inchikey == "")
data.frame(
mtox_keys = length(mtox_keys),
shared_keys = length(shared_keys),
percent_shared = round(100 * length(shared_keys) / length(mtox_keys), 1),
mtox_records_without_key = n_missing
)
#> mtox_keys shared_keys percent_shared mtox_records_without_key
#> 1 722 379 52.5 0
Of the 722 distinct MTox700+ InChIKeys, 379 (52.5%) are also in the PathBank primary-pathway dataset; 0 MTox700+ records have no InChIKey. An unmatched key means that it was not found in this PathBank subset; it does not show that the metabolite has no known pathway associations.
To combine the pathway information in PathBank with the MTox700+ database we
can use PathBank as a lookup table based on inchikeys. To do this we use the
database_lookup object.
We request a number of columns from PathBank, including pathway information and additional identifiers such as HMDB ID and KEGG ID.
# prepare object
X <- database_lookup(
query_column = "inchikey",
database = P$data,
database_column = "InChI.Key",
include = c(
"PathBank.ID", "Pathway.Name", "Pathway.Subject", "Species",
"HMDB.ID", "KEGG.ID", "ChEBI.ID", "DrugBank.ID", "SMILES"
),
suffix = ""
)
# apply
X <- model_apply(X, predicted(M))
We can now visualise e.g. the subject of the pathways captured by the MTox700+ database. This chart counts records after the PathBank join. One metabolite can contribute more than once if it matches several pathways.
C <- annotation_bar_chart(
factor_name = "Pathway.Subject",
label_rotation = TRUE,
label_location = "outside",
label_type = "percent",
legend = TRUE
)
chart_plot(C, predicted(X)) + ylim(c(0, 17500))
The percentages describe these joined records, including any unmatched records retained by the lookup. This example includes all returned species, so the chart alone does not establish coverage of human toxicology.
Metabolites can appear in multiple pathways. The PathBank database therefore contains multiple records for the same metabolite, and the relationship between MTox700+ and PathBank is one-to-many.
After obtaining pathway information from PathBank the new table has many more rows than the original MTox700+ database, as each MTox700+ record has been replicated for each match in the PathBank database.
e.g. after importing MTox700+ the number of records was:
# Number in MTox700+
nrow(predicted(M)$data)
#> [1] 722
After combing with PathBank the number of records is:
# Number after PathBank lookup
nrow(predicted(X)$data)
#> [1] 25092
Sometimes it is useful to collapse this information into a single record per
metabolite. We can use the combine_records object and its helper functions to
do this in a MetMashR workflow. Before collapsing the records we save the
full table of matches, which has one row per metabolite and pathway.
# save the full table of matches before collapsing
pathway_matches <- predicted(X)$data
# prepare object
X <- database_lookup(
query_column = "inchikey",
database = P$data,
database_column = "InChI.Key",
include = c(
"PathBank.ID", "Pathway.Name", "Pathway.Subject", "Species",
"HMDB.ID", "KEGG.ID", "ChEBI.ID", "DrugBank.ID", "SMILES"
),
suffix = ""
) +
combine_records(
group_by = "inchikey",
default_fcn = fuse_unique(" || ")
)
# apply
X <- model_apply(X, predicted(M))
We have used the fuse_unique() helper function so that records for each inchikey are
combined into a single record by only retaining unique values in each field
(column). If there are multiple unique values for a field then they are
combined into a single string using the ” || ” separator.
We can now extract the pathways associated with a particular metabolite. For example Glycolic acid:
# get index of metabolite
w <- which(predicted(X)$data$metabolite_name == "Glycolic acid")
The pathways associated with Glycolic acid are:
# print list of pathways
predicted(X)$data$Pathway.Name[w]
#> [1] "Inner Membrane Transport || Glycolate and Glyoxylate Degradation || D-Arabinose Degradation I || Ethylene Glycol Degradation"
The summary lists unique values from each column separately. Values in the same position in those lists may not come from the same original row. Use the saved full table to see which pathway ID, pathway name and species belong together:
glycolic <- pathway_matches[
pathway_matches$metabolite_name == "Glycolic acid",
c("PathBank.ID", "Pathway.Name", "Species")
]
.DT(glycolic)
Combining MTox700+ with PathBank shows which MTox700+ metabolites have pathway information in the PathBank primary-pathway dataset, and which pathway subjects those matches belong to. The comparison has limits. Matching uses full InChIKeys, so a metabolite recorded in a different form (e.g. a different charge state or stereoisomer) is not matched. Only the primary-pathway subset of PathBank is used, and all species are included. The pathway charts count joined records rather than metabolites, so metabolites in many pathways have more weight.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplotify_0.1.4 DT_0.34.0
#> [3] dplyr_1.2.1 structToolbox_1.25.0
#> [5] ggplot2_4.0.3 metabolomicsWorkbenchR_1.23.0
#> [7] MetMashR_1.7.3 struct_1.25.0
#> [9] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 farver_2.1.2
#> [3] S7_0.2.2 fastmap_1.2.0
#> [5] digest_0.6.39 lifecycle_1.0.5
#> [7] magrittr_2.0.5 compiler_4.6.1
#> [9] rlang_1.3.0 sass_0.4.10
#> [11] tools_4.6.1 yaml_2.3.12
#> [13] data.table_1.18.6.1 knitr_1.52
#> [15] labeling_0.4.3 S4Arrays_1.13.2
#> [17] htmlwidgets_1.6.4 sp_2.2-4
#> [19] DelayedArray_0.39.8 plyr_1.8.9
#> [21] xml2_1.6.0 RColorBrewer_1.1-3
#> [23] aplot_0.3.2 abind_1.4-8
#> [25] withr_3.0.3 purrr_1.2.2
#> [27] BiocGenerics_0.59.12 grid_4.6.1
#> [29] stats4_4.6.1 scales_1.4.0
#> [31] tinytex_0.61 MultiAssayExperiment_1.39.1
#> [33] dichromat_2.0-1 SummarizedExperiment_1.43.0
#> [35] cli_3.6.6 rmarkdown_2.32
#> [37] generics_0.1.4 otel_0.3.0
#> [39] httr_1.4.9 cachem_1.1.0
#> [41] ggthemes_7.0.0 BiocManager_1.30.27
#> [43] XVector_0.53.0 matrixStats_1.5.0
#> [45] vctrs_0.7.3 yulab.utils_0.2.5
#> [47] Matrix_1.7-6 jsonlite_2.0.0
#> [49] bookdown_0.48 patchwork_1.3.2
#> [51] gridGraphics_0.5-1 IRanges_2.47.5
#> [53] S4Vectors_0.51.10 magick_2.9.1
#> [55] crosstalk_1.2.2 jquerylib_0.1.4
#> [57] tidyr_1.3.2 ggVennDiagram_1.5.7
#> [59] glue_1.8.1 cowplot_1.2.0
#> [61] stringi_1.8.9 gtable_0.3.6
#> [63] GenomicRanges_1.65.4 tibble_3.3.1
#> [65] pillar_1.11.1 rappdirs_0.3.4
#> [67] htmltools_0.5.9 Seqinfo_1.3.2
#> [69] R6_2.6.1 httr2_1.3.0
#> [71] evaluate_1.0.5 lattice_0.23-1
#> [73] Biobase_2.73.2 openxlsx_4.2.9
#> [75] ggfun_0.2.2 bslib_0.12.0
#> [77] Rcpp_1.1.2 zip_3.0.2
#> [79] gridExtra_2.3.1 SparseArray_1.13.4
#> [81] xfun_0.62 forcats_1.0.1
#> [83] MatrixGenerics_1.25.0 fs_2.1.0
#> [85] pkgconfig_2.0.3