Introduction to AlpsNMR
Sergio Oller, Institute for Bioengineering of Catalonia
2021-03-31
The AlpsNMR package was written with two purposes in mind:
- to help data analysts and NMR scientists to work with NMR samples.
- to help IT pipeline builders implement automated methods for preprocessing.
Functions from this package written for data analysts and NMR scientists are prefixed with nmr_, while higher level functions written for IT pipeline builders are prefixed with pipe_. The main reason why all exported functions have a prefix is to make it easy for the user to discover the functions from the package. By typing nmr_ RStudio will return the list of exported functions. In the R terminal, nmr_ followed by the tab key (⇥) twice will have the same effect. Other popular packages, follow similar approaches (e.g: forcats: fct_*, stringr: str_*).
This vignette is written for the first group. It assumes some prior basic knowledge of NMR and data analysis, as well as some basic R programming. In case you are interested in building pipelines with this package, you may want to open the file saved in this directory (run it on your computer):
pipeline_example <- system.file("pipeline-rmd", "pipeline_example.R", package = "AlpsNMR")
pipeline_example
library(AlpsNMR)
#> Loading required package: dplyr
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
#> Loading required package: future
#> Loading required package: magrittr
library(ggplot2)
Enable parallellization
This package is able to parallellize several functions through the use of the future and furrr packages. Whether to parallelize or not is left to the user that can control the parallellization by setting “plans”.
#plan(sequential()) # disable parallellization (default)
plan(multiprocess(workers = 4)) # enable parallellization with 4 workers
Loading samples
The function to read samples is called nmr_read_samples. It expects a character vector with the samples to load that can be paths to directories of Bruker format samples or paths to JDX files.
Additionally, this function can filter by pulse sequences (e.g. load only NOESY samples) or loading only metadata.
zip_files <- fs::dir_ls(MeOH_plasma_extraction_dir, glob = "*.zip")
zip_files
#> /tmp/RtmpfjPlKu/Rinst397d1394ac4c/AlpsNMR/dataset-demo/10.zip
#> /tmp/RtmpfjPlKu/Rinst397d1394ac4c/AlpsNMR/dataset-demo/20.zip
#> /tmp/RtmpfjPlKu/Rinst397d1394ac4c/AlpsNMR/dataset-demo/30.zip
dataset <- nmr_read_samples(sample_names = zip_files)
#> Warning: `future_options()` was deprecated in furrr 0.2.0.
#> Please use `furrr_options()` instead.
dataset
#> An nmr_dataset (3 samples)
As we have not added any metadata to this dataset, the only column we see is the NMRExperiment:
nmr_meta_get(dataset, groups = "external")
#> # A tibble: 3 x 1
#> NMRExperiment
#> <chr>
#> 1 10
#> 2 20
#> 3 30
Interpolation
1D NMR samples can be interpolated together, in order to arrange all the spectra into a matrix, with one row per sample. The main parameters we would need is the range of ppm values that we want to interpolate and the resolution.
We can see the ppm resolution by looking at the ppm axis of one sample:
ppm_res <- nmr_ppm_resolution(dataset)[[1]]
message("The ppm resolution is: ", format(ppm_res, digits = 2), " ppm")
#> The ppm resolution is: 0.00023 ppm
We can interpolate the dataset, obtaining an nmr_dataset_1D object:
dataset <- nmr_interpolate_1D(dataset, axis = c(min = -0.5, max = 10, by = 2.3E-4))
This operation changes the class of the object, as now the data is on a matrix. The dataset is now of class nmr_dataset_1D. The axis element is now a numeric vector and the data_1r element is a matrix.
Plotting samples
The AlpsNMR package offers the possibility to plot nmr_dataset_1D objects. Plotting many spectra with so many points is quite expensive so it is possible to include only some regions of the spectra or plot only some samples.
Use ?plot.nmr_dataset_1D to check the parameters, among them:
NMRExperiment: A character vector with the NMR experiments to plotchemshift_range: A ppm range to plot only a small region, or to reduce the resolutioninteractive: To make the plot interactive - ...: Can be used to pass additional parameters such as color = "SubjectID" that are passed as aesthetics to ggplot.
plot(dataset, NMRExperiment = c("10", "30"), chemshift_range = c(2.2, 2.8))
Creating interactive plots
The option interactive = TRUE described above has some performance limitations. As high performance workaround, you can make many plots interactive with the function plot_interactive.
This function will use WebGL technologies to create a webpage that, once opened, allows you to interact with the plot.
Due to technical limitations, these plots need to be opened manually and can’t be embedded in RMarkdown documents. Therefore, the function saves the plot in the directory for further exploration. Additionally, some old web browsers may not be able to display these interactive plots correctly.
plt <- plot(dataset, NMRExperiment = c("10", "30"), chemshift_range = c(2.2, 2.8))
plot_interactive(plt, "plot_region.html")
Exclude regions
Some regions can easily be excluded from the spectra with nmr_exclude_region. Note that the regions are fully removed and not zeroed, as using zeros complicates a lot the implementation and has little advantages.
regions_to_exclude <- list(water = c(4.6, 5), methanol = c(3.33, 3.39))
dataset <- nmr_exclude_region(dataset, exclude = regions_to_exclude)
plot(dataset, chemshift_range = c(4.2, 5.5))
Filter samples
Maybe we just want to analyze a subset of the data, e.g., only a class group or a particular gender. We can filter some samples according to their metadata as follows:
samples_10_20 <- filter(dataset, SubjectID == "Ana")
nmr_meta_get(samples_10_20, groups = "external")
#> # A tibble: 2 x 4
#> NMRExperiment SubjectID TimePoint Age
#> <chr> <chr> <chr> <dbl>
#> 1 10 Ana baseline 29
#> 2 20 Ana 3 months 29
Robust PCA for outlier detection
The AlpsNMR package includes robust PCA analysis for outlier detection. With such a small demo dataset, it is not practical to use, but check out the documentation of nmr_pca_outliers_* functions.
pca_outliers_rob <- nmr_pca_outliers_robust(dataset, ncomp = 3)
nmr_pca_outliers_plot(dataset, pca_outliers_rob)
Baseline removal
Spectra may display an unstable baseline, specially when processing blood/fecal blood/fecal samples. If so, nmr_baseline_removal subtract the baseline by means of Asymmetric Least Squares method.
See before:
plot(dataset, chemshift_range = c(3.5,3.8))
And after:
dataset = nmr_baseline_removal(dataset, lambda = 6, p = 0.01)
plot(dataset, chemshift_range = c(3.5,3.8))
Peak detection
The peak detection is performed on short spectra segments using a continuous wavelet transform. See ?nmr_detect_peaks for more information.
Our current approach relies on the use of the baseline threshold (baselineThresh) automatic calculated (see ?nmr_baseline_threshold) and the Signal to Noise Threshold (SNR.Th) to discriminate valid peaks from noise.
The combination of the baselineThresh and the SNR.Th optimizes the number of actual peaks from noise.
The advantage of the SNR.Th method is that it estimates the noise level on each spectra region independently, so in practice it can be used as a dynamic baseline threshold level.
peak_table <- nmr_detect_peaks(dataset,
nDivRange_ppm = 0.1,
scales = seq(1, 16, 2),
baselineThresh = NULL, SNR.Th = 3)
#> Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
#> numbers without specifying argument 'seed'. There is a risk that those random
#> numbers are not statistically sound and the overall results might be invalid.
#> To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
#> numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
#> 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
#> Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
#> numbers without specifying argument 'seed'. There is a risk that those random
#> numbers are not statistically sound and the overall results might be invalid.
#> To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
#> numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
#> 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
#> Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
#> numbers without specifying argument 'seed'. There is a risk that those random
#> numbers are not statistically sound and the overall results might be invalid.
#> To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
#> numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
#> 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
NMRExp_ref <- nmr_align_find_ref(dataset, peak_table)
message("Your reference is NMRExperiment ", NMRExp_ref)
#> Your reference is NMRExperiment 20
nmr_detect_peaks_plot(dataset, peak_table, NMRExperiment = "20", chemshift_range = c(3.5,3.8))
Spectra alignment
To align the sample, we use the nmr_align function, which in turn uses a hierarchical clustering method (see ?nmr_align for further details).
The maxShift_ppm limits the maximum shift allowed for the spectra.
nmr_exp_ref <- nmr_align_find_ref(dataset, peak_table)
dataset_align <- nmr_align(dataset, peak_table, nmr_exp_ref, maxShift_ppm = 0.0015, acceptLostPeak = FALSE)
plot(dataset, chemshift_range = c(3.025, 3.063))
plot(dataset_align, chemshift_range = c(3.025, 3.063))
Normalization
There are multiple normalization techniques available. The most strongly recommended is the pqn normalization, but it may not be fully reliable when the number of samples is small, as it needs a computation of the median spectra. Nevertheless, it is possible to compute it:
dataset_norm <- nmr_normalize(dataset_align, method = "pqn")
#> Warning in norm_pqn(samples[["data_1r"]]): The Probabalistic Quotient
#> Normalization requires several samples to compute the median spectra. Your
#> number of samples is low
The AlpsNMR package offers the possibility to extract additional normalization information with nmr_normalize_extra_info(dataset), to explore the normalization factors applied to each sample:
The plot shows the dispersion with respect to the median of the normalization factors, and can highlight samples with abnormaly large or small normalization factors.
diagnostic <- nmr_normalize_extra_info(dataset_norm)
diagnostic$norm_factor
#> NMRExperiment norm_factor norm_factor_norm
#> 1 10 520968669 0.8708718
#> 2 20 690154843 1.1536901
#> 3 30 598215122 1.0000000
diagnostic$plot
Peak integration
1. Integration based on peak center and width
If we want to integrate the whole spectra, we need ppm from the peak_table. See Peak detection section. The function nmr_integrate_peak_positions generates a new nmr_dataset_1D object containing the integrals from the peak_table (ppm values corresponding to detected peaks).
peak_table_integration = nmr_integrate_peak_positions(
samples = dataset_norm,
peak_pos_ppm = peak_table$ppm,
peak_width_ppm = 0.006)
#> New names:
#> * `ppm_-0.0002` -> `ppm_-0.0002...1`
#> * ppm_0.9032 -> ppm_0.9032...4
#> * ppm_0.9085 -> ppm_0.9085...5
#> * ppm_0.9196 -> ppm_0.9196...6
#> * ppm_0.9437 -> ppm_0.9437...7
#> * ...
peak_table_integration = get_integration_with_metadata(peak_table_integration)
We can also integrate with a specific peak position and some arbitrary width:
nmr_data(
nmr_integrate_peak_positions(samples = dataset_norm,
peak_pos_ppm = c(4.1925, 4.183, 4.1775, 4.17),
peak_width_ppm = 0.006)
)
#> ppm_4.1925 ppm_4.1830 ppm_4.1775 ppm_4.1700
#> 10 3.199480e-08 3.307910e-08 2.319829e-08 1.343858e-08
#> 20 2.543387e-08 1.775268e-08 1.808552e-08 9.949101e-09
#> 30 2.571741e-08 2.550170e-08 1.911236e-08 8.217628e-09
2. Integration based on peak boundaries
Imagine we only want to integrate the four peaks corresponding to the pyroglutamic acid:
pyroglutamic_acid_region <- c(4.15, 4.20)
plot(dataset_norm, chemshift_range = pyroglutamic_acid_region) +
ggplot2::ggtitle("Pyroglutamic acid region")
We define the peak regions and integrate them. Note how we can correct the baseline or not. If we correct the baseline, the limits of the integration will be connected with a straight line and that line will be used as the baseline, that will be subtracted.
pyroglutamic_acid <- list(pyroglutamic_acid1 = c(4.19, 4.195),
pyroglutamic_acid2 = c(4.18, 4.186),
pyroglutamic_acid3 = c(4.175, 4.18),
pyroglutamic_acid4 = c(4.165, 4.172))
regions_basel_corr_ds <- nmr_integrate_regions(dataset_norm, pyroglutamic_acid, fix_baseline = TRUE)
regions_basel_corr_matrix <- nmr_data(regions_basel_corr_ds)
regions_basel_corr_matrix
#> pyroglutamic_acid1 pyroglutamic_acid2 pyroglutamic_acid3 pyroglutamic_acid4
#> 10 2.679967e-08 3.307910e-08 2.441353e-08 5.588898e-08
#> 20 2.152805e-08 1.775268e-08 1.201631e-08 5.132786e-08
#> 30 2.299145e-08 2.550170e-08 1.822286e-08 5.187043e-08
regions_basel_not_corr_ds <- nmr_integrate_regions(dataset_norm, pyroglutamic_acid, fix_baseline = FALSE)
regions_basel_not_corr_matrix <- nmr_data(regions_basel_not_corr_ds)
regions_basel_not_corr_matrix
#> pyroglutamic_acid1 pyroglutamic_acid2 pyroglutamic_acid3 pyroglutamic_acid4
#> 10 2.784051e-07 3.337131e-07 2.691285e-07 3.776341e-07
#> 20 2.720955e-07 3.286697e-07 2.667187e-07 3.788178e-07
#> 30 2.698978e-07 3.248928e-07 2.624279e-07 3.671806e-07
We may plot the integral values to explore variation based on the baseline subtraction.
dplyr::bind_rows(
regions_basel_corr_matrix %>%
as.data.frame() %>%
tibble::rownames_to_column("NMRExperiment") %>%
tidyr::gather("metabolite_peak", "area", -NMRExperiment) %>%
dplyr::mutate(BaselineCorrected = TRUE),
regions_basel_not_corr_matrix %>%
as.data.frame() %>%
tibble::rownames_to_column("NMRExperiment") %>%
tidyr::gather("metabolite_peak", "area", -NMRExperiment) %>%
dplyr::mutate(BaselineCorrected = FALSE)
) %>% ggplot() + geom_point(aes(x = NMRExperiment, y = area, color = metabolite_peak)) +
facet_wrap(~BaselineCorrected)
Identification
After applying any feature selection or machine learning, Alps allows the identification of features of interest through nmr_identify_regions_blood. The function gives 3 posibilities sorted by the most probable metabolite (see nmr_identify_regions_blood for details).
ppm_to_assign <- c(4.060960203, 3.048970634,2.405935596,0.990616851,0.986520147, 1.044258467)
identification <- nmr_identify_regions_blood (ppm_to_assign)
Free experimentation
Getting the spectra and manipulating it manually
Besides all those techniques, you can easily implement your own. You can extract the raw matrix and manipulate it at will. As long as you don’t permute the rows, you can always replace the raw matrix of the nmr_dataset_1D object through the nmr_data function:
full_spectra_matrix <- nmr_data(dataset)
full_spectra_matrix[1:3, 1:6] # change it as you wish
#> -0.5 -0.49977 -0.49954 -0.49931 -0.49908 -0.49885
#> 10 1672.044 1433.898 1072.3294 1401.1787 1634.7485 1344.6377
#> 20 1927.617 2056.491 1848.7379 1584.5820 1970.5705 2469.3993
#> 30 1337.001 1029.662 645.6248 641.0098 825.0757 799.6986
nmr_data(dataset) <- full_spectra_matrix # Rewrite the matrix
Creating an nmr_dataset_1D object from a matrix
You can also create an nmr_dataset_1D object from scratch with the new_nmr_dataset_1D function:
nsamp <- 12
npoints <- 20
# Create a random spectra matrix
dummy_ppm_axis <- seq(from = 0.2, to = 10, length.out = npoints)
dummy_spectra_matrix <- matrix(runif(nsamp*npoints), nrow = nsamp, ncol = npoints)
metadata <- list(external = data.frame(NMRExperiment = paste0("Sample", 1:12),
DummyClass = c("a", "b"),
stringsAsFactors = FALSE))
your_custom_nmr_dataset_1D <- new_nmr_dataset_1D(ppm_axis = dummy_ppm_axis,
data_1r = dummy_spectra_matrix,
metadata = metadata)
your_custom_nmr_dataset_1D
#> An nmr_dataset_1D (12 samples)
plot(your_custom_nmr_dataset_1D) +
ggtitle("Of course those random values don't make much sense...")
Final thoughts
This vignette shows many of the features of the package, some features have room for improvement, others are not fully described, and the reader will need to browse the documentation. Hopefully it is a good starting point for using the package.
sessionInfo()
#> R version 4.0.4 (2021-02-15)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
#>
#> 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
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_3.3.3 AlpsNMR_3.0.6 magrittr_2.0.1 future_1.21.0 dplyr_1.0.5
#>
#> loaded via a namespace (and not attached):
#> [1] colorspace_2.0-0 ellipsis_0.3.1
#> [3] speaq_2.6.1 corpcor_1.6.9
#> [5] XVector_0.30.0 GenomicRanges_1.42.0
#> [7] fs_1.5.0 rstudioapi_0.13
#> [9] listenv_0.8.0 furrr_0.2.2
#> [11] farver_2.1.0 ggrepel_0.9.1
#> [13] RSpectra_0.16-0 fansi_0.4.2
#> [15] mvtnorm_1.1-1 xml2_1.3.2
#> [17] codetools_0.2-18 impute_1.64.0
#> [19] knitr_1.31 mixOmics_6.14.0
#> [21] itertools_0.1-3 jsonlite_1.7.2
#> [23] cluster_2.1.1 missForest_1.4
#> [25] compiler_4.0.4 httr_1.4.2
#> [27] assertthat_0.2.1 RcppZiggurat_0.1.6
#> [29] Matrix_1.3-2 cli_2.3.1
#> [31] htmltools_0.5.1.1 tools_4.0.4
#> [33] igraph_1.2.6 qtl_1.48-1
#> [35] gtable_0.3.0 glue_1.4.2
#> [37] GenomeInfoDbData_1.2.4 reshape2_1.4.4
#> [39] Rcpp_1.0.6 limSolve_1.5.6
#> [41] Biobase_2.50.0 cellranger_1.1.0
#> [43] jquerylib_0.1.3 vctrs_0.3.7
#> [45] baseline_1.3-1 debugme_1.1.0
#> [47] iterators_1.0.13 xfun_0.22
#> [49] stringr_1.4.0 globals_0.14.0
#> [51] ps_1.6.0 rvest_1.0.0
#> [53] lpSolve_5.6.15 lifecycle_1.0.0
#> [55] zlibbioc_1.36.0 MASS_7.3-53.1
#> [57] scales_1.1.1 doSNOW_1.0.19
#> [59] MatrixGenerics_1.2.1 parallel_4.0.4
#> [61] SummarizedExperiment_1.20.0 MassSpecWavelet_1.56.0
#> [63] SparseM_1.81 RColorBrewer_1.1-2
#> [65] yaml_2.2.1 gridExtra_2.3
#> [67] sass_0.3.1 stringi_1.5.3
#> [69] highr_0.8 S4Vectors_0.28.1
#> [71] pcaPP_1.9-73 foreach_1.5.1
#> [73] randomForest_4.6-14 BiocGenerics_0.36.0
#> [75] GenomeInfoDb_1.26.4 rlang_0.4.10
#> [77] pkgconfig_2.0.3 matrixStats_0.58.0
#> [79] bitops_1.0-6 evaluate_0.14
#> [81] lattice_0.20-41 purrr_0.3.4
#> [83] labeling_0.4.2 Rfast_2.0.1
#> [85] tidyselect_1.1.0 parallelly_1.24.0
#> [87] plyr_1.8.6 R6_2.5.0
#> [89] IRanges_2.24.1 snow_0.4-3
#> [91] generics_0.1.0 DelayedArray_0.16.3
#> [93] DBI_1.1.1 pillar_1.5.1
#> [95] withr_2.4.1 RCurl_1.98-1.3
#> [97] mQTL_1.0 tibble_3.1.0
#> [99] crayon_1.4.1 rARPACK_0.11-0
#> [101] utf8_1.2.1 ellipse_0.4.2
#> [103] rmarkdown_2.7 grid_4.0.4
#> [105] readxl_1.3.1 data.table_1.14.0
#> [107] digest_0.6.27 tidyr_1.1.3
#> [109] outliers_0.14 signal_0.7-6
#> [111] stats4_4.0.4 munsell_0.5.0
#> [113] bslib_0.2.4 quadprog_1.5-8