diff --git a/DESCRIPTION b/DESCRIPTION
index 208e2d4..1a3d9b9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: SpectriPy
 Title: Enhancing Cross-Language Mass Spectrometry Data Analysis with R and Python
-Version: 1.1.5
+Version: 1.1.7
 Description: The SpectriPy package allows integration of Python-based MS
     analysis code with the Spectra package. Spectra objects can be converted
     into Python MS data structures. In addition, SpectriPy integrates
@@ -61,7 +61,8 @@ Imports:
     MsCoreUtils,
     ProtGenerics,
     methods,
-    data.table
+    data.table,
+    snakecase
 Suggests:
     testthat,
     quarto,
diff --git a/NAMESPACE b/NAMESPACE
index e9de6f5..93c6a16 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -20,6 +20,7 @@ export(select_by_intensity)
 export(select_by_mz)
 export(setSpectraVariableMapping)
 exportMethods("spectraData<-")
+exportMethods("spectraNames<-")
 exportMethods("spectraVariableMapping<-")
 exportMethods(compareSpectriPy)
 exportMethods(filterSpectriPy)
@@ -53,10 +54,12 @@ importFrom(reticulate,py_run_string)
 importFrom(reticulate,py_set_attr)
 importFrom(reticulate,py_to_r)
 importFrom(reticulate,r_to_py)
+importFrom(snakecase,to_snake_case)
 importMethodsFrom(ProtGenerics,"intensity<-")
 importMethodsFrom(ProtGenerics,"mz<-")
 importMethodsFrom(ProtGenerics,"peaksData<-")
 importMethodsFrom(ProtGenerics,"spectraData<-")
+importMethodsFrom(ProtGenerics,"spectraNames<-")
 importMethodsFrom(ProtGenerics,acquisitionNum)
 importMethodsFrom(ProtGenerics,backendInitialize)
 importMethodsFrom(ProtGenerics,centroided)
diff --git a/NEWS.md b/NEWS.md
index bf28f92..1af781e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,26 @@
 # SpectriPy 1.1
 
+## Changes in 1.1.7
+
+- Improved `$<-` operation: only replace the values for the selected spectra
+  variable.
+- Improved `peaksData<-`, `mz<-` and `intensity<-` implementations: only
+  replace the respective data, but not the full spectra data (including
+  metadata).
+- Improved `lengths()` implementation: retrieves the number of peaks directly in
+  Python.
+- Documentation updates.
+
+## Changes in 1.1.6
+
+- Automatic renaming and remapping of upper case or *camelCase* spectra
+  variables to *snake_case* metadata fields in `backendInitialize()`,
+  `spectraData<-` and `$<-`.
+- Full support of the *Spectra* test suite.
+- Add `spectraNames<-` method for `MsBackendPy` and support getting/setting
+  spectrum names. They are stored in a metadata field (spectra variable)
+  *spectrum_name*.
+
 ## Changes in 1.1.5
 
 - Fix `$<-`, `spectraData()<-` and related replacement methods for
diff --git a/R/MsBackendPython.R b/R/MsBackendPython.R
index cd8ea7f..e94446c 100644
--- a/R/MsBackendPython.R
+++ b/R/MsBackendPython.R
@@ -36,11 +36,12 @@
 #' *matchms* and *Spectra* use different names for the same metadata/variables.
 #' To support this, *MsBackendPy* implements `spectraVariableMapping` which
 #' allows to link spectra variable names to *matchms* metadata names. This
-#' mapping can be provided with the `spectraVariableMapping` to the
-#' `backendInitialize()` function. Also, be aware that *matchms* automatically
-#' **renames** certain metadata fields: a metadata `"NAME"` will be
-#' automatically renamed to `"compound_name"`. Also, *matchms* generally
-#' supports only lower-case metadata names.
+#' mapping can be provided with the `spectraVariableMapping` parameter to the
+#' `backendInitialize()` function. Note that *matchms* only supports metadata
+#' names in lower case or *snake_case*. Spectra variables are thus
+#' automatically renamed and mapped to the lower case variants if needed. Also,
+#' be aware that *matchms* automatically **renames** certain metadata fields:
+#' a metadata `"NAME"` will be automatically renamed to `"compound_name"`.
 #' See also [setSpectraVariableMapping()] for more information.
 #'
 #' @details
@@ -49,10 +50,20 @@
 #' name of the variable in Python) as well as an index pointing to the
 #' individual spectra in Python but no other data. Any data requested from
 #' the `MsBackendPy` is accessed and translated on-the-fly from the Python
-#' variable. The `MsBackendPy` is thus an interface to the MS data, but not
-#' a data container. All changes to the MS data in the Python variable
-#' (performed e.g. in Python) immediately affect any `MsBackendPy` instances
-#' pointing to this variable.
+#' variable. The `MsBackendPy` is thus an interface to the Python MS data
+#' structure, but not a data container. All changes to the MS data in the
+#' Python variable (performed e.g. in Python) immediately affect any
+#' `MsBackendPy` instances pointing to this variable.
+#'
+#' All **Subset** operations on a `MsBackendPy` with e.g. `[` are *delayed*
+#' meaning that the data in Python is not immediately changed, but that
+#' only the index to the individual spectrum objects in Python (which is stored
+#' within the backend) is updated. Importantly, however, any replacement
+#' operation such as `$<-` or `rtime<-` will realize the subset also in Python
+#' (i.e., subset and change the order of the data in Python according to the
+#' index in R). Be aware that this might cause data corruption if two
+#' `MsBackendPy` instances use to the same data in Python (i.e., refer to the
+#' same Python variable).
 #'
 #' Special care must be taken if the MS data structure in Python is subset or
 #' its order is changed (e.g. by another process). In that case it might be
@@ -67,13 +78,7 @@
 #' the variable in Python. Thus, each time MS data is requested from the
 #' backend, it is retrieved in its **current** state.
 #' If for example data was transformed or metadata added or removed in the
-#' Python object, it immediately affects the `Spectra`/backend.
-#'
-#' Any replacement operation uses internally the `spectraData()<-` method,
-#' thus replacing/updating values for individual spectra variables or peaks
-#' variables will first load the current data from Python to R, update or
-#' replace the values and then store the full MS data again to the
-#' referenced Python attribute.
+#' Python object, it immediately affects the `Spectra` object in R.
 #'
 #' @section `MsBackendPy` methods:
 #'
@@ -125,9 +130,7 @@
 #'   structure as input parameter. Both the number of spectra and the number of
 #'   peaks must match the length of the spectra and the number of existing mass
 #'   peaks. To change the number of peaks use the `peaksData()<-` method
