Analysis on E. ulmoides dataset
10 April 2023
Contents
1 Abstract
Untargeted mass spectrometry is a robust tool for biological research, but researchers universally time consumed by dataset parsing. We developed MCnebula, a novel visualization strategy proposed with multidimensional view, termed multi-chemical nebulae, involving in scope of abundant classes, classification, structures, sub-structural characteristics and fragmentation similarity. Many state-of-the-art technologies and popular methods were incorporated in MCnebula workflow to boost chemical discovery. Notably, MCnebula can be applied to explore classification and structural characteristics of unknown compounds that beyond the limitation of spectral library. MCnebula was integrated in R package and public available for custom R statistical pipeline analysis. Now, MCnebula2 (R object-oriented programming with S4 system) is further available for more friendly applications.
2 Introduction
We know that the analysis of untargeted LC-MS/MS dataset generally begin with feature detection. It detects ‘peaks’ as features in MS1 data. Each feature may represents a compound, and assigned with MS2 spectra. The MS2 spectra was used to find out the compound identity. The difficulty lies in annotating these features to discover their compound identity, mining out meaningful information, so as to serve further biological research. Herein, a classified visualization method, called MCnebula, was used for addressing this difficulty. MCnebula utilizes the state-of-the-art computer prediction technology, SIRIUS workflow (SIRIUS, ZODIAC, CSI:fingerID, CANOPUS)1–5, for compound formula prediction, structure retrieve and classification prediction. MCnebula integrates an abundance-based classes (ABC) selection algorithm into features annotation: depending on the user, MCnebula focuses chemical classes with more or less features in the dataset (the abundance of classes), visualizes them, and displays the features they involved; these classes can be dominant structural classes or sub-structural classes. With MCnebula, we can switch from untargeted to targeted analysis, focusing precisely on the compound or chemical class of interest to the researcher.
Eucommia ulmoides Oliv. (E. ulmoides), as a traditional Chinese medicine (TCM), after being processed with saline water, was applied to the treatment of renal diseases for a long time in China. Due to its complex composition, discovering chemical changes during processing (such as processed with saline water) is challenging. We would next demonstrate the addressing of this challenge with MCnebula, which may be enlightening for the study of phytopharmaceuticals.
3 Set-up
Load the R package used for analysis. In the following analysis
process, to illustrate the source of the function, we use the symbol :: to
mark the functions, e.g., dplyr::filter. The functions that were not marked
may source from MCnebula2 or the packages that R (version 4.2) loaded by
default.
library(MCnebula2)
library(exMCnebula2)4 Integrate data and Create Nebulae
4.1 Initialize analysis
Set SIRIUS project path and its version to initialize mcnebula object.
mcn <- mcnebula()
mcn <- initialize_mcnebula(mcn, "sirius.v4", ".")
ion_mode(mcn) <- "neg"Create a temporary folder to store the output data.
tmp <- paste0(tempdir(), "/temp_data")
dir.create(tmp, F)In order to demonstrate the process of analyzing data with MCnebula2,
we provide a ‘mcnebula’ object that was extracted in advance using the
collate_used function, which means that all the data
used in the subsequent analysis has already stored in this ‘mcnebula’
object, without the need to obtain it from the original Project
directory. This avoids the hassle of downloading and storing a dozen
GB of raw files. The following, we
use the collated dataset containing 1612 features
with chemical formula identification.
exfiles <- system.file("extdata", package = "exMCnebula2")Load the ‘.rdata’ file.
load(paste0(exfiles, "/mcn_herbal1612.rdata"))
mcn <- mcn_herbal1612
rm(mcn_herbal1612)
export_path(mcn) <- tmp4.2 Filter candidates
Suppose we predicted a potential compound represented by LC-MS/MS spectrum, and obtained the candidates of chemical molecular formula, structure and chemical class. These candidates include both positive and negative results: for chemical molecular formula and chemical structure, the positive prediction was unique; for chemical class, multiple positive predictions that belong to various classification were involved. We did not know the exact negative and positive. Normally, we ranked and filtered these according to the scores. There were numerious scores, for isotopes, for mass error, for structural similarity, for chemical classes… Which score selected to rank candidates depends on the purpose of research. Such as:
- To find out the chemical structure mostly be positive, ranking the candidates by structural score.
- To determine whether the potential compound may be of a certain chemical classes, ranking the candidates by the classified score.
Ether by filter_formula(), filter_structure() or filter_ppcp(), the
candidate with top score can be obtained. However, for the three module
(formula, structure, classes), sometimes thier top score candidates were not in
line with each other. That is, thier top score towards different chemical
molecular formulas. To find out the corresponding data in other modules,
create_reference() should be performed to establish the ‘specific_candidate’
for subsequent filtering.
mcn <- filter_structure(mcn)
mcn <- create_reference(mcn)
mcn <- filter_formula(mcn, by_reference = T)4.3 Filter chemical classes
The PPCP (Posterior Probability of Classification Prediction) data for each ‘feature’ contains the prediction of thousands of classes for the potential compound (even if the chemical structure was unknown). See http://www.nature.com/articles/s41587-020-0740-8 for details about the prediction. The data contains attributes of:
class.name: name of classes.pp.value: value of posterior probability.hierarchy: hierarchy of classes in the taxonomy. See https://jcheminf.biomedcentral.com/articles/10.1186/s13321-016-0174-y for details about hierarchy and taxonomy of chemical classification.- …
The method create_stardust_classes() use these inner attributes to filter
classes candidates for each ‘feature’.
Compared to the chemical class filtering within PPCP data by
create_stardust_classes(), the filtering within ‘stardust_classes’ data by
cross_filter_stardust() is fundamentally different.
- For
create_stardust_classes(), the PPCP data belongs to each ‘feature’. When performing the filtering, only simple threshold conditions or absolute conditions are set to filter the chemical classes; there is no crossover between the different attributes and no crossover between the ‘features’. Therefore, we consider this as ‘inner’ filtering. - For
cross_filter_stardust(), the data of the chemical classes and their classified ‘features’, i.e. ‘stardust_classes’ data, were combined and then grouped upon the chemical classes. After grouping, each chemical class has a certain quantity of “features”. When filtering, statistics may be performed on ‘features’ data within a group; statistics may be performed on these data in conjunction with ‘features_annotation’ data; and statistics may be performed to compare groups with each other. As its crossover, we consider this as ‘cross’ filtering.
Use help(cross_filter_stardust) to get more details about the algorithm.
mcn <- create_stardust_classes(mcn)
mcn <- create_features_annotation(mcn)
mcn <- cross_filter_stardust(mcn,
max_ratio = 0.1, cutoff = 0.4,
identical_factor = 0.6
)
classes <- unique(stardust_classes(mcn)$class.name)
table.filtered.classes <- backtrack_stardust(mcn)Manually filter some repetitive classes or sub-structural classes. By means of Regex matching, we obtained a number of recurring name of chemical classes that would contain manay identical compounds as their sub-structure.
classes## [1] "Lactones"
## [2] "Ketones"
## [3] "Benzoic acids and derivatives"
## [4] "Methoxyphenols"
## [5] "Sugar acids and derivatives"
## [6] "Phenylpropanoids and polyketides"
## [7] "Benzoyl derivatives"
## [8] "Fatty acid esters"
## [9] "Unsaturated fatty acids"
## [10] "Hydroxy fatty acids"
## [11] "Hydroxy acids and derivatives"
## [12] "Lineolic acids and derivatives"
## [13] "Dialkyl ethers"
## [14] "Cyclic alcohols and derivatives"
## [15] "Benzoic acid esters"
## [16] "Alpha hydroxy acids and derivatives"
## [17] "Lignans, neolignans and related compounds"
## [18] "Lignan glycosides"
## [19] "Disaccharides"
## [20] "Monoterpenoids"
## [21] "Bicyclic monoterpenoids"
## [22] "Iridoids and derivatives"
## [23] "Tertiary alcohols"
## [24] "Beta hydroxy acids and derivatives"
## [25] "Fatty acyl glycosides"
## [26] "Furofurans"
## [27] "Furofuran lignans"
## [28] "Terpene glycosides"
## [29] "Alkyl glycosides"
## [30] "Methoxybenzoic acids and derivatives"
## [31] "Long-chain fatty acids"
## [32] "Medium-chain fatty acids"
## [33] "Furanoid lignans"
## [34] "Fatty acyl glycosides of mono- and disaccharides"
## [35] "Vinylogous esters"
## [36] "Iridoid O-glycosides"
## [37] "Dimethoxybenzenes"pattern <- c("fatty acid", "hydroxy")
dis <- unlist(lapply(pattern, grep, x = classes, ignore.case = T))
dis <- classes[dis]
dis## [1] "Fatty acid esters" "Unsaturated fatty acids"
## [3] "Hydroxy fatty acids" "Long-chain fatty acids"
## [5] "Medium-chain fatty acids" "Hydroxy fatty acids"
## [7] "Hydroxy acids and derivatives" "Alpha hydroxy acids and derivatives"
## [9] "Beta hydroxy acids and derivatives"Remove these classes.
mcn <- backtrack_stardust(mcn, dis, remove = T)4.4 Create Nebulae
Create Nebula-Index data. This data created based on ‘stardust_classes’ data.
mcn <- create_nebula_index(mcn)Whether it is all filtered by the algorithm provided by MCnebula2’s function or custom filtered for some chemical classes, we now have a data called ‘nebula_index’. This data records a number of chemical classes and the ‘features’ attributed to them. The subsequent analysis process or visualization will be based on it. Each chemical class is considered as a ‘nebula’ and its classified ‘features’ are the components of these ‘nebulae’. In the visualization, these ‘nebulae’ will be visualized as networks. Formally, we call these ‘nebulae’ formed on the basis of ‘nebula_index’ data as Child-Nebulae. In comparison, when we put all the ‘features’ together to form a large network, then this ‘nebula’ is called Parent-Nebulae.
mcn <- compute_spectral_similarity(mcn)
mcn <- create_parent_nebula(mcn)
mcn <- create_child_nebulae(mcn)4.5 Visualize Nebulae
Create layouts for Parent-Nebula or Child-Nebulae visualizations.
mcn <- create_parent_layout(mcn)
mcn <- create_child_layouts(mcn)
mcn <- activate_nebulae(mcn)The available chemical classes for visualization and its sequence in storage.
table.nebulae <- visualize(mcn)## [INFO] MCnebula2: visualize## Specify item as following to visualize:table.nebulae## # A tibble: 29 × 3
## seq hierarchy class.name
## <int> <dbl> <chr>
## 1 1 5 Alkyl glycosides
## 2 2 5 Benzoic acid esters
## 3 3 4 Benzoic acids and derivatives
## 4 4 4 Benzoyl derivatives
## 5 5 5 Bicyclic monoterpenoids
## 6 6 5 Cyclic alcohols and derivatives
## 7 7 5 Dialkyl ethers
## 8 8 5 Dimethoxybenzenes
## 9 9 5 Disaccharides
## 10 10 4 Fatty acyl glycosides
## # … with 19 more rowsDraw and save as .png or .pdf image files.
p <- visualize(mcn, "parent")
ggsave(f5.61 <- paste0(tmp, "/parent_nebula.png"), p)
pdf(f5.62 <- paste0(tmp, "/child_nebula.pdf"), 12, 14)
visualize_all(mcn)
dev.off()
Figure 1: Parent-Nebula
Figure 2: Child-Nebulae
In general, Parent-Nebulae (Fig. 1) is too informative to show, so Child-Nebulae (Fig. 2) was used to dipict the abundant classes of features (metabolites) in a grid panel, intuitively. In a bird’s eye view of Child-Nebulae, we can obtain many characteristics of features, involving classes distribution, structure identified accuracy, as well as spectral similarity within classes.
5 Nebulae for Downstream analysis
5.1 Statistic analysis
Next we perform a statistical analysis with quantification data of the
features. Note that the SIRIUS project does not contain quantification
data of features, so our object mcn naturally does not contain
that either. We need to get it from elsewhere.
utils::untar(paste0(exfiles, "/herbal.tar.gz"), exdir = tmp)
origin <- data.table::fread(paste0(tmp, "/features.csv"))
origin <- tibble::as_tibble(origin)Now, let’s check the columns in the table.
origin## # A tibble: 2,579 × 25
## `row ID` `row m/z` `row retention time` `EU-BlANK.mzML Peak …` `EU-BlANK.mzML…` `EU-BlANK.mzML…`
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 591. 15.7 15.6 15.8 175173005.
## 2 2 194. 22.4 22.2 22.5 142181248.
## 3 3 635. 16.1 16.1 16.3 133667589.
## 4 4 547. 15.0 14.9 15.1 127952206.
## 5 5 250. 24.8 24.8 24.9 93773479.
## 6 6 299. 24.9 24.7 25.1 184652820.
## 7 7 293. 22.1 22.0 22.2 98564860.
## 8 8 242. 20.5 20.6 20.7 1671972.
## 9 9 293. 22.8 22.6 22.8 74335303.
## 10 10 191. 22.2 22.0 22.3 93057814.
## # … with 2,569 more rows, and 19 more variables: `EU-Raw3.mzML Peak RT start` <dbl>,
## # `EU-Raw3.mzML Peak RT end` <dbl>, `EU-Raw3.mzML Peak area` <dbl>,
## # `EU-Pro2.mzML Peak RT start` <dbl>, `EU-Pro2.mzML Peak RT end` <dbl>,
## # `EU-Pro2.mzML Peak area` <dbl>, `EU-Pro3.mzML Peak RT start` <dbl>,
## # `EU-Pro3.mzML Peak RT end` <dbl>, `EU-Pro3.mzML Peak area` <dbl>,
## # `EU-Raw1.mzML Peak RT start` <dbl>, `EU-Raw1.mzML Peak RT end` <dbl>,
## # `EU-Raw1.mzML Peak area` <dbl>, `EU-Raw2.mzML Peak RT start` <dbl>, …Remove the rest of the columns and keep only the columns for ID, m/z, retention time, and quantification.
quant <- dplyr::select(origin, id = 1, dplyr::contains("Peak area"))
colnames(quant) <- gsub("\\.mzML Peak area", "", colnames(quant))
quant <- dplyr::mutate(quant, .features_id = as.character(id))Create the metadata table and store it in the mcn object
along with the quantification data.
gp <- c(Blank = "EU-BlANK", Raw = "EU-Raw", Pro = "EU-Pro")
metadata <- MCnebula2:::group_strings(colnames(quant), gp, "sample")
metadata$annotation <- vapply(metadata$group, switch,
FUN.VALUE = character(1),
Blank = "methanol/water (1:1, v/v)", Raw = "Raw bark of Eucommia ulmoides Oliv.",
Pro = "Precessed (with saline water) bark of Eucommia ulmoides Oliv."
)
features_quantification(mcn) <- dplyr::select(quant, -id)
sample_metadata(mcn) <- metadataVariance analysis was used as a way to detect whether there were differences between the experimental and control groups and whether the differences were significant. Linear models are an effective tool for variance analysis, and it permit very general analyses. The ‘limma’ package6 integrates a number of functions for creating linear models and regression analysis. The statistical analysis provided in MCnebula2 is mainly built around the functions in the ‘limma’ package.
In the following we use the
binary_comparison function for variance analysis.
To accommodate the downstream analysis of gene
expression that the limma package was originally used for, we
should log2-transform and centralize this data
(use the default parameter ‘fun_norm’ of binary_comparison()).
mcn <- binary_comparison(mcn, Pro - Raw)
top.list <- top_table(statistic_set(mcn))Check the results.
top.list[[1]]## # A tibble: 2,579 × 7
## .features_id logFC AveExpr t P.Value adj.P.Val B
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1623 24.9 -9.42 196. 1.96e-11 0.00000000708 10.6
## 2 1980 24.1 -9.75 192. 2.19e-11 0.00000000708 10.5
## 3 353 -25.8 -9.47 -191. 2.21e-11 0.00000000708 10.5
## 4 2318 23.8 -9.91 190. 2.30e-11 0.00000000708 10.5
## 5 1977 24.1 -9.74 189. 2.37e-11 0.00000000708 10.5
## 6 1629 24.2 -9.71 187. 2.53e-11 0.00000000708 10.5
## 7 1644 23.9 -9.86 182. 2.88e-11 0.00000000708 10.5
## 8 1641 24.0 -9.82 181. 2.94e-11 0.00000000708 10.5
## 9 2323 23.2 -10.2 181. 2.96e-11 0.00000000708 10.5
## 10 2339 23.2 -10.1 180. 3.03e-11 0.00000000708 10.5
## # … with 2,569 more rows5.2 Set tracer in Child-Nebulae
Tracking top features obtained by Variance analysis in Nebulae provides insight not only into the chemical classes of these top features, but also into other features (may be analogous metabolites). Other features are not among the top ranked features, but they may contain key features that were missed due to algorithmic specificity. By tracking top features, it is possible to revisit all features at the overall data level.
n <- 20
tops <- select_features(mcn2,
tani.score_cutoff = 0.5, order_by_coef = 1,
togather = T
)
top20 <- tops[1:n]
palette_set(melody(mcn)) <- colorRampPalette(palette_set(mcn))(n)
mcn2 <- set_tracer(mcn, top20)
mcn2 <- create_child_nebulae(mcn2)
mcn2 <- create_child_layouts(mcn2)
mcn2 <- activate_nebulae(mcn2)
mcn2 <- set_nodes_color(mcn2, use_tracer = T)Draw and save the image.
pdf(f8.2 <- paste0(tmp, "/tracer_child_nebula.pdf"), 12, 14)
visualize_all(mcn2)
dev.off()