This user’s guide provides an overview of the package ASICS
. ASICS
is a
fully automated procedure to identify and quantify metabolites in \(^1\)H 1D-NMR
spectra of biological mixtures (Tardivel et al., 2017). It will enable
empowering NMR-based metabolomics by quickly and accurately helping experts to
obtain metabolic profiles. In addition to the quantification method, several
functions allowing spectrum preprocessing or statistical analyses of quantified
metabolites are available.
library(ASICS)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo
library(ASICSdata)
In this user’s guide, a subset of the public datasets from Salek et al. (2007) is used. The experiment has been designed to improve the understanding of early stage of type 2 diabetes mellitus (T2DM) development. In the dataset used, \(^1\)H-NMR human metabolome was obtained from 25 healthy volunteers and 25 T2DM patients. Raw 1D Bruker spectral data files were found in the MetaboLights database (https://www.ebi.ac.uk/metabolights/, study MTBLS1).
For most time consumming functions, a parallel implementation is available for
unix-like OS using the BiocParallel package of Bioconductor. The number of used
cores is set with the option ncores
of the corresponding functions (default
to 1
, no parallel environment).
An object of class PureLibrary
with spectra of pure metabolites is required to
perform the quantification. Such a reference library is provided in ASICS
with
191 pure metabolite spectra. These spectra are metabolite spectra used as
references for quantification: only metabolites that are present in the library
object can be identified and quantified with ASICS
.
The default library is automatically loaded at package start. Available metabolites are displayed with:
head(getSampleName(pure_library), n = 8)
## [1] "1,3-Diaminopropane" "Levoglucosan" "1-Methylhydantoin"
## [4] "1-Methyl-L-Histidine" "QuinolinicAcid" "2-AminoAdipicAcid"
## [7] "2-AminobutyricAcid" "2-Deoxyadenosine"
This library can be complemented or another library can be created with new spectra of pure metabolites. These spectra are imported from Bruker files and a new library can be created with:
pure_spectra <- importSpectraBruker(system.file("extdata", "example_library",
package = "ASICS"))
new_pure_library <- createPureLibrary(pure_spectra,
nb.protons = c(5, 4))
A new library can also be created from txt or csv files, with samples in columns
and chemical shifts in rows (see help page of createPureLibrary
function for
all details).
The newly created library can be used for quantification or merged with another one:
merged_pure_library <- c(pure_library[1:10], new_pure_library)
The PureLibrary
merged_pure_library
contains the first ten spectra of the
default library and the two newly imported spectra.
First, data are imported in a data frame from Bruker files with the
importSpectraBruker
function. These spectra are baseline corrected
(Wang et al, 2013) and normalised by the area under the curve.
spectra_data <- importSpectraBruker(system.file("extdata",
"Human_diabetes_example",
package = "ASICSdata"))
Data can also be imported from other file types with importSpectra
function.
The only constraint is to have a data frame with spectra in columns
(column names are sample names) and chemical shifts in rows (row names
correspond to the ppm grid).
diabetes <- system.file("extdata", package = "ASICSdata")
spectra_data_txt <- importSpectra(name.dir = diabetes,
name.file = "spectra_diabetes_example.txt",
type = "txt")
Several functions for the preprocessing of spectra are also available: normalisation and alignment on a reference spectrum (based on Vu et al. (2011)).
Many types of normalisation are available. By default, spectra are normalised
to a constant sum (type.norm = "CS"
). Otherwise, a normalisation method
implemented in the PepsNMR
package could be used. For example:
spectra_norm <- normaliseSpectra(spectra_data_txt, type.norm = "pqn")
## Normalisation...
The alignment algorithm is based on Vu et al. (2011). To find the reference spectrum, the LCSS similarity is used (Vlachos et al. (2002)). Then the alignment is performed using the FFT cross-correlation and a hierarchical classification.
spectra_align <- alignSpectra(spectra_norm)
Finally, from the data frame, a Spectra
object is created. This is a required
step for the quantification.
spectra_obj <- createSpectra(spectra_norm)
Identification and quantification of metabolites can now be carried out using
only the function ASICS
. This function takes approximately 2 minutes per
spectrum to run. To control randomness in the algorithm (used in the estimation
of the significativity of a given metabolite concentration), the set.seed
parameter can be used.
# part of the spectrum to exclude (water and urea)
to_exclude <- matrix(c(4.5, 5.1, 5.5, 6.5), ncol = 2, byrow = TRUE)
ASICS_results <- ASICS(spectra_obj, exclusion.areas = to_exclude)
Summary of ASICS results:
ASICS_results
## An object of class ASICSResults
## It contains 50 spectra of 31087 points.
##
## ASICS results:
## 141 metabolites are identified for this set of spectra.
## Most concentrated metabolites are: Creatinine, AceticAcid, Citrate, L-Cysteine, L-Glycine, L-GlutamicAcid
The quality of the results can be assessed by stacking the original and the reconstructed spectra on one plot. A pure metabolite spectrum can also be added for visual comparison. For example, the first spectrum with Creatinine:
plot(ASICS_results, idx = 1, xlim = c(2.8, 3.3), add.metab = "Creatinine")
Relative concentrations of identified metabolites are saved in a data frame
accessible via the get_quantification
function:
head(getQuantification(ASICS_results), 10)[, 1:2]
## ADG10003u_007 ADG10003u_008
## Creatinine 0.007767615 0.001952661
## AceticAcid 0.005654296 0.000000000
## Citrate 0.004692770 0.002339383
## L-Cysteine 0.004257320 0.002966719
## L-Glycine 0.002948796 0.000000000
## L-GlutamicAcid 0.002862505 0.000000000
## PyroglutamicAcid 0.002456795 0.000000000
## L-Proline 0.002441877 0.002791684
## L-Arabitol 0.002318606 0.001741284
## 2-AminoAdipicAcid 0.002236440 0.001156904
Some analysis functions are also available in ASICS
.
First, a design data frame is imported. In this data frame, the first column needs to correspond to sample names of all spectra.
design <- read.table(system.file("extdata", "design_diabete_example.txt",
package = "ASICSdata"), header = TRUE)
Then, a preprocessing is performed on relative quantifications: metabolites with more than 75% of null quantifications are removed as well as two samples that are considered as outliers.
analysis_data <- formatForAnalysis(getQuantification(ASICS_results),
design = design, zero.threshold = 75,
zero.group = "condition",
outliers = c("ADG10003u_007",
"ADG19007u_163"))
To explore results of ASICS quantification, a PCA can be performed on results of preprocessing with:
resPCA <- pca(analysis_data)
plot(resPCA, graph = "ind", col.ind = "condition")
plot(resPCA, graph = "var")
It is also possible to find differences between two conditions with an OPLS-DA (Thevenot et al, 2015) or with Kruskall-Wallis tests:
resOPLSDA <- oplsda(analysis_data, condition = "condition", orthoI = 1)
## OPLS-DA
## 24 samples x 100 variables and 1 response
## standard scaling of predictors and response(s)
## R2X(cum) R2Y(cum) Q2(cum) RMSEE RMSEP pre ort
## Total 0.304 0.822 0.355 0.225 0.294 1 1
resOPLSDA
## OPLS-DA performed on quantifications
## Cross validation error: 0.04
##
## Variable with the higher VIP:
## Control Group diabetes mellitus VIP influential
## 5-AminoValericAcid 0.0005837748 8.194722e-05 2.417431 TRUE
## L-Citrulline 0.0019762890 9.520208e-04 2.355666 TRUE
## L-Anserine 0.0011222063 5.337960e-04 2.338824 TRUE
## Galactitol 0.0012304333 2.774101e-04 2.135719 TRUE
## Glycerol 0.0005507758 1.200991e-03 1.985294 TRUE
## Trigonelline 0.0007805405 1.214517e-04 1.850452 TRUE
## N-(2-Furoyl)Glycine 0.0008571650 1.260515e-04 1.786568 TRUE
## PyroglutamicAcid 0.0022103241 1.682716e-03 1.731453 TRUE
## Levoglucosan 0.0010271393 6.190129e-04 1.729835 TRUE
## N-Acetyl-L-AsparticAcid 0.0003549810 1.190362e-04 1.705624 TRUE
## [...]
plot(resOPLSDA)
Results of Kruskall-Wallis tests and Benjamini-Hochberg correction:
resTests <- kruskalWallis(analysis_data, "condition")
resTests
## Kruskal-Wallis tests performed on quantifications
## Variable with the lower adjusted p-value:
##
## Feature Adjusted.p.value
## 1 Galactitol 8.330429e-05
## 2 L-Citrulline 1.820263e-04
## 3 5-AminoValericAcid 2.347467e-04
## 4 2-HydroxyphenylAceticAcid 1.099853e-03
## 5 L-Anserine 5.052104e-03
## 6 PyroglutamicAcid 7.790578e-03
## 7 Trigonelline 1.325028e-02
## 8 Glycerol 2.330794e-02
## 9 HomovanillicAcid 7.468038e-02
## 10 HippuricAcid 1.802918e-01
## [...]
plot(resTests)
An analysis on buckets can also be performed. An alignment is required before the spectrum bucketing:
spectra_align <- alignSpectra(spectra_norm)
spectra_bucket <- binning(spectra_align)
Alignment visualization:
spectra_obj_align <- createSpectra(spectra_align)
plotAlignment(spectra_obj, xlim = c(3.5,4))
plotAlignment(spectra_obj_align, xlim = c(3.5,4))
Then, a SummarizedExperiment
object is created with the formatForAnalysis
function as for quantification:
analysis_data_bucket <- formatForAnalysis(spectra_bucket, design = design,
zero.threshold = 75)
Finally, all analyses can be carried out on this object with the parameter
type.data
set to buckets
. For example, the OPLS-DA is performed with:
resOPLSDA_buckets <- oplsda(analysis_data_bucket, condition = "condition",
type.data = "buckets")
## OPLS-DA
## 26 samples x 881 variables and 1 response
## standard scaling of predictors and response(s)
## R2X(cum) R2Y(cum) Q2(cum) RMSEE RMSEP pre ort
## Total 0.369 0.988 0.616 0.0622 0.413 1 3
resOPLSDA_buckets
## OPLS-DA performed on buckets
## Cross validation error: 0.12
##
## Variable with the higher VIP:
## Control Group diabetes mellitus VIP influential
## 8.335 2.397862e-04 3.861426e-04 2.169500 TRUE
## 8.935 -1.256751e-05 -7.185148e-07 2.146483 TRUE
## 4.115 1.029309e-03 1.313057e-03 2.123301 TRUE
## 3.685 8.065244e-03 5.194614e-03 2.086658 TRUE
## 9.735 -8.507668e-06 2.194277e-06 2.069487 TRUE
## 1.035 3.999484e-04 2.995090e-04 2.060325 TRUE
## 8.625 -1.276220e-05 8.068675e-06 2.030505 TRUE
## 1.415 7.189269e-04 1.317042e-03 2.016880 TRUE
## 3.695 1.000146e-02 6.760782e-03 2.014882 TRUE
## 8.025 9.299961e-05 1.361142e-04 2.013879 TRUE
## [...]
Moreover, another plot with the median spectrum and OPLS-DA results can be
produced with the option graph = "buckets"
:
plot(resOPLSDA_buckets, graph = "buckets")
Tardivel P., Canlet C., Lefort G., Tremblay-Franco M., Debrauwer L., Concordet D., Servien R. (2017). ASICS: an automatic method for identification and quantification of metabolites in complex 1D 1H NMR spectra. Metabolomics, 13(10), 109. https://doi.org/10.1007/s11306-017-1244-5
Salek, R. M., Maguire, M. L., Bentley, E., Rubtsov, D. V., Hough, T., Cheeseman, M., … & Connor, S. C. (2007). A metabolomic comparison of urinary changes in type 2 diabetes in mouse, rat, and human. Physiological genomics, 29(2), 99-108.
Wang, K. C., Wang, S. Y., Kuo, C. H., Tseng, Y. J. (2013). Distribution-based classification method for baseline correction of metabolomic 1D proton nuclear magnetic resonance spectra. Analytical Chemistry, 85(2), 1231–1239.
Vlachos, M., Kollios, G., & Gunopulos, D. (2002). Discovering similar multidimensional trajectories. In Data Engineering, 2002. Proceedings. 18th International Conference on (pp. 673-684). IEEE.
Vu, T. N., Valkenborg, D., Smets, K., Verwaest, K. A., Dommisse, R., Lemiere, F., … & Laukens, K. (2011). An integrated workflow for robust alignment and simplified quantitative analysis of NMR spectrometry data. BMC bioinformatics, 12(1), 405.
Thevenot, E.A., Roux, A., Xu, Y., Ezan, E., Junot, C. 2015. Analysis of the human adult urinary metabolome variations with age, body mass index and gender by implementing a comprehensive workflow for univariate and OPLS statistical analyses. Journal of Proteome Research. 14, 3322-3335.
This user’s guide has been created with the following system configuration:
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ASICSdata_1.4.0 ASICS_1.99.1 BiocStyle_2.12.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 matrixStats_0.54.0
## [3] xts_0.11-2 RColorBrewer_1.1-2
## [5] webshot_0.5.1 GenomeInfoDb_1.20.0
## [7] tools_3.6.0 R6_2.4.0
## [9] KernSmooth_2.23-15 pdc_1.0.3
## [11] lazyeval_0.2.2 BiocGenerics_0.30.0
## [13] colorspace_1.4-1 manipulateWidget_0.10.0
## [15] tidyselect_0.2.5 gridExtra_2.3
## [17] compiler_3.6.0 Biobase_2.44.0
## [19] DelayedArray_0.10.0 labeling_0.3
## [21] bookdown_0.11 scales_1.0.0
## [23] ifultools_2.0-5 quadprog_1.5-7
## [25] proxy_0.4-23 dtw_1.20-1
## [27] stringr_1.4.0 digest_0.6.19
## [29] clv_0.3-2.1 rmarkdown_1.13
## [31] XVector_0.24.0 pkgconfig_2.0.2
## [33] htmltools_0.3.6 TSdist_3.5
## [35] htmlwidgets_1.3 rlang_0.3.4
## [37] shiny_1.2.0 bindr_0.1.1
## [39] zoo_1.8-6 jsonlite_1.6
## [41] crosstalk_1.0.0 BiocParallel_1.18.0
## [43] ropls_1.16.0 dplyr_0.7.8
## [45] RCurl_1.95-4.12 magrittr_1.5
## [47] GenomeInfoDbData_1.2.1 Matrix_1.2-17
## [49] Rcpp_1.0.1 munsell_0.5.0
## [51] S4Vectors_0.22.0 stringi_1.4.3
## [53] yaml_2.2.0 MASS_7.3-51.1
## [55] SummarizedExperiment_1.14.0 zlibbioc_1.30.0
## [57] plyr_1.8.4 grid_3.6.0
## [59] misc3d_0.8-4 parallel_3.6.0
## [61] promises_1.0.1 crayon_1.3.4
## [63] miniUI_0.1.1.1 lattice_0.20-38
## [65] TSclust_1.2.4 knitr_1.23
## [67] pillar_1.4.1 wmtsa_2.0-3
## [69] GenomicRanges_1.36.0 codetools_0.2-16
## [71] reshape2_1.4.3 stats4_3.6.0
## [73] glue_1.3.1 evaluate_0.14
## [75] BiocManager_1.30.4 splus2R_1.2-2
## [77] nloptr_1.2.1 httpuv_1.4.5
## [79] gtable_0.3.0 purrr_0.3.2
## [81] assertthat_0.2.1 ggplot2_3.1.1
## [83] xfun_0.7 mime_0.6
## [85] xtable_1.8-4 later_0.8.0
## [87] class_7.3-15 tibble_1.4.2
## [89] IRanges_2.18.1 locpol_0.7-0
## [91] ptw_1.9-13 longitudinalData_2.4.1
## [93] bindrcpp_0.2.2 cluster_2.0.9
## [95] rgl_0.100.19 PepsNMR_1.2.1