-#'   instead that replaces the *m/z* and intensity values at the same time.
-#'   Calling `intensity()<-` will replace the full MS data (spectra variables
-#'   as well as peaks variables) of the associated Python variable.
+#'   instead that replaces the *m/z* **and** intensity values at the same time.
 #'
 #' - `mz()`, `mz()<-`: get or replace the *m/z* values. `mz()` returns a
 #'   `NumericList` of length equal to the number of spectra with each element
@@ -135,10 +138,8 @@
 #'   `mz()<-` takes the same list-like structure as input parameter. Both the
 #'   number of spectra and the number of peaks must match the length of the
 #'   spectra and the number of existing mass peaks. To change the number of
-#'   peaks use the `peaksData()<-` method instead that replaces the *m/z* and
-#'   intensity values at the same time.
-#'   Calling `mz()<-` will replace the full MS data (spectra variables
-#'   as well as peaks variables) of the associated Python variable.
+#'   peaks use the `peaksData()<-` method instead that replaces the *m/z*
+#'   **and** intensity values at the same time.
 #'
 #' - `peaksData()`: extracts the peaks data matrices from the backend. Python
 #'   code is applied to the data structure in Python to
@@ -151,9 +152,7 @@
 #'   values) for all spectra. Parameter `value` has to be a `list`-like
 #'   structure with each element being a `numeric` matrix with one column
 #'   (named `"mz"`) containing the spectrum's *m/z* and one column (named
-#'   `"intensity"`) with the intensity values. This method will replace the
-#'   full data of the associated Python variable (i.e., both the spectra as
-#'   well as the peaks data).
+#'   `"intensity"`) with the intensity values.
 #'
 #' - `spectraData()`: extracts the spectra data from the backend. Which spectra
 #'   variables are translated and retrieved from the Python objects depends on
@@ -161,11 +160,18 @@
 #'   retrieved and added to the returned `DataFrame` (with eventually missing
 #'   *core* spectra variables filled with `NA`).
 #'
-#' - `spectraData()<-`: replaces the full spectra (+ peaks) data of the backend
-#'   with the values provided with the submitted `DataFrame`. The number of
-#'   rows of this `DataFrame` has to match the number of spectra of `object`
-#'   (i.e., being equal to `length(object)`) and the `DataFrame` must also
-#'   contain the spectras' *m/z* and intensity values.
+#' - `spectraData()<-`: replaces the **full** spectra (+ peaks) data of the
+#'   backend with the values provided with the submitted `DataFrame`. The
+#'   number of rows of this `DataFrame` has to match the number of spectra of
+#'   `object` (i.e., being equal to `length(object)`) and the `DataFrame` must
+#'   also contain the spectras' *m/z* and intensity values (in columns named
+#'   `"mz"` and `"intensity"`).
+#'
+#' - `spectraNames()`, `spectraNames()<-`: extracts, respectively, adds
+#'   (replaces) names for the individual spectra. These are stored as spectra
+#'   variable/metadata field `"spectrum_name"` in the Python representation.
+#'   This is only supported if the Python *matchms* library is used, the
+#'   *spectrum_utils* library does not support arbitrary spectra names.
 #'
 #' - `spectraVariables()`: retrieves available spectra variables, which include
 #'   the names of all metadata attributes in the `matchms.Spectrum` objects
@@ -179,10 +185,9 @@
 #'   of the expected format).
 #'
 #' - `$`, `$<-`: extract or add/replace values for a spectra variable from/in
-#'   the backend. Replacing or adding values for a spectra variable cause the
-#'   full data to be replaced. In detail, first the full data is retrieved from
-#'   Python, then the values are added/replaced and then the data is again
-#'   transferred to Python.
+#'   the backend. If the *spectrum_utils* Python library is used, only a
+#'   restricted set of spectra variables are supported, i.e., `precursorMz`,
+#'   `precursorCharge`, `rtime` and `scanIndex`.
 #'
 #' @section Additional helper and utility functions:
 #'
@@ -442,16 +447,24 @@ setMethod("backendInitialize", "MsBackendPy",
                   data <- .drop_na_core_spectra_variables(data)
                   spectraVariableMapping <- spectraVariableMapping[
                       names(spectraVariableMapping) %in% colnames(data)]
+                  snake_svm <- to_snake_case(spectraVariableMapping)
+                  if (length(idx <- which(snake_svm != spectraVariableMapping)))
+                      warning(
+                          "Spectra variables ",
+                          paste0("\"", spectraVariableMapping[idx], "\"",
+                                 collapse = ", "),
+                          " have been remapped to metadata fields ",
+                          paste0("\"", snake_svm[idx], "\"", collapse = ", "))
                   py_set_attr(
                       py, pythonVariableName,
                       switch(pythonLibrary,
                              matchms = .rspec_to_matchms_pyspec(
-                                 data, mapping = spectraVariableMapping),
+                                 data, mapping = snake_svm),
                              spectrum_utils = .rspec_to_spectrum_utils_pyspec(
                                  data, mapping = spectraVariableMapping)))
                   backendInitialize(
                       object, pythonVariableName = pythonVariableName,
-                      spectraVariableMapping = spectraVariableMapping,
+                      spectraVariableMapping = snake_svm,
                       pythonLibrary = pythonLibrary)
               }
           })
@@ -563,24 +576,6 @@ setReplaceMethod("spectraData", "MsBackendPy", function(object, value) {
                       data = value)
 })
 
-#' Helper function to choose Python names that match the core ones.
-#'
-#' @param x `character` with spectra variables.
-#'
-#' @return named `character` with names corresponding to the spectra variables
-#'     and values the respective names in Python.
-#'
-#' @noRd
-.py_variable_map <- function(x = character(), library = "matchms") {
-    cm <- spectraVariableMapping(library)
-    names(x) <- x
-    idx <- match(x, names(cm))
-    nna <- !is.na(idx)
-    if (any(nna))
-        x[nna] <- cm[idx[nna]]
-    x
-}
-
 #' @importMethodsFrom ProtGenerics peaksData
 #'
 #' @importFrom reticulate iterate
