Analysis on serum 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.
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) <- "pos"Create a temporary folder to store the output data.
tmp <- paste0(tempdir(), "/temp_data")
dir.create(tmp, F)
export_path(mcn) <- tmpIn 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 6501 features
with chemical formula identification.
This dataset was origin and processed from the research in article:
https://doi.org/10.1016/j.cell.2020.07.040
exfiles <- system.file("extdata", package = "exMCnebula2")Load the ‘.rdata’ file.
load(paste0(exfiles, "/mcn_serum6501.rdata"))
mcn <- mcn_serum6501
rm(mcn_serum6501)4.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.05, 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] "Pyrans" "Pyridines and derivatives"
## [3] "Pyrroles" "Ketones"
## [5] "Benzopyrans" "Glycerolipids"
## [7] "Indoles and derivatives" "Sugar acids and derivatives"
## [9] "Glycerophospholipids" "Steroids and steroid derivatives"
## [11] "Prenol lipids" "Branched fatty acids"
## [13] "Unsaturated fatty acids" "Hydroxy fatty acids"
## [15] "Hydroxy acids and derivatives" "Pyranones and derivatives"
## [17] "Lineolic acids and derivatives" "Steroidal glycosides"
## [19] "Acyl carnitines" "Glycinated bile acids and derivatives"
## [21] "Diacylglycerols" "Carboxylic acid salts"
## [23] "Oxosteroids" "Hydroxysteroids"
## [25] "Phosphatidylcholines" "Lysophosphatidylcholines"
## [27] "Bile acids, alcohols and derivatives" "Steroid glucuronide conjugates"
## [29] "Diterpenoids" "Bilirubins"
## [31] "Tertiary alcohols" "Beta hydroxy acids and derivatives"
## [33] "Terpene glycosides" "Hydroxy bile acids, alcohols and derivatives"
## [35] "Glycerophosphocholines" "Substituted pyrroles"
## [37] "Indoles" "Long-chain fatty acids"
## [39] "Organic cations" "Organic salts"
## [41] "Vinylogous acids"pattern <- c("stero", "fatty acid", "pyr", "hydroxy", "^orga")
dis <- unlist(lapply(pattern, grep, x = classes, ignore.case = T))
dis <- classes[dis]
dis## [1] "Steroids and steroid derivatives" "Steroidal glycosides"
## [3] "Oxosteroids" "Hydroxysteroids"
## [5] "Steroid glucuronide conjugates" "Branched fatty acids"
## [7] "Unsaturated fatty acids" "Hydroxy fatty acids"
## [9] "Long-chain fatty acids" "Pyrans"
## [11] "Pyridines and derivatives" "Pyrroles"
## [13] "Benzopyrans" "Pyranones and derivatives"
## [15] "Substituted pyrroles" "Hydroxy fatty acids"
## [17] "Hydroxy acids and derivatives" "Hydroxysteroids"
## [19] "Beta hydroxy acids and derivatives" "Hydroxy bile acids, alcohols and derivatives"
## [21] "Organic cations" "Organic salts"dis <- dis[-1]4.4 Create Nebulae
Create Nebula-Index data. This data created based on ‘stardust_classes’ data.
mcn <- backtrack_stardust(mcn, dis, remove = T)
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: 22 × 3
## seq hierarchy class.name
## <int> <dbl> <chr>
## 1 1 5 Acyl carnitines
## 2 2 4 Bile acids, alcohols and derivatives
## 3 3 4 Bilirubins
## 4 4 5 Carboxylic acid salts
## 5 5 5 Diacylglycerols
## 6 6 4 Diterpenoids
## 7 7 3 Glycerolipids
## 8 8 4 Glycerophosphocholines
## 9 9 3 Glycerophospholipids
## 10 10 5 Glycinated bile acids and derivatives
## # … with 12 more rowsDraw and save as .png or .pdf image files.
p <- visualize(mcn, "parent")
ggsave(f4.61 <- paste0(tmp, "/parent_nebula.png"), p)
pdf(f4.62 <- paste0(tmp, "/child_nebula.pdf"), 12, 14)
visualize_all(mcn)
dev.off()
Figure 1: Parent-Nebula
Figure 2: Child-Nebulae
5 Nebulae for Downstream analysis
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.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, "/serum.tar.gz"), exdir = tmp)
origin <- data.table::fread(paste0(tmp, "/serum_origin_mmc3.tsv"),
skip = 1
)
origin <- tibble::as_tibble(origin)Its original data can obtained from: https://www.cell.com/cms/10.1016/j.cell.2020.07.040/attachment/f13178d1-d1ee-4179-9d33-227a02e604f1/mmc3.xlsx. Now, let’s check the columns in the table.
origin## # A tibble: 5,280 × 225
## Unique_ID `m/z` RT Subnetwork Molecular_Informati… Spectral_Librar… Superclass Class Subclass
## <int> <dbl> <dbl> <int> <chr> <chr> <chr> <chr> <chr>
## 1 349 540. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 2 228 548. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 3 1963 1092. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 4 971 547. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 5 13 546. 6.57 45 <NA> <NA> <NA> <NA> <NA>
## 6 4146 541. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 7 854 289. 4.12 -1 <NA> <NA> <NA> <NA> <NA>
## 8 4046 1036. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 9 1374 788. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## 10 2304 789. 6.57 -1 <NA> <NA> <NA> <NA> <NA>
## # … with 5,270 more rows, and 216 more variables: Direct_Parent <chr>, Infection_pvalue <dbl>,
## # Infection_FC <dbl>, Infection_BH_CV <dbl>, Infection_BH_Sig <chr>, Mortality_pvalue <chr>,
## # Mortality_FC <chr>, Mortality_BH_CV <dbl>, Mortality_BH_Sig <chr>, EFS_Rank <int>,
## # MW_Rank <int>, Avg_Rank <dbl>, ANOVA_pvalue <dbl>, `Cluster_#` <int>, ANOVA_CV <dbl>,
## # ANOVA_BH_Sig <chr>, NN1 <dbl>, NN10 <dbl>, NN11 <dbl>, NN12 <dbl>, NN13 <dbl>, NN14 <dbl>,
## # NN15 <dbl>, NN2 <dbl>, NN3 <dbl>, NN4 <dbl>, NN5 <dbl>, NN6 <dbl>, NN7 <dbl>, NN8 <dbl>,
## # NN9 <dbl>, HN1 <dbl>, HN10 <dbl>, HN2 <dbl>, HN3 <dbl>, HN4 <dbl>, HN5 <dbl>, HN6 <dbl>, …Remove the rest of the columns and keep only the columns for ID, m/z, retention time, and quantification.
keep <- grep("^[A-Z]{2}[0-9]{1,3}$", colnames(origin))
quant <- dplyr::select(origin, oid = 1, mz = 2, rt = 3, dplyr::all_of(keep))The IDs in the data quant are
different from the IDs in the object mcn, so we need to align them
first, according to mz and rt (they originate from the same batch of samples).
In the following, we have allowed the two sets of data to be merged as
accurately as possible in the form of evaluation of score:
- Score = (1 - rt.difference / rt.tolerance) * rt.weight + (1 - mz.defference / mz.tolerance) * mz.weight
meta_col <- dplyr::select(
features_annotation(mcn), .features_id,
mz, rt.secound
)
meta_col$rt.min <- meta_col$rt.secound / 60
merged <- align_merge(meta_col, quant, ".features_id",
rt.main = "rt.min",
rt.sub = "rt"
)
merged <- dplyr::select(
merged, -mz.main, -mz.sub, -rt.min, -rt,
-rt.secound
)Due to the differences in feature detection algorithms, some of the features inevitably do not get matched.
merged## # A tibble: 3,680 × 202
## .features_id oid NN1 NN10 NN11 NN12 NN13 NN14 NN15 NN2
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2435 0.0000298 0.0000176 7.3 e-6 8.34e-6 4.62e-6 NA NA 4.99e-6
## 2 100 1614 0.000101 0.0000681 6.34e-5 3.94e-5 9.5 e-5 9.34e-5 6.68e-5 1.06e-4
## 3 1000 372 0.000623 0.000787 3.50e-4 1.85e-4 5.79e-4 4.27e-4 7.19e-4 6.70e-4
## 4 1001 1001 0.00110 0.000138 1.66e-4 1.57e-4 2.27e-4 1.87e-4 2.40e-4 6.67e-4
## 5 1002 2364 NA NA NA NA NA NA NA NA
## 6 1003 452 NA NA NA NA NA 8.49e-5 NA NA
## 7 1004 2962 0.0000214 NA NA NA NA NA NA 1.5 e-5
## 8 1005 1313 0.00000545 0.0000462 3.43e-5 1.46e-5 3.36e-5 1.36e-4 2.65e-5 6.82e-5
## 9 1006 883 0.000119 0.0000681 1.89e-4 5.92e-5 5.06e-5 8.37e-5 3.21e-5 1.47e-4
## 10 1007 2889 0.0000567 0.000283 2.35e-4 1.91e-4 4.64e-4 2.35e-4 3.43e-4 7.47e-5
## # … with 3,670 more rows, and 192 more variables: NN3 <dbl>, NN4 <dbl>, NN5 <dbl>, NN6 <dbl>,
## # NN7 <dbl>, NN8 <dbl>, NN9 <dbl>, HN1 <dbl>, HN10 <dbl>, HN2 <dbl>, HN3 <dbl>, HN4 <dbl>,
## # HN5 <dbl>, HN6 <dbl>, HN7 <dbl>, HN8 <dbl>, HN9 <dbl>, HS1 <dbl>, HS10 <dbl>, HS11 <dbl>,
## # HS12 <dbl>, HS13 <dbl>, HS14 <dbl>, HS15 <dbl>, HS16 <dbl>, HS17 <dbl>, HS18 <dbl>,
## # HS19 <dbl>, HS2 <dbl>, HS20 <dbl>, HS21 <dbl>, HS22 <dbl>, HS23 <dbl>, HS24 <dbl>, HS25 <dbl>,
## # HS26 <dbl>, HS27 <dbl>, HS28 <dbl>, HS29 <dbl>, HS3 <dbl>, HS30 <dbl>, HS31 <dbl>, HS32 <dbl>,
## # HS33 <dbl>, HS34 <dbl>, HS35 <dbl>, HS36 <dbl>, HS37 <dbl>, HS38 <dbl>, HS39 <dbl>, …Create the metadata table and store it in the mcn object
along with the quantification data.
gp <- c(NN = "^NN", HN = "^HN", HS = "^HS", HM = "^HM")
metadata <- MCnebula2:::group_strings(colnames(merged), gp, "sample")
metadata$group_name <- vapply(metadata$group, switch,
FUN.VALUE = character(1),
NN = "non-hospital & non-infected", HN = "hospital & non-infected",
HS = "hospital & survival", HM = "hospital & mortality"
)
metadata$supergroup <- vapply(metadata$group, switch,
FUN.VALUE = character(1),
NN = "control groups", HN = "control groups", HS = "infection groups",
HM = "infection groups"
)
features_quantification(mcn) <- dplyr::select(merged, -oid)
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. Note that
the quantification data in origin has been normalized.
To accommodate the downstream analysis of gene
expression that the limma package was originally used for, we
should log2-transform and centralize this data.
mcn <- binary_comparison(mcn, (HS + HM) - (NN + HN), HM - HS,
fun_norm = function(x) scale(log2(x), scale = F)
)
top.list <- top_table(statistic_set(mcn))To verify the validity of the above variance analysis, the data columns were merged to obtain the IDs from the original analysis.
top.list <- lapply(top.list, merge, y = dplyr::select(
merged,
.features_id, oid
), by = ".features_id", all.x = T, sort = F)
top.list <- lapply(top.list, tibble::as_tibble)Verify with the EFS_Rank and MW_Rank column in the origin data.
(The original authors used the two methods to rank the features.)
origin_top50 <- dplyr::filter(origin, EFS_Rank <= 50 | MW_Rank <=
50)
inter. <- lapply(top.list, function(df) {
match <- head(df, n = 50)$oid %in% origin_top50$Unique_ID
oid <- head(df, n = 50)$oid[match]
list(table.match = table(match), oid = oid)
})
lapply(inter., function(x) x$table.match)## $`(HS + HM) - (NN + HN)`
## match
## FALSE TRUE
## 48 2
##
## $`HM - HS`
## match
## FALSE TRUE
## 13 37Let’s see which compounds were identified that intersected our ranking and the original ranking of features.
inter.compound <- dplyr::filter(origin, Unique_ID %in% inter.[[2]]$oid)
table(inter.compound$Spectral_Library_Match, useNA = "if")##
## Decanoyl-L-carnitine L-THYROXINE <NA>
## 1 1 27Interestingly, these two compounds were critical compounds in the original study.
5.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 <- 50
tops <- unique(unlist(lapply(top.list, function(df) df$.features_id[1:n])))
palette_set(melody(mcn)) <- colorRampPalette(palette_set(mcn))(length(tops))
mcn2 <- set_tracer(mcn, tops)
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(f6.2 <- paste0(tmp, "/tracer_child_nebula.pdf"), 12, 14)
visualize_all(mcn2)
dev.off()