diff --git a/DESCRIPTION b/DESCRIPTION
index 1a3d9b9..4872233 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.7
+Version: 1.1.8
 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
diff --git a/NAMESPACE b/NAMESPACE
index 93c6a16..6431060 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ export(compareSpectriPy)
 export(defaultSpectraVariableMapping)
 export(filterSpectriPy)
 export(normalize_intensities)
+export(pyspec_copy_on_replace)
 export(pyspec_to_rspec)
 export(reindex)
 export(remove_peaks_around_precursor_mz)
@@ -49,6 +50,7 @@ importFrom(reticulate,py_dict)
 importFrom(reticulate,py_eval)
 importFrom(reticulate,py_get_attr)
 importFrom(reticulate,py_has_attr)
+importFrom(reticulate,py_list_attributes)
 importFrom(reticulate,py_require)
 importFrom(reticulate,py_run_string)
 importFrom(reticulate,py_set_attr)
diff --git a/NEWS.md b/NEWS.md
index 1af781e..3849d60 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,10 @@
 # SpectriPy 1.1
 
+## Changes in 1.1.8
+
+- Implement `pyspec_copy_on_replace()` to enable copying/cloning the MS data
+  in Python for any data replacement operation performed in R (issue #91).
+
 ## Changes in 1.1.7
 
 - Improved `$<-` operation: only replace the values for the selected spectra
diff --git a/R/MsBackendPython.R b/R/MsBackendPython.R
index e94446c..106313e 100644
--- a/R/MsBackendPython.R
+++ b/R/MsBackendPython.R
@@ -63,7 +63,8 @@
 #' (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).
+#' same Python variable). See [pyspec_copy_on_replace()] for a strategy to
+#' avoid such data corruption.
 #'
 #' 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
@@ -249,6 +250,9 @@
 #'
 #' @name MsBackendPy
 #'
+#' @seealso [pyspec_copy_on_replace()] to ensure MS data in Python gets copied
+#'     on any data replacement operation.
+#'
 #' @return
 #'
 #' See description of individual functions for their return values.
@@ -565,6 +569,7 @@ setReplaceMethod("spectraData", "MsBackendPy", function(object, value) {
     if (nrow(value) != length(object))
         stop("The number of rows of the provided 'DataFrame' has to match ",
              "the number of spectra (=", length(object), ").")
+    object <- .backend_copy_on_replace(object)
     svm <- object@spectraVariableMapping
     ## keep only variable mapping for which we have values in 'value'
     svm <- svm[names(svm) %in% colnames(value)]
@@ -612,6 +617,7 @@ setMethod(
 #' @rdname MsBackendPy
 setReplaceMethod("peaksData", "MsBackendPy", function(object, value) {
     Spectra:::.check_peaks_data_value(value, length(object))
+    object <- .backend_copy_on_replace(object)
     object <- .py_realize_subset(object)
     switch(object@py_lib,
            matchms = .py_matchms_replace(object@py_var, value, "peaks"),
@@ -639,6 +645,7 @@ setMethod("$", "MsBackendPy", function(x, name) {
 #' @importMethodsFrom Spectra peaksVariables
 setReplaceMethod("$", "MsBackendPy", function(x, name, value) {
     if (!length(x)) return(x)
+    x <- .backend_copy_on_replace(x)
     x <- .py_realize_subset(x)
     if (x@py_lib == "matchms") {
         if (name %in% c("mz", "intensity")) {
@@ -723,6 +730,7 @@ setMethod("intensity", "MsBackendPy", function(object) {
 #' @rdname MsBackendPy
 setReplaceMethod("intensity", "MsBackendPy", function(object, value) {
     .check_mz_intensity(value, length(object), lengths(object))
+    object <- .backend_copy_on_replace(object)
     object <- .py_realize_subset(object)
     switch(object@py_lib,
            matchms = .py_matchms_replace(
@@ -762,6 +770,7 @@ setMethod("mz", "MsBackendPy", function(object) {
 #' @rdname MsBackendPy
 setReplaceMethod("mz", "MsBackendPy", function(object, value) {
     .check_mz_intensity(value, length(object), lengths(object))
+    object <- .backend_copy_on_replace(object)
     object <- .py_realize_subset(object)
     switch(object@py_lib,
            matchms = .py_matchms_replace(
@@ -830,7 +839,6 @@ setMethod("spectraNames", "MsBackendPy", function(object) {
 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
 })
@@ -981,7 +989,8 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
 #'
 #' @noRd
 .py_realize_subset <- function(x) {
-    if (length(x) && !identical(unname(x@i), seq_along(x))) {
+    if (length(x) && (.py_var_length(x) != 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"),
@@ -992,14 +1001,6 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
     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.
@@ -1247,6 +1248,10 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
                            "intensity = _tmp_var_[i], "))
 }
 
+#' Gets the number of peaks per spectrum iterating through the individual
+#' spectra
+#'
+#' @noRd
 .py_lengths <- function(s) {
     py_run_string(
         paste0("_l_ = []\n",
@@ -1255,12 +1260,198 @@ setMethod("setBackend", c("Spectra", "MsBackendPy"),
         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)))
-## }
+################################################################################
+##    Functionality related to COPY ON REPLACE
+################################################################################
+
+#' @title Copy Python MS data structure on MS data replacement operations
+#'
+#' @description
+#'
+#' Assigning a `MsBackendPy` to another variable name in R with e.g. `a <- b`
+#' with `b` being a `MsBackendPy` object results in two different `MsBackendPy`
+#' instances in R that point however to the **same** MS data structure in
+#' Python. Changing thus one of the two `MsBackendPy` either by subset or data
+#' replacement operations changes the shared MS data in Python and thus
+#' inadvertedly also the data from any other `MsBackendPy` instance in R
+#' pointing to the same variable in Python.
+#'
+#' Setting `pyspec_copy_on_replace(TRUE)` **before** any data replacement
+#' operation on a `MsBackendPy` instance causes the full MS data in Python to
+#' be copied to a **new** Phython variable before replacing any values. Thus,
+#' in situations described above, `a` and `b` would then point to two different
+#' Python variables. After the subset or replacement operation has been
+#' performed, the *copy-on-replace* setting can again be disabled with
+#' `pyspec_copy_on_replace(FALSE)`.
+#'
+#' In general, if there are multiple `MsBackendPy` instances (or `Spectra`
+#' objects) pointing to the same Python variable, then
+#' `pyspec_copy_on_replace(TRUE)` should be called before one of the replacement
+#' methods `$<-`, `intensity<-`, `mz<-`, `peaksData<-`, `spectraData<-` or
+#' `applyProcessing()` is called on one of them. After the operation,
+#' `pyspec_copy_on_replace(FALSE)` should be called.
+#'
+#' See examples below for details.
+#'
+#' Calling `pyspec_copy_on_replace()` without `TRUE` or `FALSE` returns a
+#' `logical(1)` whether *copy-on-replace* is enabled or not.
+#'
+#' The *copy-on-replace* approach avoids data corruptions by accidental
+#' modification of common Python MS data shared by different `Spectra` objects/
+#' `MsBackendPy` instances. On the downside, it comes with a slighlty lower
+#' performance (as the copy has to be cloned) and can result in potentially
+#' multiple copies of the (same) MS data in Python. *Copy-on-replace* should
+#' thus be used with care.
+#'
+#' @details
+#'
+#' For the name of the Python MS data variable a number is appended
+#' to the original variable name, ensuring that no other Python variable with
+#' the same name exists (to avoid replacing any other variable in Python). For
+#' example, if the original Python variable's name is `"my_data"`, `"my_data_1"`
+#' is used if there is no other variable with that name. If there is already
+#' another Python variable `"my_data_1"`, then `"my_data_2"` is used instead.
+#'
+#' *copy-on-replace* can also be enabled using an R option (e.g.
+#' `options(pyspec_copy_on_replace = TRUE)`) or System environment variable
+#' (e.g. `Sys.setenv(pyspec_copy_on_replace = TRUE)`).
+#'
+#' @param x `logical()`
+#'
+#' @return If `pyspec_copy_on_replace()` is called without providing an input
+#'     parameter, a `logical(1)` is returned whether the *copy-on-replace*
+#'     option is enabled or not.
+#'
+#' @author Johannes Rainer
+#'
+#' @export
+#'
+#' @examples
+#'
+#' ## Is copy-on-replace enabled?
+#' pyspec_copy_on_replace()
+#'
+#' ## Setting *copy-on-replace* to FALSE (the default)
+#' pyspec_copy_on_replace(FALSE)
+#'
+#' ## Load an example Spectra object and change to use a `MsBackendPy`
+#' library(Spectra)
+#' library(MsBackendMgf)
+#' fl <- system.file("extdata", "mgf", "test.mgf", package = "SpectriPy")
+#' s <- Spectra(fl, source = MsBackendMgf())
+#' s <- setBackend(s, MsBackendPy(), pythonVariableName = "mgf_data")
+#'
+#' ## `s` is now a `Spectra` object with the data in Python, in a variable
+#' ## `"mgf_data"`.
+#' s
+#'
+#' ## DATA CORRUPTION
+#'
+#' ## To show how data corruption can happen, we next create a subset of
+#' ## `s` and assign that to new variable `s_sub`. Both point to the same
+#' ## data in Python
+#' s_sub <- s[1:10]
+#' s_sub
+#'
+#' ## Without any data manipulation, both variables are valid
+#' intensity(s)
+#' intensity(s_sub)
+#'
+#' ## However, any data manipulation on any of the two variables will affect
+#' ## the shared MS data in Python. Below we assign retention times to the
+#' ## data subset `s_sub`. This will cause the associated Python data in
+#' ## `"mgf_data"` to be subset to the first 10 spectra and the retention times
+#' ## be added.
+#' s_sub$rtime <- 1:10 + 0.1
+#' rtime(s_sub)
+#'
+#' ## Calling `s` throws an error; the data in Python was restricted to the
+#' ## first 10 spectra, so, `s` is no longer *valid*.
+#'
+#' ## AVOID DATA COPRRUPTION WIHT COPY ON REPLACE
+#'
+#' ## If any MS data needs to be updated, replaced or added during an analysis,
+#' ## it is suggested to enable the *copy-on-replace* option to ensure that
+#' ## also the MS data in Python gets copied. We repeat the above analysis
+#' ## after enabling *copy-on-replace*:
+#' s <- Spectra(fl, source = MsBackendMgf())
+#' s <- setBackend(s, MsBackendPy(), pythonVariableName = "mgf_data")
+#'
+#' pyspec_copy_on_replace(TRUE)
+#'
+#' ## Subset `s` and add/replace retention times in the subset
+#' s_sub <- s[1:10]
+#' s_sub
+#'
+#' s_sub$rtime <- 1:10 + 0.5
+#'
+#' ## `s` and `s_sub` are now pointing to two **different** variables in Python
+#' s
+#' s_sub
+#'
+#' ## Thus, the data from `s` was not changed
+#' rtime(s)
+#'
+#' ## While the one from `s_sub` was
+#' rtime(s_sub)
+#'
+#' ## Disable *copy-on-replace* after the data replacement
+#' pyspec_copy_on_replace(FALSE)
+pyspec_copy_on_replace <- function(x = logical()) {
+    if (!length(x))
+        return(.do_copy_on_replace())
+    if (!is.logical(x) || length(x) != 1L)
+        stop("'x' must be a logical of length 1")
+    options(pyspec_copy_on_replace = x)
+}
+
+#' Whether copy-on-replace should be done or not.
+#'
+#' @noRd
+.do_copy_on_replace <- function() {
+    getOption("pyspec_copy_on_replace", default = FALSE)
+}
+
+#' If *copy-on-replace* is enabled, `.backend_copy_on_replace()` copies the
+#' MS data structure of the Python attribyte `x@py_var` to a new Python
+#' attribute and assigns its name to `x@py_var`.
+#'
+#' @details
+#'
+#' The function uses `copy.deepcopy()` since with a simple assignment only
+#' a new list with references to the original values are created.
+#'
+#' @param x `MsBackendPy`
+#'
+#' @return `MsBackendPy` with the `@py_var` updated to the new copy of the
+#'     MS data in Python (if `.do_copy_on_replace()` returns `TRUE`)
+#'
+#' @noRd
+#'
+#' @importFrom reticulate py_list_attributes
+.backend_copy_on_replace <- function(x) {
+    if (.do_copy_on_replace()) {
+        nname <- .make_unique_name(x@py_var, py_list_attributes(py))
+        py_run_string(
+            paste0("import copy\n",
+                   nname, " = copy.deepcopy(", x@py_var, ")"),
+            local = FALSE, convert = FALSE)
+        x@py_var <- nname
+    }
+    x
+}
+
+#' Checks `x` against `names` and, if it already exists in `names`, tries to
+#' rename it appending a number to make a unique/distinct name.
+#'
+#' @noRd
+.make_unique_name <- function(x, names = character()) {
+    if (x %in% names) {
+        xl <- strsplit(x, split = "_", fixed = TRUE)[[1L]]
+        if (suppressWarnings(!is.na(i <- as.integer(xl[length(xl)])))) {
+            xl[length(xl)] <- i + 1
+        } else xl <- c(xl, 1)
+        x <- paste0(xl, collapse = "_")
+        .make_unique_name(x, names)
+    } else x
+}
diff --git a/R/conversion.R b/R/conversion.R
index ba12c40..e03e4c5 100644
--- a/R/conversion.R
+++ b/R/conversion.R
@@ -696,18 +696,3 @@ 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/R/zzz.R b/R/zzz.R
index c0e4ed3..c98168b 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -39,6 +39,14 @@ spectrum_utils <- NULL
     if ((res <- Sys.getenv("RETICULATE_PYTHON_ENV")) != "")
         packageStartupMessage("Using Python environment defined by ",
                               "'RETICULATE_PYTHON_ENV': ", res)
+    ## set pyspec_copy_on_replace
+    cor <- Sys.getenv("pyspec_copy_on_replace")
+    if (cor %in% c("0", "1"))
+        cor <- as.integer(cor)
+    cor <- as.logical(cor)
+    ## Only set if not set in options.
+    if (!is.na(cor) && !length(getOption("pyspec_copy_on_replace")))
+        options(pyspec_copy_on_replace = cor)
 }
 
 .is_spectripy_use_system <- function() {
diff --git a/man/MsBackendPy.Rd b/man/MsBackendPy.Rd
index 57b76d8..6c062a9 100644
--- a/man/MsBackendPy.Rd
+++ b/man/MsBackendPy.Rd
@@ -152,10 +152,21 @@ The \code{MsBackendPy} keeps only a reference to the MS data in Python (i.e. the
 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 \code{MsBackendPy} is accessed and translated on-the-fly from the Python
-variable. The \code{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 \code{MsBackendPy} instances
-pointing to this variable.
+variable. The \code{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
+\code{MsBackendPy} instances pointing to this variable.
+
+All \strong{Subset} operations on a \code{MsBackendPy} with e.g. \code{[} are \emph{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 \verb{$<-} or \verb{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
+\code{MsBackendPy} instances use to the same data in Python (i.e., refer to the
+same Python variable). See \code{\link[=pyspec_copy_on_replace]{pyspec_copy_on_replace()}} for a strategy to
+avoid such data corruption.
 
 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
@@ -414,6 +425,10 @@ spectraData(be)
 ## Extract the precursor m/z
 be$precursorMz
 }
+\seealso{
+\code{\link[=pyspec_copy_on_replace]{pyspec_copy_on_replace()}} to ensure MS data in Python gets copied
+on any data replacement operation.
+}
 \author{
 Johannes Rainer and the EuBIC hackathon team
 }
diff --git a/man/pyspec_copy_on_replace.Rd b/man/pyspec_copy_on_replace.Rd
new file mode 100644
index 0000000..113a202
--- /dev/null
+++ b/man/pyspec_copy_on_replace.Rd
@@ -0,0 +1,139 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/MsBackendPython.R
+\name{pyspec_copy_on_replace}
+\alias{pyspec_copy_on_replace}
+\title{Copy Python MS data structure on MS data replacement operations}
+\usage{
+pyspec_copy_on_replace(x = logical())
+}
+\arguments{
+\item{x}{\code{logical()}}
+}
+\value{
+If \code{pyspec_copy_on_replace()} is called without providing an input
+parameter, a \code{logical(1)} is returned whether the \emph{copy-on-replace}
+option is enabled or not.
+}
+\description{
+Assigning a \code{MsBackendPy} to another variable name in R with e.g. \code{a <- b}
+with \code{b} being a \code{MsBackendPy} object results in two different \code{MsBackendPy}
+instances in R that point however to the \strong{same} MS data structure in
+Python. Changing thus one of the two \code{MsBackendPy} either by subset or data
+replacement operations changes the shared MS data in Python and thus
+inadvertedly also the data from any other \code{MsBackendPy} instance in R
+pointing to the same variable in Python.
+
+Setting \code{pyspec_copy_on_replace(TRUE)} \strong{before} any data replacement
+operation on a \code{MsBackendPy} instance causes the full MS data in Python to
+be copied to a \strong{new} Phython variable before replacing any values. Thus,
+in situations described above, \code{a} and \code{b} would then point to two different
+Python variables. After the subset or replacement operation has been
+performed, the \emph{copy-on-replace} setting can again be disabled with
+\code{pyspec_copy_on_replace(FALSE)}.
+
+In general, if there are multiple \code{MsBackendPy} instances (or \code{Spectra}
+objects) pointing to the same Python variable, then
+\code{pyspec_copy_on_replace(TRUE)} should be called before one of the replacement
+methods \verb{$<-}, \verb{intensity<-}, \verb{mz<-}, \verb{peaksData<-}, \verb{spectraData<-} or
+\code{applyProcessing()} is called on one of them. After the operation,
+\code{pyspec_copy_on_replace(FALSE)} should be called.
+
+See examples below for details.
+
+Calling \code{pyspec_copy_on_replace()} without \code{TRUE} or \code{FALSE} returns a
+\code{logical(1)} whether \emph{copy-on-replace} is enabled or not.
+
+The \emph{copy-on-replace} approach avoids data corruptions by accidental
+modification of common Python MS data shared by different \code{Spectra} objects/
+\code{MsBackendPy} instances. On the downside, it comes with a slighlty lower
+performance (as the copy has to be cloned) and can result in potentially
+multiple copies of the (same) MS data in Python. \emph{Copy-on-replace} should
+thus be used with care.
+}
+\details{
+For the name of the Python MS data variable a number is appended
+to the original variable name, ensuring that no other Python variable with
+the same name exists (to avoid replacing any other variable in Python). For
+example, if the original Python variable's name is \code{"my_data"}, \code{"my_data_1"}
+is used if there is no other variable with that name. If there is already
+another Python variable \code{"my_data_1"}, then \code{"my_data_2"} is used instead.
+
+\emph{copy-on-replace} can also be enabled using an R option (e.g.
+\code{options(pyspec_copy_on_replace = TRUE)}) or System environment variable
+(e.g. \code{Sys.setenv(pyspec_copy_on_replace = TRUE)}).
+}
+\examples{
+
+## Is copy-on-replace enabled?
+pyspec_copy_on_replace()
+
+## Setting *copy-on-replace* to FALSE (the default)
+pyspec_copy_on_replace(FALSE)
+
+## Load an example Spectra object and change to use a `MsBackendPy`
+library(Spectra)
+library(MsBackendMgf)
+fl <- system.file("extdata", "mgf", "test.mgf", package = "SpectriPy")
+s <- Spectra(fl, source = MsBackendMgf())
+s <- setBackend(s, MsBackendPy(), pythonVariableName = "mgf_data")
+
+## `s` is now a `Spectra` object with the data in Python, in a variable
+## `"mgf_data"`.
+s
+
+## DATA CORRUPTION
+
+## To show how data corruption can happen, we next create a subset of
+## `s` and assign that to new variable `s_sub`. Both point to the same
+## data in Python
+s_sub <- s[1:10]
+s_sub
+
+## Without any data manipulation, both variables are valid
+intensity(s)
+intensity(s_sub)
+
+## However, any data manipulation on any of the two variables will affect
+## the shared MS data in Python. Below we assign retention times to the
+## data subset `s_sub`. This will cause the associated Python data in
+## `"mgf_data"` to be subset to the first 10 spectra and the retention times
+## be added.
+s_sub$rtime <- 1:10 + 0.1
+rtime(s_sub)
+
+## Calling `s` throws an error; the data in Python was restricted to the
+## first 10 spectra, so, `s` is no longer *valid*.
+
+## AVOID DATA COPRRUPTION WIHT COPY ON REPLACE
+
+## If any MS data needs to be updated, replaced or added during an analysis,
+## it is suggested to enable the *copy-on-replace* option to ensure that
+## also the MS data in Python gets copied. We repeat the above analysis
+## after enabling *copy-on-replace*:
+s <- Spectra(fl, source = MsBackendMgf())
+s <- setBackend(s, MsBackendPy(), pythonVariableName = "mgf_data")
+
+pyspec_copy_on_replace(TRUE)
+
+## Subset `s` and add/replace retention times in the subset
+s_sub <- s[1:10]
+s_sub
+
+s_sub$rtime <- 1:10 + 0.5
+
+## `s` and `s_sub` are now pointing to two **different** variables in Python
+s
+s_sub
+
+## Thus, the data from `s` was not changed
+rtime(s)
+
+## While the one from `s_sub` was
+rtime(s_sub)
+
+## Disable *copy-on-replace* after the data replacement
+pyspec_copy_on_replace(FALSE)
+}
+\author{
+Johannes Rainer
+}
diff --git a/tests/testthat/test_MsBackendPython.R b/tests/testthat/test_MsBackendPython.R
index 4d6937f..7fe8a79 100644
--- a/tests/testthat/test_MsBackendPython.R
+++ b/tests/testthat/test_MsBackendPython.R
@@ -1055,6 +1055,217 @@ test_that(".py_realize_subset works", {
     expect_equal(intensity(a), intensity(ref))
 })
 
+test_that("pyspec_copy_on_replace and .do_copy_on_replace work", {
+    expect_false(.do_copy_on_replace())
+    expect_false(pyspec_copy_on_replace())
+    pyspec_copy_on_replace(TRUE)
+    expect_true(.do_copy_on_replace())
+    expect_true(pyspec_copy_on_replace())
+    pyspec_copy_on_replace(FALSE)
+    expect_false(.do_copy_on_replace())
+    expect_false(pyspec_copy_on_replace())
+
+    expect_error(pyspec_copy_on_replace(c(FALSE, TRUE)), "logical")
+    expect_error(pyspec_copy_on_replace("4"), "logical")
+})
+
+test_that(".backend_copy_on_replace works", {
+    ## matchms
+    pyspec_copy_on_replace(TRUE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+    b <- .backend_copy_on_replace(a)
+    expect_equal(rtime(a), rtime(b))
+    expect_equal(precursorMz(a), precursorMz(b))
+    expect_equal(precursorMz(a), precursorMz(s))
+    expect_equal(b@py_var, "tmp_1")
+    d <- .backend_copy_on_replace(a)
+    expect_equal(rtime(a), rtime(d))
+    expect_equal(precursorMz(a), precursorMz(d))
+    expect_equal(d@py_var, "tmp_2")
+    d$precursorMz <- d$precursorMz + 10
+    expect_equal(precursorMz(a), precursorMz(s))
+    expect_equal(precursorMz(a) + 10, precursorMz(d))
+    pyspec_copy_on_replace(FALSE)
+    py_del_attr(py, "tmp")
+    py_del_attr(py, "tmp_1")
+    py_del_attr(py, "tmp_2")
+
+    ## spectrum_utils
+    pyspec_copy_on_replace(TRUE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    b <- .backend_copy_on_replace(a)
+    expect_equal(rtime(a), rtime(b))
+    expect_equal(precursorMz(a), precursorMz(b))
+    expect_equal(precursorMz(a), precursorMz(s))
+    expect_equal(b@py_var, "tmp_1")
+    d <- .backend_copy_on_replace(a)
+    expect_equal(rtime(a), rtime(d))
+    expect_equal(precursorMz(a), precursorMz(d))
+    expect_equal(d@py_var, "tmp_2")
+    d$precursorMz <- d$precursorMz + 10
+    expect_equal(precursorMz(a), precursorMz(s))
+    expect_equal(precursorMz(a) + 10, precursorMz(d))
+    pyspec_copy_on_replace(FALSE)
+    py_del_attr(py, "tmp")
+    py_del_attr(py, "tmp_1")
+    py_del_attr(py, "tmp_2")
+})
+
+test_that("deep copies work in real use cases with matchms", {
+    if (py_has_attr(py, "tmp")) py_del_attr(py, "tmp")
+    if (py_has_attr(py, "tmp_1")) py_del_attr(py, "tmp_1")
+    if (py_has_attr(py, "tmp_2")) py_del_attr(py, "tmp_2")
+    if (py_has_attr(py, "tmp_3")) py_del_attr(py, "tmp_3")
+    if (py_has_attr(py, "tmp_4")) py_del_attr(py, "tmp_4")
+    if (py_has_attr(py, "tmp_5")) py_del_attr(py, "tmp_5")
+    if (py_has_attr(py, "tmp_6")) py_del_attr(py, "tmp_6")
+
+    ## matchms
+    ## First errors and problems without deepcopy.
+    pyspec_copy_on_replace(FALSE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+    b <- a[1:10]
+    b$rtime <- 1:10 + 0.1
+    expect_error(show(a), "Indices are out of bound")
+    expect_equal(length(b), .py_var_length(b))
+    expect_equal(rtime(b), 1:10 + 0.1)
+    expect_equal(b@py_var, "tmp")
+
+    ## These have to work. Changing any values should not affect other *copies*
+    ## of the data.
+    pyspec_copy_on_replace(TRUE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp")@backend
+    expect_equal(a@py_var, "tmp")
+    b <- a[1:10]
+    expect_equal(b@py_var, "tmp")
+    ## $<-
+    b$rtime <- 1:10 + 0.1
+    expect_equal(a@py_var, "tmp")
+    expect_equal(b@py_var, "tmp_1")
+    expect_equal(length(b), .py_var_length(b))
+    expect_equal(rtime(b), 1:10 + 0.1)
+    expect_equal(rtime(a), rtime(s))
+    ## intensity<-
+    intensity(b) <- intensity(b) / 2
+    expect_equal(b@py_var, "tmp_2")
+    expect_equal(intensity(b), intensity(a)[1:10] / 2)
+    ## mz<-
+    mz(b) <- mz(b) + 100
+    expect_equal(b@py_var, "tmp_3")
+    expect_equal(mz(b), mz(a)[1:10] + 100)
+    ## peaksData<-
+    peaksData(b) <- peaksData(a)[1:10]
+    expect_equal(b@py_var, "tmp_4")
+    expect_equal(peaksData(b), peaksData(a)[1:10])
+    ## spectraData<-
+    sd <- spectraData(b)
+    sd$msLevel <- 3L
+    spectraData(b) <- sd
+    expect_equal(b@py_var, "tmp_5")
+    expect_equal(msLevel(b), rep(3L, 10))
+
+    ## Works also if backend is in a Spectra object
+    s_b <- Spectra(b)
+    rtime(s_b) <- 1:10 + 0.5
+    expect_equal(s_b@backend@py_var, "tmp_6")
+    expect_equal(rtime(s_b), 1:10 + 0.5)
+
+    pyspec_copy_on_replace(FALSE)
+
+    if (py_has_attr(py, "tmp")) py_del_attr(py, "tmp")
+    py_del_attr(py, "tmp_1")
+    py_del_attr(py, "tmp_2")
+    py_del_attr(py, "tmp_3")
+    py_del_attr(py, "tmp_4")
+    py_del_attr(py, "tmp_5")
+    py_del_attr(py, "tmp_6")
+})
+
+test_that("deep copies work in real use cases with spectrum_utils", {
+    if (py_has_attr(py, "tmp")) py_del_attr(py, "tmp")
+    if (py_has_attr(py, "tmp_1")) py_del_attr(py, "tmp_1")
+    if (py_has_attr(py, "tmp_2")) py_del_attr(py, "tmp_2")
+    if (py_has_attr(py, "tmp_3")) py_del_attr(py, "tmp_3")
+    if (py_has_attr(py, "tmp_4")) py_del_attr(py, "tmp_4")
+    if (py_has_attr(py, "tmp_5")) py_del_attr(py, "tmp_5")
+    if (py_has_attr(py, "tmp_6")) py_del_attr(py, "tmp_6")
+
+    ## spectrum_utils
+    ## First errors and problems without deepcopy.
+    pyspec_copy_on_replace(FALSE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    b <- a[1:10]
+    b$rtime <- 1:10 + 0.1
+    expect_error(show(a), "Indices are out of bound")
+    expect_equal(length(b), .py_var_length(b))
+    expect_equal(rtime(b), 1:10 + 0.1, tolerance = SPECTRUM_UTILS_TOLERANCE)
+    expect_equal(b@py_var, "tmp")
+
+    ## These have to work. Changing any values should not affect other *copies*
+    ## of the data.
+    pyspec_copy_on_replace(TRUE)
+    a <- setBackend(s, MsBackendPy(), pythonVariableName = "tmp",
+                    pythonLibrary = "spectrum_utils")@backend
+    expect_equal(a@py_var, "tmp")
+    b <- a[1:10]
+    expect_equal(b@py_var, "tmp")
+    ## $<-
+    b$rtime <- 1:10 + 0.1
+    expect_equal(a@py_var, "tmp")
+    expect_equal(b@py_var, "tmp_1")
+    expect_equal(length(b), .py_var_length(b))
+    expect_equal(rtime(b), 1:10 + 0.1, tolerance = SPECTRUM_UTILS_TOLERANCE)
+    expect_equal(rtime(a), rtime(s))
+    ## intensity<-
+    intensity(b) <- intensity(b) / 2
+    expect_equal(b@py_var, "tmp_2")
+    expect_equal(intensity(b), intensity(a)[1:10] / 2)
+    ## mz<-
+    mz(b) <- mz(b) + 100
+    expect_equal(b@py_var, "tmp_3")
+    expect_equal(mz(b), mz(a)[1:10] + 100)
+    ## peaksData<-
+    peaksData(b) <- peaksData(a)[1:10]
+    expect_equal(b@py_var, "tmp_4")
+    expect_equal(peaksData(b), peaksData(a)[1:10])
+    ## spectraData<-
+    sd <- spectraData(b)
+    sd$rtime <- 1:10 + 0.7
+    spectraData(b) <- sd
+    expect_equal(b@py_var, "tmp_5")
+    expect_equal(rtime(b), 1:10 + 0.7, tolerance = SPECTRUM_UTILS_TOLERANCE)
+
+    ## Works also if backend is in a Spectra object
+    s_b <- Spectra(b)
+    rtime(s_b) <- 1:10 + 0.5
+    expect_equal(s_b@backend@py_var, "tmp_6")
+    expect_equal(rtime(s_b), 1:10 + 0.5, tolerance = SPECTRUM_UTILS_TOLERANCE)
+
+    pyspec_copy_on_replace(FALSE)
+
+    py_del_attr(py, "tmp")
+    py_del_attr(py, "tmp_1")
+    py_del_attr(py, "tmp_2")
+    py_del_attr(py, "tmp_3")
+    py_del_attr(py, "tmp_4")
+    py_del_attr(py, "tmp_5")
+    py_del_attr(py, "tmp_6")
+})
+
+test_that(".make_unique_name works", {
+    n <- c("a", "b_b", "c_3", "d_3_1", "a_2", "a_1")
+    expect_equal(.make_unique_name("a", n), "a_3")
+    expect_equal(.make_unique_name("z", n), "z")
+    expect_equal(.make_unique_name("a_2", n), "a_3")
+    expect_equal(.make_unique_name("b", n), "b")
+    expect_equal(.make_unique_name("b_b", n), "b_b_1")
+    expect_equal(.make_unique_name("d_3", n), "d_3")
+    expect_equal(.make_unique_name("d_3_1", n), "d_3_2")
+    expect_equal(.make_unique_name("c_3", n), "c_4")
+    expect_equal(.make_unique_name("c_3.1", n), "c_3.1")
+})
 
 ## Comments, thoughts TODO
 ## DONE spectraData()<-: replaces the full data and allows adding/removing
@@ -1063,17 +1274,17 @@ test_that(".py_realize_subset works", {
 ##      spectraVariableMapping: we can use spectraDataMapping now to add
 ##      additional columns - missing data will be dropped - and not registered
 ##      in @spectraVariableMapping
-## 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
+## DONE 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/.
+## TODO applyProcessing(): replace the full data? is that needed? should
+##      internally call peaksData()<-
diff --git a/tests/testthat/test_zzz.R b/tests/testthat/test_zzz.R
index 11005d2..67e67e7 100644
--- a/tests/testthat/test_zzz.R
+++ b/tests/testthat/test_zzz.R
@@ -142,4 +142,42 @@ test_that(".check_environment works", {
     if (ret_env == "")
         Sys.unsetenv("RETICULATE_PYTHON_ENV")
     else Sys.setenv(RETICULATE_PYTHON_ENV=ret_env)
+
+    ## pyspec_copy_on_replace
+    cor_opt <- getOption("pyspec_copy_on_replace")
+    cor_opt_sys <- Sys.getenv("pyspec_copy_on_replace")
+    Sys.unsetenv("pyspec_copy_on_replace")
+    options(pyspec_copy_on_replace = NULL)
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_null(res)
+    Sys.setenv(pyspec_copy_on_replace = "false")
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_false(res)
+    options(pyspec_copy_on_replace = NULL)
+    Sys.setenv(pyspec_copy_on_replace = "TRUE")
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_true(res)
+    options(pyspec_copy_on_replace = NULL)
+    Sys.setenv(pyspec_copy_on_replace = "true")
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_true(res)
+    options(pyspec_copy_on_replace = NULL)
+    Sys.setenv(pyspec_copy_on_replace = "1")
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_true(res)
+    options(pyspec_copy_on_replace = FALSE)
+    Sys.setenv(pyspec_copy_on_replace = "1")
+    .check_environment()
+    res <- getOption("pyspec_copy_on_replace")
+    expect_false(res)
+
+    if (length(cor_opt))
+        options(pyspec_copy_on_replace = cor_opt)
+    if (cor_opt_sys != "")
+        Sys.setenv(pyspec_copy_on_replace = cor_opt_sys)
 })
diff --git a/vignettes/SpectriPy.qmd b/vignettes/SpectriPy.qmd
index f72b041..52afd0b 100644
--- a/vignettes/SpectriPy.qmd
+++ b/vignettes/SpectriPy.qmd
@@ -631,16 +631,109 @@ sps$SMILES |> head()
 See also the next section for more information on the mapping between
 `Spectra`'s spectra variables and *matchms* metadata.
 
+### Replacing data and ensuring data consistency
+
 The `MsBackendPy` has full *read/write* support, i.e., it allows to add new
 spectra variables or change existing spectra and/or peaks variables through the
 available replacement methods `spectraData()<-`, `peaksData()<-`,
-`intensity()<-`, `mz()<-` and `$<-`. Currently, this operations change however
-the full MS data of the associated Python variable. Thus, even if for example
-only the retention times are replaced using `$rtime<-`, the (full!) data is
-first loaded from Python to R, the retention time values are then replaced and
-finally the full data (spectra and peaks variables) are again stored to the
-associated Python variable. Replacement operations are thus quite memory
-demanding.
+`intensity()<-`, `mz()<-` and `$<-`. These operations directly change the MS
+data in the associated Python variable/data structure and special care is
+advised if multiple copies of a `Spectra` object pointing to the **same** Python
+variable exist. To avoid inadvertently changing data in other copies of a
+`Spectra` object present in R it is possible to (temporarily) enable a
+*copy-on-replace* strategy that copies (clones) the MS data of the MS data
+structure in Python to another variable before replacing the values.
+
+In the example below we first load an example data from an MGF file and change
+the backend to a `MsBackendPy` object hence translating the full MS data to
+Python.
+
+```{r}
+#' R session:
+
+fl <- system.file("extdata", "mgf", "test.mgf", package = "SpectriPy")
+s_mgf <- Spectra(fl, source = MsBackendMgf())
+s_mgf <- setBackend(s_mgf, MsBackendPy(), pythonVariableName = "mgf_data")
+s_mgf
+```
+
+The data is now stored in a Python variable with the name `"mgf_data"`. We next
+enable the *copy-on-replace* option of *SpectriPy* and make a subset of the
+data, still keeping the original `Spectra` object `s_mgf`.
+
+```{r}
+#' R session:
+
+#' Enable copy-on-replace for MsBackendPy
+pyspec_copy_on_replace(TRUE)
+
+#' Make a subset of the original data and assign that to a different Spectra
+s_mgf_sub <- s_mgf[1:10]
+```
+
+*SpectriPy* uses a *delayed* subset strategy keeping the index to the
+individual spectra in the `MsBackendPy`. Both `Spectra` objects however still
+reference the same original data in Python:
+
+```{r}
+#' R session:
+
+#' Get the name of the asscociated Python variable
+s_mgf@backend@py_var
+s_mgf_sub@backend@py_var
+
+#' The index to the spectrum objects in the Python list
+s_mgf@backend@i
+s_mgf_sub@backend@i
+```
+
+Next we assign/replace the retention times of the `s_mgf_sub` `Spectra` object.
+
+```{r}
+#' R session:
+
+s_mgf_sub$rtime <- 1:10 + 0.1
+```
+
+After this replacement operation, the `s_mgf_sub` `Spectra` object points to a
+**different** Python variable:
+
+```{r}
+#' R session:
+
+s_mgf_sub@backend@py_var
+```
+
+The *copy-on-replace* strategy of *SpectriPy* copies (clones) the MS data
+associated to the `MsBackendPy` to a new variable before replacing or updating
+any of its data. Thus, the *original* `Spectra` object `s_mgf` points to a
+Python variable with the MS data in its original state:
+
+```{r}
+#' R session:
+
+s_mgf@backend@py_var
+```
+
+Note that this strategy is only needed when copies (or subsets) of a `Spectra`
+object are created in R and data replacements are performed on only one of them,
+or if different operations are applied to each. If no longer needed, i.e., if
+each `Spectra` object present in R points to its own Python variable,
+*copy-on-replace* should be disabled again to avoid multiple unnecessary copies
+of the data.
+
+```{r}
+#' R session:
+
+pyspec_copy_on_replace(FALSE)
+```
+
+**To summarize**: use `pyspec_copy_on_replace(TRUE)` if you have multiple
+`Spectra` objects pointing to the same Python variable **before** you apply one
+of the replacement operations `$<-`, `spectraData<-`, `mz<-`, `intensity<-`,
+`peaksData<-`, or `applyProcessing()`. Use `pyspec_copy_on_replace(FALSE)`
+**after** the operation was performed.
+
 
 ## Conversion of spectra variables