@@ -617,16 +612,11 @@ setMethod(
 #' @rdname MsBackendPy
 setReplaceMethod("peaksData", "MsBackendPy", function(object, value) {
     Spectra:::.check_peaks_data_value(value, length(object))
-    spd <- spectraData(object, union(names(spectraVariableMapping(object)),
-                                     peaksVariables(object)))
-    cns <- colnames(value[[1L]])
-    for (cn in cns) {
-        vals <- lapply(value, "[", , cn)
-        if (cn %in% c("mz", "intensity"))
-            vals <- NumericList(vals, compress = FALSE)
-        spd[[cn]] <- vals
-    }
-    spectraData(object) <- spd
+    object <- .py_realize_subset(object)
+    switch(object@py_lib,
+           matchms = .py_matchms_replace(object@py_var, value, "peaks"),
+           spectrum_utils = .py_spectrum_utils_replace(
+               object@py_var, value, "peaks"))
     object
 })
 
@@ -644,17 +634,58 @@ setMethod("$", "MsBackendPy", function(x, name) {
 
 #' @rdname MsBackendPy
 #'
+#' @importFrom snakecase to_snake_case
+#'
 #' @importMethodsFrom Spectra peaksVariables
 setReplaceMethod("$", "MsBackendPy", function(x, name, value) {
-    spd <- spectraData(
-        x, union(spectraVariables(x), peaksVariables(x)))
-    spd[[name]] <- value
-    spectraData(x) <- spd
+    if (!length(x)) return(x)
+    x <- .py_realize_subset(x)
+    if (x@py_lib == "matchms") {
+        if (name %in% c("mz", "intensity")) {
+            .check_mz_intensity(value, length(x), lengths(x))
+            .py_matchms_replace(x@py_var, value, name)
+            return(x)
+        }
+        if (length(value) == 1L) value <- rep(value, length(x))
+        if (length(value) && length(value) != length(x))
+            stop("'value' has to be of length 0, 1 or ", length(x))
+        snake_name <- to_snake_case(name)
+        names(snake_name) <- name
+        if (length(value)) {
+            .py_matchms_update_metadata(x@py_var, snake_name, value)
+            if (!any(names(x@spectraVariableMapping) == name))
+                x@spectraVariableMapping <-
+                    c(x@spectraVariableMapping, snake_name)
+        } else {
+            .py_matchms_delete_metadata(x@py_var, snake_name)
+            x@spectraVariableMapping <-
+                x@spectraVariableMapping[names(x@spectraVariableMapping)!=name]
+        }
+    } else {
+        if (name %in% c("mz", "intensity")) {
+            .check_mz_intensity(value, length(x), lengths(x))
+            .py_spectrum_utils_replace(x@py_var, value, name)
+            return(x)
+        }
+        su_name <- spectraVariableMapping("spectrum_utils")[name]
+        if (is.na(su_name))
+            stop("no metadata '", name, "'  supported by spectrum_utils. ",
+                 "Allowed are ",
+                 paste0("\"", names(spectraVariableMapping("spectrum_utils")),
+                        "\"", collapse = ", "))
+        if (length(value) == 1L) value <- rep(value, length(x))
+        if (length(value) != length(x))
+            stop("'value' has to be of length 1 or ", length(x))
+        if (!name %in% names(x@spectraVariableMapping))
+            x@spectraVariableMapping <- c(x@spectraVariableMapping, su_name)
+        .py_spectrum_utils_update_metadata(x@py_var, unname(su_name), value)
+    }
     x
 })
 
 setMethod("lengths", "MsBackendPy", function(x, use.names = FALSE) {
-    vapply(peaksData(x), nrow, NA_integer_)
+    if (!length(x)) return(integer(0))
+    .py_lengths(x@py_var)
 })
 
 #' @importMethodsFrom S4Vectors isEmpty
@@ -692,10 +723,12 @@ setMethod("intensity", "MsBackendPy", function(object) {
 #' @rdname MsBackendPy
 setReplaceMethod("intensity", "MsBackendPy", function(object, value) {
     .check_mz_intensity(value, length(object), lengths(object))
-    spd <- spectraData(object, union(spectraVariables(object),
-                                     peaksVariables(object)))
-    spd[["intensity"]] <- value
-    spectraData(object) <- spd
+    object <- .py_realize_subset(object)
+    switch(object@py_lib,
+           matchms = .py_matchms_replace(
+               object@py_var, value, "intensity"),
+           spectrum_utils = .py_spectrum_utils_replace(
+               object@py_var, value, "intensity"))
     object
 })
 
@@ -729,10 +762,12 @@ setMethod("mz", "MsBackendPy", function(object) {
 #' @rdname MsBackendPy
 setReplaceMethod("mz", "MsBackendPy", function(object, value) {
     .check_mz_intensity(value, length(object), lengths(object))
-    spd <- spectraData(object, union(names(spectraVariableMapping(object)),
-                                     peaksVariables(object)))
-    spd[["mz"]] <- value
-    spectraData(object) <- spd
+    object <- .py_realize_subset(object)
+    switch(object@py_lib,
+           matchms = .py_matchms_replace(
+               object@py_var, value, "mz"),
+           spectrum_utils = .py_spectrum_utils_replace(
+               object@py_var, value, "mz"))
     object
 })
 
@@ -782,7 +817,22 @@ setMethod("smoothed", "MsBackendPy", function(object) {
 
 #' @importMethodsFrom ProtGenerics spectraNames
 setMethod("spectraNames", "MsBackendPy", function(object) {
-    rownames(spectraData(object))
+    if (any(spectraVariables(object) == "spectrum_name"))
+        object$spectrum_name
+    else NULL
+})
+
+#' @importMethodsFrom ProtGenerics spectraNames<-
+#'
+#' @rdname MsBackendPy
+#'
+#' @export
+setReplaceMethod("spectraNames", "MsBackendPy", function(object, value) {
+    if (length(value) != length(object))
+        stop("'value' has to be of length ", length(object))
+    object <- .py_realize_subset(object)
+    object$spectrum_name <- value
+    object
 })
 
 #' @importMethodsFrom ProtGenerics tic
@@ -920,6 +970,54 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
     TRUE
 }
 
+#' Subset the Python list of spectrum objects using the index of the
+#' MsBackendPy if needed (i.e., if the backend was in any way subset before).
+#' This function replaces the original attribute/Python variable name with the
+#' data subset.
+#'
+#' @param x `MsBackendPy` to subset
+#'
+#' @return the subset `MsBackendPy`
+#'
+#' @noRd
+.py_realize_subset <- function(x) {
+    if (length(x) && !identical(unname(x@i), seq_along(x))) {
+        py_set_attr(py, "_tmp_i_", as.list(x@i - 1L))
+        py_run_string(paste0(
+            x@py_var, " = list(map(lambda i: ", x@py_var, "[i], _tmp_i_))\n"),
+            local = FALSE, convert = FALSE)
+        x@i <- seq_along(x@i)
+        on.exit(py_del_attr(py, "_tmp_i_"))
+    }
+    x
+}
+
+#' Copy a Python variable to another Python variable - in R using reticulate.
+#' This is faster than using `py_run_string()`
+#'
+#' @noRd
+.r_copy_py_attr <- function(from, to) {
+    py_set_attr(py, to, py_get_attr(py, from))
+}
+
+#' Helper function to choose Python names that match the core ones.
+#'
+#' @param x `character` with spectra variables.
+#'
+#' @return named `character` with names corresponding to the spectra variables
+#'     and values the respective names in Python.
+#'
+#' @noRd
+.py_variable_map <- function(x = character(), library = "matchms") {
+    cm <- spectraVariableMapping(library)
+    names(x) <- x
+    idx <- match(x, names(cm))
+    nna <- !is.na(idx)
+    if (any(nna))
+        x[nna] <- cm[idx[nna]]
+    x
+}
+
 #' Helper function to get the variable with the Python data. Can be either
 #' in R or in Python.
 #'
@@ -1017,3 +1115,152 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
     if (!all(vapply1l(x, is.numeric)))
         stop("elements of 'value' are expected to be numeric vectors.")
 }
+
+#' Replace or add a metadata key/field to a matchms `Spectrum` object.
+#'
+#' @param s `character(1)` with the name of the Python Spectrum list variable
+#'
+#' @param key `character(1)` with the metadata name to replace/add
+#'
+#' @param vname `character(1)` with the name of the variable containing the
+#'     values
+#'
+#' @noRd
+.py_matchms_update_metadata <- function(s, key, value) {
+    py_set_attr(py, "_tmp_var_", value)
+    py_run_string(paste0("for i in range(len(", s, ")):\n",
+                         "  ", s, "[i].set('", key, "', _tmp_var_[i])\n"),
+                  local = TRUE, convert = FALSE)
+    on.exit(py_del_attr(py, "_tmp_var_"))
+}
+
+#' Delete a metadata field from a matchms `Spectrum` object
+#'
+#' @param s `character(1)` with the name of Python Spectrum list variable
+#'
+#' @param key `character(1)` with the name of the metadata to remove
+#'
+#' @noRd
+.py_matchms_delete_metadata <- function(s, key) {
+    py_run_string(paste0("for i in range(len(", s, ")):\n",
+                         "  m = ", s, "[i].metadata\n",
+                         "  del m['", key, "']\n",
+                         "  ", s, "[i].metadata = m\n"),
+                  local = TRUE, convert = FALSE)
+}
+
+#' Replace or add a metadata key/field to a spectrum_utils spectrum object.
+#'
+#' @param s `character(1)` with the name of the Python list variable
+#'
+#' @param key `character(1)` with the metadata name to replace/add. Has to be
+#'     one of `spectraVariableMapping("spectrum_utils")`.
+#'
+#' @param vname `character(1)` with the name of the variable containing the
+#'     values
+#'
+#' @noRd
+.py_spectrum_utils_update_metadata <- function(s, key, value) {
+    py_set_attr(py, "_tmp_var_", value)
+    py_run_string(paste0("for i in range(len(", s, ")):\n",
+                         "  ", s, "[i].", key, " = _tmp_var_[i]\n"),
+                  local = TRUE, convert = FALSE)
+    on.exit(py_del_attr(py, "_tmp_var_"))
+}
+
+#' Replace the full peaks data, mz or intensity values from a list of
+#' matchms Spectrum classes
+#'
+#' @param s `character(1)` with the name of the Python list
+#'
+#' @param value `list` of two-column `numeric` matrices with the m/z and
+#'     intensity values
+#'
+#' @param what `character(1)` what to replace; can be `"peaks"`, `"mz"` or
+#'     `"intensity"`
+#'
+#' @noRd
+.py_matchms_replace <- function(s, value, what) {
+    if (!is.list(value)) value <- as.list(value)
+    py_set_attr(py, "_tmp_var_", value)
+    py_run_string(
+        paste0("import numpy as np\n",
+               "from matchms import Fragments\n",
+               "for i in range(len(", s, ")):\n",
+               "  ", s, "[i].peaks = Fragments(",
+               .py_matchms_replace_cmd(s, what),
+               ")\n"),
+        local = TRUE, convert = FALSE)
+    on.exit(py_del_attr(py, "_tmp_var_"))
+}
+
+.py_matchms_replace_cmd <- function(s, what) {
+    switch(
+        what,
+        peaks = paste0("mz = _tmp_var_[i][:, 0], ",
+                       "intensities = _tmp_var_[i][:, 1]"),
+        mz = paste0("mz = np.array(_tmp_var_[i], dtype = 'float'), ",
+                    "intensities = ", s, "[i].intensities"),
+        intensity = paste0("mz = ", s, "[i].mz, intensities = ",
+                           "np.array(_tmp_var_[i], dtype = 'float')"))
+}
+
+#' Replace the full peaks data from a list of spectrum_utils classes.
+#' `spectrum_utils.MsmsSpectrum` classes don't have setter methods, thus
+#' we need to create new objects and move metadata information from the
+#' original ones over.
+#'
+#' @param s `character(1)` with the name of the Python list
+#'
+#' @param value `list` of two-column `numeric` matrices with the m/z and
+#'     intensity values
+#'
+#' @param what `character(1)` whether to replease `"peaks"`, `"mz"` or
+#'     `"intensity"` values
+#'
+#' @noRd
+.py_spectrum_utils_replace <- function(s, value,
+                                       what = c("peaks", "mz", "intensity")) {
+    if (!is.list(value)) value <- as.list(value)
+    py_set_attr(py, "_tmp_var_", value)
+    py_run_string(
+        paste0("from spectrum_utils.spectrum import MsmsSpectrum\n",
+               "for i in range(len(", s, ")):\n",
+               "  s = ", s, "[i]\n",
+               .py_spectrum_utils_replace_cmd(s, what),
+               "identifier = s.identifier, ",
+               "precursor_mz = s.precursor_mz, ",
+               "precursor_charge = s.precursor_charge, ",
+               "retention_time = s.retention_time)\n"),
+        local = TRUE, convert = FALSE)
+    on.exit(py_del_attr(py, "_tmp_var_"))
+}
+
+.py_spectrum_utils_replace_cmd <- function(s, what) {
+    switch(
+        what,
+        peaks = paste0("  ", s, "[i] = MsmsSpectrum(mz = _tmp_var_[i][:, 0], ",
+                       "intensity = _tmp_var_[i][:, 1], "),
+        mz = paste0("  ", s, "[i] = MsmsSpectrum(mz = _tmp_var_[i], ",
+                    "intensity = s.intensity, "),
+        intensity = paste0("  ", s, "[i] = MsmsSpectrum(mz = s.mz, ",
+                           "intensity = _tmp_var_[i], "))
+}
+
+.py_lengths <- function(s) {
+    py_run_string(
+        paste0("_l_ = []\n",
+               "for i in range(len(", s, ")):\n",
+               "  _l_.append(len(", s,"[i].mz))"),
+        local = TRUE, convert = TRUE)[["_l_"]]
+}
+
+## #' This function gets the **full** data (except empty core variables) from the
+## #' backend.
+## #'
+## #' @param x `MsBackendPy` object
+## #'
+## #' @return `DataFrame` with the full data.
+## .full_data <- function(x) {
+##     spectraData(x, union(names(SpectriPy:::.py_get_metadata_names(x)), peaksVariables(x)))
+## }
diff --git a/R/conversion.R b/R/conversion.R
index e03e4c5..ba12c40 100644
--- a/R/conversion.R
+++ b/R/conversion.R
@@ -696,3 +696,18 @@ pyspec_to_rspec <- function(x, mapping = spectraVariableMapping(),
            "for i in _i_:\n",
            "  _res_.append(np.column_stack((",x,"[i].mz,",x,"[i].intensity)))")
 }
+
+## .var_name <- function(x) {
+##     i <- 0
+##     if (any(ls(envir = .GlobalEnv) == x)) {
+##         xs <- strsplit(x, "_", fixed = TRUE)[[1L]]
+##         if (!is.na(suppressWarnings(i2 <- as.integer(xs[length(xs)]))))
+##             i <- i2
+##         x <- .var_name(paste0(c(xs[1:(length(xs) - 1)], i + 1L),
+##                               collapse = "_"))
+##     }
+##     x
+## }
+
+## .var_name("a")
+## a_1 <- 2
diff --git a/man/MsBackendPy.Rd b/man/MsBackendPy.Rd
index a4b9629..57b76d8 100644
--- a/man/MsBackendPy.Rd
+++ b/man/MsBackendPy.Rd
@@ -13,6 +13,7 @@
 \alias{$<-,MsBackendPy-method}
 \alias{intensity<-,MsBackendPy-method}
 \alias{mz<-,MsBackendPy-method}
+\alias{spectraNames<-,MsBackendPy-method}
 \alias{spectraVariableMapping<-,MsBackendPy-method}
 \alias{spectraVariableMapping,MsBackendPy-method}
 \alias{reindex}
@@ -48,6 +49,8 @@
 
 \S4method{mz}{MsBackendPy}(object) <- value
 
+\S4method{spectraNames}{MsBackendPy}(object) <- value
+
 \S4method{spectraVariableMapping}{MsBackendPy}(object) <- value
 
 \S4method{spectraVariableMapping}{MsBackendPy}(object, value)
@@ -166,13 +169,7 @@ Python and the backend only references to this data through the name of
 the variable in Python. Thus, each time MS data is requested from the
 backend, it is retrieved in its \strong{current} state.
 If for example data was transformed or metadata added or removed in the
-Python object, it immediately affects the \code{Spectra}/backend.
-
-Any replacement operation uses internally the \verb{spectraData()<-} method,
-thus replacing/updating values for individual spectra variables or peaks
-variables will first load the current data from Python to R, update or
-replace the values and then store the full MS data again to the
-referenced Python attribute.
+Python object, it immediately affects the \code{Spectra} object in R.
 }
 \section{Python metadata names \emph{vs} spectra variables}{
 
@@ -182,11 +179,12 @@ a spectrum, similar to \emph{Spectra}'s spectra variables concept. However,
 \emph{matchms} and \emph{Spectra} use different names for the same metadata/variables.
 To support this, \emph{MsBackendPy} implements \code{spectraVariableMapping} which
 allows to link spectra variable names to \emph{matchms} metadata names. This
-mapping can be provided with the \code{spectraVariableMapping} to the
-\code{backendInitialize()} function. Also, be aware that \emph{matchms} automatically
-\strong{renames} certain metadata fields: a metadata \code{"NAME"} will be
-automatically renamed to \code{"compound_name"}. Also, \emph{matchms} generally
-supports only lower-case metadata names.
+mapping can be provided with the \code{spectraVariableMapping} parameter to the
+\code{backendInitialize()} function. Note that \emph{matchms} only supports metadata
+names in lower case or \emph{snake_case}. Spectra variables are thus
+automatically renamed and mapped to the lower case variants if needed. Also,
+be aware that \emph{matchms} automatically \strong{renames} certain metadata fields:
+a metadata \code{"NAME"} will be automatically renamed to \code{"compound_name"}.
 See also \code{\link[=setSpectraVariableMapping]{setSpectraVariableMapping()}} for more information.
 }
 
@@ -239,19 +237,15 @@ mass peaks per spectrum. \verb{intensity()<-} takes the same list-like
 structure as input parameter. Both the number of spectra and the number of
 peaks must match the length of the spectra and the number of existing mass
 peaks. To change the number of peaks use the \verb{peaksData()<-} method
-instead that replaces the \emph{m/z} and intensity values at the same time.
-Calling \verb{intensity()<-} will replace the full MS data (spectra variables
-as well as peaks variables) of the associated Python variable.
+instead that replaces the \emph{m/z} \strong{and} intensity values at the same time.
 \item \code{mz()}, \verb{mz()<-}: get or replace the \emph{m/z} values. \code{mz()} returns a
 \code{NumericList} of length equal to the number of spectra with each element
 being the \emph{m/z} values of the individual mass peaks per spectrum.
 \verb{mz()<-} takes the same list-like structure as input parameter. Both the
 number of spectra and the number of peaks must match the length of the
 spectra and the number of existing mass peaks. To change the number of
-peaks use the \verb{peaksData()<-} method instead that replaces the \emph{m/z} and
-intensity values at the same time.
-Calling \verb{mz()<-} will replace the full MS data (spectra variables
-as well as peaks variables) of the associated Python variable.
+peaks use the \verb{peaksData()<-} method instead that replaces the \emph{m/z}
+\strong{and} intensity values at the same time.
 \item \code{peaksData()}: extracts the peaks data matrices from the backend. Python
 code is applied to the data structure in Python to
 extract the \emph{m/z} and intensity values as a list of (numpy) arrays. These
@@ -262,19 +256,23 @@ loop in R is required to set the column names to \code{"mz"} and \code{"intensit
 values) for all spectra. Parameter \code{value} has to be a \code{list}-like
 structure with each element being a \code{numeric} matrix with one column
 (named \code{"mz"}) containing the spectrum's \emph{m/z} and one column (named
-\code{"intensity"}) with the intensity values. This method will replace the
-full data of the associated Python variable (i.e., both the spectra as
-well as the peaks data).
+\code{"intensity"}) with the intensity values.
 \item \code{spectraData()}: extracts the spectra data from the backend. Which spectra
 variables are translated and retrieved from the Python objects depends on
 the backend's \code{spectraVariableMapping()}. All metadata names defined are
 retrieved and added to the returned \code{DataFrame} (with eventually missing
 \emph{core} spectra variables filled with \code{NA}).
-\item \verb{spectraData()<-}: replaces the full spectra (+ peaks) data of the backend
-with the values provided with the submitted \code{DataFrame}. The number of
-rows of this \code{DataFrame} has to match the number of spectra of \code{object}
-(i.e., being equal to \code{length(object)}) and the \code{DataFrame} must also
-contain the spectras' \emph{m/z} and intensity values.
+\item \verb{spectraData()<-}: replaces the \strong{full} spectra (+ peaks) data of the
+backend with the values provided with the submitted \code{DataFrame}. The
+number of rows of this \code{DataFrame} has to match the number of spectra of
+\code{object} (i.e., being equal to \code{length(object)}) and the \code{DataFrame} must
+also contain the spectras' \emph{m/z} and intensity values (in columns named
+\code{"mz"} and \code{"intensity"}).
+\item \code{spectraNames()}, \verb{spectraNames()<-}: extracts, respectively, adds
+(replaces) names for the individual spectra. These are stored as spectra
+variable/metadata field \code{"spectrum_name"} in the Python representation.
+This is only supported if the Python \emph{matchms} library is used, the
+\emph{spectrum_utils} library does not support arbitrary spectra names.
 \item \code{spectraVariables()}: retrieves available spectra variables, which include
 the names of all metadata attributes in the \code{matchms.Spectrum} objects
 and the \emph{core} spectra variables \code{\link[Spectra:spectraData]{Spectra::coreSpectraVariables()}}.
@@ -284,10 +282,9 @@ and the \emph{core} spectra variables \code{\link[Spectra:spectraData]{Spectra::
 backend (see \code{\link[=setSpectraVariableMapping]{setSpectraVariableMapping()}} for details and description
 of the expected format).
 \item \code{$}, \verb{$<-}: extract or add/replace values for a spectra variable from/in
-the backend. Replacing or adding values for a spectra variable cause the
-full data to be replaced. In detail, first the full data is retrieved from
-Python, then the values are added/replaced and then the data is again
-transferred to Python.
+the backend. If the \emph{spectrum_utils} Python library is used, only a
+restricted set of spectra variables are supported, i.e., \code{precursorMz},
+\code{precursorCharge}, \code{rtime} and \code{scanIndex}.
 }
 }
 
diff --git a/tests/testthat/test_MsBackendPython.R b/tests/testthat/test_MsBackendPython.R
index 0b57d9d..4d6937f 100644
--- a/tests/testthat/test_MsBackendPython.R
+++ b/tests/testthat/test_MsBackendPython.R
@@ -115,6 +115,17 @@ test_that("backendInitialize,MsBackendPy works with providing data", {
     expect_equal(res$SMILES, s$SMILES)
     expect_equal(res$INCHI, s$INCHI)
 
+    ## spectraVariableMapping using UPPERCASE map
+    svm <- c(defaultSpectraVariableMapping(), INCHI = "INCHI",SMILES = "SMILES")
+    expect_warning(
+        res <- backendInitialize(MsBackendPy(), pythonVariableName = "ttt",
+                                 data = d, spectraVariableMapping = svm),
+        "have been remapped to")
+    expect_true(any(res@spectraVariableMapping == "inchi"))
+    expect_true(any(res@spectraVariableMapping == "smiles"))
+    expect_true(all(c("INCHI", "SMILES") %in% spectraVariables(res)))
+    expect_equal(res$INCHI, s$INCHI)
+
     ## Set core variables to NA
     d$rtime <- NA
     d$precursorMz <- NA
@@ -172,7 +183,7 @@ test_that(".py_var_length, length,MsBackendPy works", {
 
 test_that("spectraVariables and .py_get_metadata_names works", {
     be <- backendInitialize(MsBackendPy(), "s_p")
-    vars <- SpectriPy:::.py_get_metadata_names(be)
+    vars <- .py_get_metadata_names(be)
     expect_true(is.character(vars))
     expect_true(
         all(c("charge", "collision_energy", "ms_level",
@@ -540,6 +551,15 @@ test_that("spectraNames,MsBackendPy works", {
     be <- backendInitialize(MsBackendPy(), "s_p")
     res <- spectraNames(be)
     expect_equal(res, spectraNames(s@backend))
+
+    tmp <- backendInitialize(MsBackendPy(), data = spectraData(be),
+                             pythonVariableName = "delete_me")
+    expect_equal(spectraNames(tmp), NULL)
+    expect_error(spectraNames(tmp) <- 1, "'value'")
+    spectraNames(tmp) <- as.character(seq_along(tmp))
+    expect_true(any(spectraVariables(tmp) == "spectrum_name"))
+    expect_equal(tmp$spectrum_name, as.character(seq_along(tmp)))
+    expect_equal(spectraNames(tmp), as.character(seq_along(tmp)))
 })
 
 test_that("tic,MsBackendPy works", {
@@ -648,7 +668,14 @@ test_that("$<-,MsBackendPy works", {
     a$other_thing <- NULL
     expect_false(any(names(a@backend@spectraVariableMapping) == "other_thing"))
 
-    expect_error(a$other_thing <- c("a", "b", "c"), "replace")
+    expect_error(a$other_thing <- c("a", "b", "c"), "length 0, 1")
+
+    ## camelCase variable
+    a <- a@backend
+    a$myCoolVar <- 3
+    expect_true("my_cool_var" %in% a@spectraVariableMapping)
+    expect_true("myCoolVar" %in% spectraVariables(a))
+    expect_true(all(a$myCoolVar == 3))
 })
 
 test_that(".py_variable_map works", {
@@ -782,10 +809,13 @@ test_that("setBackend,Spectra,MsBackendPy works", {
 
     ## mapping.
     mp <- c(defaultSpectraVariableMapping(), TITLE = "TITLE", NAME = "NAME")
-    smod_1 <- setBackend(smod, MsBackendPy(), pythonVariableName = "smod_1",
-                         applyProcessing = TRUE, spectraVariableMapping = mp)
-    expect_true(all(c("compound_name", "title") %in% spectraVariables(smod_1)))
-    expect_equal(smod_1$title, smod$TITLE)
+    expect_warning(
+        smod_1 <- setBackend(
+            smod, MsBackendPy(), pythonVariableName = "smod_1",
+            applyProcessing = TRUE, spectraVariableMapping = mp),
+        "remapped to")
+    expect_true(all(c("compound_name", "TITLE") %in% spectraVariables(smod_1)))
+    expect_equal(smod_1$TITLE, smod$TITLE)
     expect_equal(smod_1$compound_name, smod$NAME)
 })
 
@@ -811,7 +841,7 @@ test_that("MsBackendPy of length 1 works", {
     ## Subset
     a_3 <- a[c(3L, 1L, 5L)]
     sd <- spectraData(a_3)
-    expect_equal(sd$SMILES, a$SMILES[c(3L, 1L, 5L)])
+    expect_equal(sd$SMILES, s$SMILES[c(3L, 1L, 5L)])
     a_3$other_var <- "a"
     res <- spectraData(a_3)
     expect_true(all(colnames(sd) %in% colnames(res)))
@@ -819,7 +849,7 @@ test_that("MsBackendPy of length 1 works", {
                  res[, intersect(colnames(sd), colnames(res))])
     expect_true(all(res$other_var == "a"))
 
-    a_1 <- a[3L]
+    a_1 <- a_3[3L]
     sd <- spectraData(a_1)
     a_1$other_var_2 <- "a"
     res <- spectraData(a_1)
@@ -827,8 +857,205 @@ test_that("MsBackendPy of length 1 works", {
     expect_equal(sd[, intersect(colnames(sd), colnames(res))],
                  res[, intersect(colnames(sd), colnames(res))])
     expect_equal(res$other_var_2, "a")
+    expect_equal(a_1$precursorMz, s[5]$precursorMz)
+})
+
+test_that("efficient $<- works with matchms", {
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+
+    ## .py_matchms_update_metadata
+    .py_matchms_update_metadata(a@py_var, key = "new_var", rep("A", length(a)))
+    expect_true("new_var" %in% spectraVariables(a))
+    expect_equal(a$new_var, rep("A", length(a)))
+
+    ## .py_matchms_delete_metadata
+    .py_matchms_delete_metadata(a@py_var, key = "new_var")
+    expect_false("new_var" %in% spectraVariables(a))
+
+    ## $<- length 1 and camelCase to snake_case
+    a$myVar <- 4
+    expect_true("myVar" %in% spectraVariables(a))
+    expect_equal(a$myVar, rep(4, length(a)))
+    expect_true("my_var" %in% a@spectraVariableMapping)
+
+    ## overwriting works too
+    a$myVar <- seq_along(a)
+    expect_true("myVar" %in% spectraVariables(a))
+    expect_equal(a$myVar, seq_along(a))
+    expect_true("my_var" %in% a@spectraVariableMapping)
+
+    ## removing works
+    a$myVar <- NULL
+    expect_false("myVar" %in% spectraVariables(a))
+    expect_false("my_var" %in% a@spectraVariableMapping)
+
+    expect_error(a$new_var <- c(1, 2), "0, 1")
+
+    a <- MsBackendPy()
+    a$my_var <- 3
+    expect_equal(a, MsBackendPy())
+
+    ## spectrum_utils
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    expect_true(all(is.na(a$rtime)))
+    .py_spectrum_utils_update_metadata("tmp", "retention_time",
+                                       seq_along(a) + 0.1)
+    expect_true(all(is.na(a$rtime)))
+    a@spectraVariableMapping <- c(a@spectraVariableMapping,
+                                  rtime = "retention_time")
+    expect_equal(a$rtime, seq_along(a) + 0.1,
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
+
+    a$rtime <- a$rtime + 5
+    expect_equal(a$rtime, seq_along(a) + 5.1,
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
+
+    expect_error(a$msLevel <- 3L, "no metadata")
+    expect_error(a$rtime <- NULL, "length 1 or 100")
+    expect_error(a$rtime <- c(123.2, 13.4), "length 1 or 100")
+})
+
+test_that("replacing m/z, intensity and peaks data works for matchms", {
+    ## Helper function
+    res <- .py_matchms_replace_cmd("aaa", "peaks")
+    expect_match(res, "mz = _tmp_var")
+    expect_match(res, "intensities = _tmp_var")
+    res <- .py_matchms_replace_cmd("aaa", "mz")
+    expect_match(res, "mz = np.array(_tmp_var", fixed = TRUE)
+    expect_no_match(res, "intensities = np.array(_tmp_var", fixed = TRUE)
+    res <- .py_matchms_replace_cmd("aaa", "intensity")
+    expect_no_match(res, "mz = np.array(_tmp_var", fixed = TRUE)
+    expect_match(res, "intensities = np.array(_tmp_var", fixed = TRUE)
+
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+
+    ## peaksData<-
+    pd <- peaksData(s, return.type = "list")
+    pd <- lapply(pd, function(z) {
+        z[, 2L] <- z[, 2L] / 2
+        z
+    })
+    .py_matchms_replace("tmp", pd, "peaks")
+    res <- peaksData(a)
+    expect_equal(pd, res)
+    res <- lapply(res, function(z) {
+        z[, 2L] <- z[, 2L] * 2
+        z
+    })
+    peaksData(a) <- res
+    expect_equal(peaksData(a), peaksData(s, return.type = "list"))
+
+    ## mz<-
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+    mzs <- lapply(peaksData(s, return.type = "list"), `[`, , 1L)
+    mzs <- lapply(mzs, function(z) z + 10)
+    .py_matchms_replace("tmp", mzs, "mz")
+    res <- as.list(mz(a))
+    expect_equal(mzs, res)
+    mz(a) <- lapply(mzs, function(z) z - 5)
+    res <- as.list(mz(a))
+    expect_equal(res, lapply(mzs, function(z) z -5))
+    a$mz <- lapply(mzs, function(z) z -3)
+    res <- as.list(mz(a))
+    expect_equal(res, lapply(mzs, function(z) z -3))
+
+    ## intensity<-
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+    ints <- lapply(peaksData(s, return.type = "list"), `[`, , 2L)
+    ints <- lapply(ints, function(z) z / 2)
+    .py_matchms_replace("tmp", ints, "intensity")
+    res <- as.list(intensity(a))
+    expect_equal(ints, res)
+    intensity(a) <- lapply(ints, function(z) z * 5)
+    res <- as.list(intensity(a))
+    expect_equal(res, lapply(ints, function(z) z * 5))
+    a$intensity <- lapply(ints, function(z) z -3)
+    res <- as.list(intensity(a))
+    expect_equal(res, lapply(ints, function(z) z -3))
+})
+
+test_that("replacing m/z, intensity and peaks data works for spectrum_utils", {
+    ## Helper function
+    res <- .py_spectrum_utils_replace_cmd("aaa", "peaks")
+    expect_match(res, "mz = _tmp_var")
+    expect_match(res, "intensity = _tmp_var")
+    res <- .py_spectrum_utils_replace_cmd("aaa", "mz")
+    expect_match(res, "mz = _tmp_var")
+    expect_no_match(res, "intensity = _tmp_var")
+    res <- .py_spectrum_utils_replace_cmd("aaa", "intensity")
+    expect_no_match(res, "mz = _tmp_var")
+    expect_match(res, "intensity = _tmp_var")
+
+    ## peaksData<-
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    pd <- peaksData(s, return.type = "list")
+    pd <- lapply(pd, function(z) {
+        z[, 2L] <- z[, 2L] / 2
+        z[, 1L] <- z[, 1L] - 2
+        z
+    })
+    .py_spectrum_utils_replace("tmp", pd, "peaks")
+    res <- peaksData(a)
+    expect_equal(pd, res, tolerance = SPECTRUM_UTILS_TOLERANCE)
+    res <- lapply(res, function(z) {
+        z[, 2L] <- z[, 2L] * 2
+        z[, 1L] <- z[, 1L] + 2
+        z
+    })
+    peaksData(a) <- res
+    expect_equal(peaksData(a), peaksData(s, return.type = "list"),
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
+
+    ## mz<-
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    mzs <- as.list(mz(a))
+    .py_spectrum_utils_replace(
+                    "tmp", lapply(mzs, function(z) z - 1), "mz")
+    expect_equal(as.list(mz(a)), lapply(mzs, function(z) z - 1))
+    mz(a) <- lapply(mzs, function(z) z + 10)
+    expect_equal(as.list(mz(a)), lapply(mzs, function(z) z + 10))
+    a$mz <- lapply(mzs, function(z) z / 10)
+    expect_equal(as.list(mz(a)), lapply(mzs, function(z) z / 10))
+
+    ## intensity<-
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    ints <- as.list(intensity(a))
+    SpectriPy:::.py_spectrum_utils_replace(
+                    "tmp", lapply(ints, function(z) z + 101), "intensity")
+    expect_equal(as.list(intensity(a)), lapply(ints, function(z) z + 101),
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
+    intensity(a) <- lapply(ints, function(z) z + 10.1)
+    expect_equal(as.list(intensity(a)), lapply(ints, function(z) z + 10.1),
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
+    a$intensity <- lapply(ints, function(z) z / 10)
+    expect_equal(as.list(intensity(a)), lapply(ints, function(z) z / 10),
+                 tolerance = SPECTRUM_UTILS_TOLERANCE)
 })
 
+test_that(".py_realize_subset works", {
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+
+    a <- .py_realize_subset(a)
+    expect_equal(precursorMz(a), precursorMz(s))
+
+    idx <- c(4L, 1L, 6L, 1L, 13L)
+    a@i <- idx
+
+    a <- .py_realize_subset(a)
+    expect_equal(a@i, seq_along(idx))
+    expect_equal(length(a), 5)
+    ref <- s[idx]
+    expect_equal(precursorMz(a), precursorMz(ref))
+    expect_equal(mz(a), mz(ref))
+    expect_equal(intensity(a), intensity(ref))
+})
+
+
 ## Comments, thoughts TODO
 ## DONE spectraData()<-: replaces the full data and allows adding/removing
 ##      spectra variables. number of spectra has to match.
@@ -836,11 +1063,17 @@ test_that("MsBackendPy of length 1 works", {
 ##      spectraVariableMapping: we can use spectraDataMapping now to add
 ##      additional columns - missing data will be dropped - and not registered
 ##      in @spectraVariableMapping
-## DONE $<- : replace the full data?
-## DONE peaksData()<-: replace the full data?
-## DONE mz()<-
-## DONE intensity()<-
 ## TODO applyProcessing(): replace the full data? is that needed? should
 ##      internally call peaksData()<-
 ## DONE all other replacement methods.
+## DONE $<- to replace only part of the data (not the full).
+## DONE more efficient peaksData()<-: don't replace the full data.
+## DONE more efficient mz()<-
+## DONE more efficient intensity()<-
+## TODO more efficient lengths() function
 ## TODO support the Spectra unit test suite
+## TODO selectSpectraVariables
+## TODO add copy-on-replace: backendInitialize, if variable exists, create a
+##      new one and save to that: maybe a `SpetriPyDeepCopy(TRUE/FALSE)`
+##      to enable deep copy, and then, each set operation will trigger renaming
+##      the variable/.
diff --git a/vignettes/SpectriPy.qmd b/vignettes/SpectriPy.qmd
index 1304e9e..f72b041 100644
--- a/vignettes/SpectriPy.qmd
+++ b/vignettes/SpectriPy.qmd
@@ -863,7 +863,7 @@ and copying of the MS data. See also section @sec-comments for general comments.
 
 - Since MS data can be large, it is suggested the user converts MS data mostly
   manually, and only if/when needed, using `rspec_to_pyspec()` and
-  `pyspec_to_rspec()`.
+  `pyspec_to_rspec()` - or, ideally, use a `MsBackendPy` backend.
 
 - The `rspec_to_pyspec()` and `pyspec_to_rspec()` functions **copy** the data
   while transferring. Thus, there will be eventually two (detached) copies of
@@ -877,6 +877,14 @@ and copying of the MS data. See also section @sec-comments for general comments.
   but simply interfaces the MS data in Python, any changes to this data in
   Python affect also the `Spectra` object using that backend.
 
+- Be careful when re-assigning a `Spectra` object that uses a `MsBackendPy` to
+  a new variable (e.g. `b <- a`, where `a` is such a `Spectra` object). Both
+  variables will point to the **same** Python variable and changing data/values
+  in one of the variables will affect/change the data in the other variable
+  (e.g. changing the MS level of all spectra in `b` using `b$msLevel <- 3` will
+  change the MS levels in the shared Python objects, hence, `a$msLevel` will
+  then also be `3` for all).
+
 
 # Session information