Using CluMSID with a Publicly Available MetaboLights Data Set
Introduction
In this tutorial, we would like to demonstate the use of CluMSID with a publicly available LC-MS/MS data set deposited on MetaboLights. We chose data set MTBLS433 that can be accessed on the MetaboLights web page (https://www.ebi.ac.uk/metabolights/MTBLS433) and which has been published in the following article:
Kalogiouri, N. P., Alygizakis, N. A., Aalizadeh, R., & Thomaidis, N. S. (2016). Olive oil authenticity studies by target and nontarget LC–QTOF-MS combined with advanced chemometric techniques. Analytical and bioanalytical chemistry, 408(28), 7955-7970.
The authors analysed olive oil of various providence using reversed-phase ultra high performance liquid chromatography–electrospray ionisation quadrupole time of flight tandem mass spectrometry in negative mode with auto-MS/MS fragmentation.
As a representative pooled sample is not provided, we will combine MS2 data from several runs and use the peak picking done by the authors of the study for the merging of MS2 spectra. Some metabolite annotations are also included in the MTBLS433 data set which we will integrate into our analysis.
Extract MS2 spectra from multiple *.mzML files
For demonstration, not all files from the analysis will be included into the analysis. Four data files from the data set have been chosen that represent olive oil samples from different regions in Greece:
YH1_GA7_01_10463.mzML: YH1, from KomiAX1_GB5_01_10470.mzML: AX1, from MegaloxoriLP1_GB3_01_10467.mzML: LP1, from MoriaBR1_GB6_01_10471.mzML: BR1, from Agia Paraskevi
Note that these are mzML files that can be processed the exact same way as mzXML files.
Furthermore, we would like to use the peak picking and annotation
data from the original authors which we can read from the file
m_mtbls433_metabolite_profiling_mass_spectrometry_v2_maf.tsv.
First, we extract MS2 spectra from the respective files
separately by using extractMS2spectra(). Then, we just
combine the resulting lists into one list using base R
functionality:
YH1 <- system.file("extdata", "YH1_GA7_01_10463.mzML",
package = "CluMSIDdata")
AX1 <- system.file("extdata", "AX1_GB5_01_10470.mzML",
package = "CluMSIDdata")
LP1 <- system.file("extdata", "LP1_GB3_01_10467.mzML",
package = "CluMSIDdata")
BR1 <- system.file("extdata", "BR1_GB6_01_10471.mzML",
package = "CluMSIDdata")
YH1list <- extractMS2spectra(YH1)
AX1list <- extractMS2spectra(AX1)
LP1list <- extractMS2spectra(LP1)
BR1list <- extractMS2spectra(BR1)
raw_oillist <- c(YH1list, AX1list, LP1list, BR1list)Merge spectra with external peak list
First, we import the peak list by reading the respective table and
filtering for the relevant information. We only need the columns
metabolite_identification, mass_to_charge and
rentention_time and we would like to replace
"unknown" with an empty field in the
metabolite_identification column. Plus, the features do not
have a unique identifier in the table but we can easily generate that
from m/z and RT. Note that the retention time in the raw data
is given in seconds and in the data table it is in minutes, so we have
to convert. For the sake of consistency, we also change the column
names. We use tidyverse syntax but users can do as they
prefer.
raw_mtbls_df <- system.file("extdata",
"m_mtbls433_metabolite_profiling_mass_spectrometry_v2_maf.tsv",
package = "CluMSIDdata")
require(magrittr)
mtbls_df <- readr::read_delim(raw_mtbls_df, "\t") %>%
dplyr::mutate(metabolite_identification =
stringr::str_replace(metabolite_identification,
"unknown", "")) %>%
dplyr::mutate(id = paste0("M", mass_to_charge, "T", retention_time)) %>%
dplyr::mutate(retention_time = retention_time * 60) %>%
dplyr::select(id,
mass_to_charge,
retention_time,
metabolite_identification) %>%
dplyr::rename(mz = mass_to_charge,
rt = retention_time,
annotation = metabolite_identification)This peak list, or its first three columns, can now be used to merge spectra. We exclude spectra that do not match to any of the peaks in the peak list. As we are not very familiar with instrumental setup, we set the limits for retention time and m/z deviation a little wider. To make an educated guess on mass accuracy, we take a look at an identified metabolite, its measured m/z and its theoretical m/z. We use arachidic acid [M-H]-, whose theoretical m/z is 311.2956:
## Define theoretical m/z
th <- 311.2956
## Get measured m/z for arachidic acid data from mtbls_df
ac <- mtbls_df %>%
dplyr::filter(annotation == "Arachidic acid") %>%
dplyr::select(mz) %>%
as.numeric()
## Calculate relative m/z difference in ppm
abs(th - ac)/th * 1e6
#> [1] 28.91143So, we will work with an an m/z tolerance of 30ppm (which seems rather high for a high resolution mass spectrometer).
A 30ppm mass accuracy necessitates an m/z tolerance of 60ppm, because deviations can go both ways:
Add annotations
To add annotations, we use mtbls_df as well, as
described in the General Tutorial:
Generate distance matrix
For the generation of the distance matrix, too, we use an m/z tolerance of 30ppm:
Explore data
To explore the data, we have a look at a cluster dendrogram:
Figure 1: Circularised dendrogram as a result of agglomerative hierarchical clustering with average linkage as agglomeration criterion based on MS2 spectra similarities of the MTBLS433 LC-MS/MS example data set. Each leaf represents one feature and colours encode cluster affiliation of the features. Leaf labels display feature IDs. Distance from the central point is indicative of the height of the dendrogram.
Since it was not in the focus of their study, the authors identified only a few metabolites. If we look at the positions of these metabolites in the cluster dendrogram, we see that the poly-unsaturated fatty acids alpha-linolenic acid and alpha-linolenic acid are nicely separated from the saturated fatty acids stearic acid and arachidic acid. We would expect the latter to cluster together but a look at the spectra reveals that stearic acid barely produces any fragment ions and mainly contains the unfragmented [M-H]- parent ion:
Figure 2: Barplot for the feature M283.2642T14.62, identified as stearic acid, displaying fragment m/z on the x-axis and intensity normalised to the maximum intensity on the y-axis.
In contrast, arachidic acid produces a much richer spectrum:
Figure 3: Barplot for the feature M311.3046T15.1, identified as arachidic acid, displaying fragment m/z on the x-axis and intensity normalised to the maximum intensity on the y-axis.
Inspecting the features that cluster close to arachidic acid shows that many of them have an exact m/z that conforms with other fatty acids of different chain length or saturation (within the m/z tolerance), e.g. the neighbouring feature M339.2125T15.32 that could be arachidonic acid [M+Cl]-.
Looking at oleic acid [M-H]-, we see that it clusters very closely to M563.5254T13.93, whose m/z is consistent with oleic acid [2M-H]- and some other possible adducts.
As a last example, the only identified metabolite that does not belong to the class of fatty acids is acetosyringone, a phenolic secondary plant metabolite. It forms part of a rather dense cluster in the dendrogram, suggesting high spectral similarities to the other members of the cluster. It would be interesting to try to annotate more of these metabolite to find out if they are also phenolic compounds.
In conclusion, we demonstrated how to use CluMSID with a
publicly available data set from the MetaboLights repository and how to
include external information such as peak lists or feature annotations
into a CluMSID workflow. In doing so, we had a look on a
few example findings that could help to annotate more of the features in
the data set and thereby showed the usefulness of CluMSID
for samples very different from the ones in the other tutorials.
Session Info
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 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: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] magrittr_2.0.3 metaMSdata_1.41.0 metaMS_1.42.0
#> [4] CAMERA_1.63.0 xcms_4.4.0 BiocParallel_1.41.0
#> [7] Biobase_2.67.0 BiocGenerics_0.53.0 CluMSIDdata_1.21.0
#> [10] CluMSID_1.23.0 rmarkdown_2.28
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 rstudioapi_0.17.1
#> [3] sys_3.4.3 jsonlite_1.8.9
#> [5] MultiAssayExperiment_1.33.0 farver_2.1.2
#> [7] MALDIquant_1.22.3 fs_1.6.5
#> [9] zlibbioc_1.52.0 vctrs_0.6.5
#> [11] base64enc_0.1-3 progress_1.2.3
#> [13] htmltools_0.5.8.1 S4Arrays_1.6.0
#> [15] Formula_1.2-5 SparseArray_1.6.0
#> [17] mzID_1.45.0 sass_0.4.9
#> [19] KernSmooth_2.23-24 bslib_0.8.0
#> [21] htmlwidgets_1.6.4 plyr_1.8.9
#> [23] impute_1.81.0 plotly_4.10.4
#> [25] cachem_1.1.0 buildtools_1.0.0
#> [27] igraph_2.1.1 lifecycle_1.0.4
#> [29] iterators_1.0.14 pkgconfig_2.0.3
#> [31] Matrix_1.7-1 R6_2.5.1
#> [33] fastmap_1.2.0 GenomeInfoDbData_1.2.13
#> [35] MatrixGenerics_1.19.0 clue_0.3-65
#> [37] digest_0.6.37 pcaMethods_1.99.0
#> [39] colorspace_2.1-1 GGally_2.2.1
#> [41] S4Vectors_0.44.0 Hmisc_5.2-0
#> [43] GenomicRanges_1.59.0 Spectra_1.16.0
#> [45] fansi_1.0.6 httr_1.4.7
#> [47] abind_1.4-8 compiler_4.4.1
#> [49] bit64_4.5.2 withr_3.0.2
#> [51] doParallel_1.0.17 backports_1.5.0
#> [53] htmlTable_2.4.3 DBI_1.2.3
#> [55] ggstats_0.7.0 highr_0.11
#> [57] gplots_3.2.0 MASS_7.3-61
#> [59] MsExperiment_1.9.0 DelayedArray_0.33.1
#> [61] gtools_3.9.5 caTools_1.18.3
#> [63] mzR_2.41.0 tools_4.4.1
#> [65] foreign_0.8-87 PSMatch_1.10.0
#> [67] ape_5.8 nnet_7.3-19
#> [69] glue_1.8.0 dbscan_1.2-0
#> [71] nlme_3.1-166 QFeatures_1.16.0
#> [73] grid_4.4.1 checkmate_2.3.2
#> [75] cluster_2.1.6 reshape2_1.4.4
#> [77] generics_0.1.3 gtable_0.3.6
#> [79] tzdb_0.4.0 preprocessCore_1.68.0
#> [81] tidyr_1.3.1 hms_1.1.3
#> [83] sna_2.8 data.table_1.16.2
#> [85] MetaboCoreUtils_1.15.0 utf8_1.2.4
#> [87] XVector_0.46.0 foreach_1.5.2
#> [89] pillar_1.9.0 stringr_1.5.1
#> [91] vroom_1.6.5 limma_3.63.0
#> [93] robustbase_0.99-4-1 dplyr_1.1.4
#> [95] lattice_0.22-6 bit_4.5.0
#> [97] RBGL_1.82.0 tidyselect_1.2.1
#> [99] maketools_1.3.1 knitr_1.48
#> [101] gridExtra_2.3 IRanges_2.41.0
#> [103] ProtGenerics_1.38.0 SummarizedExperiment_1.36.0
#> [105] stats4_4.4.1 xfun_0.48
#> [107] statmod_1.5.0 MSnbase_2.33.0
#> [109] matrixStats_1.4.1 DEoptimR_1.1-3
#> [111] stringi_1.8.4 UCSC.utils_1.2.0
#> [113] statnet.common_4.10.0 lazyeval_0.2.2
#> [115] yaml_2.3.10 evaluate_1.0.1
#> [117] codetools_0.2-20 MsCoreUtils_1.19.0
#> [119] tibble_3.2.1 graph_1.85.0
#> [121] BiocManager_1.30.25 cli_3.6.3
#> [123] affyio_1.77.0 rpart_4.1.23
#> [125] munsell_0.5.1 jquerylib_0.1.4
#> [127] network_1.18.2 Rcpp_1.0.13
#> [129] GenomeInfoDb_1.43.0 MassSpecWavelet_1.73.0
#> [131] coda_0.19-4.1 XML_3.99-0.17
#> [133] parallel_4.4.1 readr_2.1.5
#> [135] ggplot2_3.5.1 prettyunits_1.2.0
#> [137] AnnotationFilter_1.31.0 bitops_1.0-9
#> [139] viridisLite_0.4.2 MsFeatures_1.15.0
#> [141] scales_1.3.0 affy_1.85.0
#> [143] ncdf4_1.23 purrr_1.0.2
#> [145] crayon_1.5.3 rlang_1.1.4
#> [147] vsn_3.74.0