diff --git a/HISTORY.md b/HISTORY.md
index 6581f1f..f84441a 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -63,6 +63,10 @@ and this project adheres to
 - `proteoform_simulation_metrics()`: after running proteoform
   inference on a simulated dataset, returns one record of metrics
   plus the per-protein results table.
+- `apply_peptide_batch_effects()`: simulates peptide-specific batch effects
+  (`flavor="subset_shift"`) by shifting `n_peptides` peptides per
+  protein by `log2fc` in selected (`batches` as list) or randomly
+  drawn (`batches` as int) batches of `batch_key`.
 
 ### Fixed
 
diff --git a/docs/sphinx/source/api/utils.rst b/docs/sphinx/source/api/utils.rst
index 64ca59d..e89088a 100644
--- a/docs/sphinx/source/api/utils.rst
+++ b/docs/sphinx/source/api/utils.rst
@@ -30,7 +30,8 @@
 .. rubric:: Proteoform Group Simulation
 
 Simulate proteoform groups to evaluate how well a proteoform inference
-method recovers them in a dataset.
+method recovers them in a dataset, optionally with simulated batch
+effects to evaluate batch-aware inference.
 
 .. autosummary::
    :toctree: generated/
@@ -38,6 +39,7 @@ method recovers them in a dataset.
 
    proteopy.utils.proteoform_group_simulator
    proteopy.utils.proteoform_simulation_metrics
+   proteopy.utils.apply_peptide_batch_effects
 
 .. rubric:: Plotting
 
diff --git a/proteopy/utils/__init__.py b/proteopy/utils/__init__.py
index d628194..f7ee5d8 100644
--- a/proteopy/utils/__init__.py
+++ b/proteopy/utils/__init__.py
@@ -8,3 +8,4 @@
     proteoform_group_simulator,
     proteoform_simulation_metrics,
 )
+from .batch_effect_simulation import apply_peptide_batch_effects
diff --git a/proteopy/utils/batch_effect_simulation.py b/proteopy/utils/batch_effect_simulation.py
new file mode 100644
index 0000000..8910f01
--- /dev/null
+++ b/proteopy/utils/batch_effect_simulation.py
@@ -0,0 +1,377 @@
+"""Simulated batch effects for evaluating batch-aware methods.
+
+`apply_peptide_batch_effects` shifts the intensities of selected peptides in
+selected batches.
+"""
+
+import warnings
+from collections.abc import Sequence
+from numbers import Integral, Real
+
+import anndata as ad
+import numpy as np
+import pandas as pd
+from scipy import sparse
+
+from proteopy.utils.anndata import check_proteodata
+from proteopy.utils.array import is_log_transformed
+
+BATCH_EFFECT_UNS_KEY = "batch_effect_simulation"
+BATCH_EFFECT_VAR_KEY = "batch_effect"
+FLAVORS = ("subset_shift",)
+
+
+def _is_int(value) -> bool:
+    return isinstance(value, Integral) and not isinstance(value, bool)
+
+
+def _batch_labels(adata, batch_key):
+    """Distinct non-missing batch labels in lexicographic order."""
+    labels = pd.unique(adata.obs[batch_key].dropna())
+    return sorted(labels, key=str)
+
+
+def _validate_batches(batches, labels, batch_key):
+    if _is_int(batches):
+        if not 1 <= batches < len(labels):
+            raise ValueError(
+                f"`batches` as an int must be between 1 and "
+                f"{len(labels) - 1} (the number of batches in "
+                f".obs['{batch_key}'] minus 1), got {batches}."
+            )
+        return
+    if isinstance(batches, str) or not isinstance(
+        batches, (Sequence, np.ndarray, pd.Index)
+    ):
+        raise TypeError(
+            "`batches` must be an int or a list of batch labels, "
+            f"got {type(batches).__name__}."
+        )
+    batches = list(batches)
+    if len(batches) == 0:
+        raise ValueError("`batches` must not be empty.")
+    if len(set(batches)) != len(batches):
+        raise ValueError("`batches` must not contain duplicates.")
+    missing = [b for b in batches if b not in labels]
+    if missing:
+        raise ValueError(
+            f"Batches {missing} not found in .obs['{batch_key}']. "
+            f"Available: {labels}."
+        )
+    if len(batches) == len(labels):
+        raise ValueError(
+            "`batches` must leave at least one batch unaffected; "
+            "shifting every batch is a global shift that does not "
+            "act as a batch effect."
+        )
+
+
+def _validate_args(
+    adata,
+    batch_key,
+    batches,
+    flavor,
+    n_peptides,
+    log2fc,
+    fraction,
+    log_space,
+    random_state,
+):
+    """Validate arguments and return the batch labels."""
+    _, level = check_proteodata(adata)
+    if level != "peptide":
+        raise ValueError(
+            "apply_peptide_batch_effects requires peptide-level data, "
+            f"got {level}-level data."
+        )
+    if flavor not in FLAVORS:
+        raise ValueError(f"`flavor` must be one of {FLAVORS}, got {flavor!r}.")
+    if batch_key not in adata.obs.columns:
+        raise KeyError(f"batch_key '{batch_key}' not found in .obs.")
+    labels = _batch_labels(adata, batch_key)
+    if len(labels) < 2:
+        raise ValueError(
+            f".obs['{batch_key}'] must contain at least 2 batches, "
+            f"found {len(labels)}."
+        )
+    _validate_batches(batches, labels, batch_key)
+
+    if not _is_int(n_peptides) or n_peptides < 1:
+        raise ValueError(
+            f"`n_peptides` must be a positive int, got {n_peptides!r}."
+        )
+    if (
+        not isinstance(log2fc, Real)
+        or isinstance(log2fc, bool)
+        or not np.isfinite(log2fc)
+    ):
+        raise ValueError(f"`log2fc` must be a finite number, got {log2fc!r}.")
+    if (
+        not isinstance(fraction, Real)
+        or isinstance(fraction, bool)
+        or not 0 < fraction <= 1
+    ):
+        raise ValueError(f"`fraction` must be in (0, 1], got {fraction!r}.")
+    if not isinstance(log_space, bool):
+        raise TypeError(
+            f"`log_space` must be a bool, got {type(log_space).__name__}."
+        )
+    if random_state is not None and (
+        not _is_int(random_state) or random_state < 0
+    ):
+        raise ValueError(
+            "`random_state` must be None or a non-negative int, "
+            f"got {random_state!r}."
+        )
+    return labels
+
+
+def _draw_peptides(rng, proteins, n_peptides, fraction):
+    """Draw affected proteins and `n_peptides` peptides of each.
+
+    Returns the positional var indices of the drawn peptides and the
+    number of drawn proteins.
+    """
+    protein_indices = proteins.groupby(proteins).indices
+    # A protein must keep at least one unshifted peptide, else the
+    # shift moves the whole protein
+    eligible = sorted(
+        p for p, idx in protein_indices.items() if len(idx) > n_peptides
+    )
+    n_affected = int(round(fraction * len(eligible)))
+    if n_affected == 0:
+        return np.array([], dtype=int), 0
+    affected = rng.choice(eligible, n_affected, replace=False)
+    var_idx = np.concatenate(
+        [
+            rng.choice(protein_indices[p], n_peptides, replace=False)
+            for p in affected
+        ]
+    )
+    return np.sort(var_idx), n_affected
+
+
+def apply_peptide_batch_effects(
+    adata: ad.AnnData,
+    batch_key: str,
+    batches: int | Sequence,
+    *,
+    flavor: str = "subset_shift",
+    n_peptides: int = 2,
+    log2fc: float = -1.0,
+    fraction: float = 1.0,
+    log_space: bool = False,
+    random_state: int | None = None,
+    inplace: bool = True,
+    verbose: bool = False,
+) -> ad.AnnData | None:
+    r"""
+    Simulate peptide-specific batch effects in peptide-level data.
+
+    With ``flavor="subset_shift"``, ``n_peptides`` peptides are drawn
+    per protein and their intensities are changed by the same fold
+    change ``2 ** log2fc`` in every sample of the affected batches.
+    The drawn peptides then move together across batches, apart from
+    the other peptides of their protein, mimicking batch effects that
+    act on single peptides (ionisation, digestion, retention-time
+    drift) rather than on whole proteins.
+
+    With :math:`B` the samples of the affected batches and :math:`P`
+    the shifted peptides, the ``"subset_shift"`` flavor computes, in
+    linear space (``log_space=False``) and in log2 space
+    (``log_space=True``) respectively:
+
+    .. math::
+
+        X'_{s,p} = X_{s,p} \cdot 2^{\mathrm{log2fc}}
+        \qquad
+        X'_{s,p} = X_{s,p} + \mathrm{log2fc}
+        \qquad
+        \text{for } s \in B,\ p \in P
+
+    All other values are unchanged.
+
+    Parameters
+    ----------
+    adata : AnnData
+        Peptide-level :class:`~anndata.AnnData`.
+    batch_key : str
+        ``.obs`` column defining the batches. Observations with a
+        missing batch label are never affected.
+    batches : int | Sequence
+        Batches receiving the effect. A list gives the batch labels
+        (values of ``.obs[batch_key]``). An int draws that many batches
+        at random. At least one batch must stay unaffected.
+    flavor : str
+        Kind of batch effect. Only ``"subset_shift"`` is available.
+    n_peptides : int
+        Number of peptides shifted per affected protein. Only proteins
+        with more than ``n_peptides`` peptides are eligible, so each
+        keeps at least one unshifted peptide. The same peptides are
+        shifted in all affected batches.
+    log2fc : float
+        Log2 fold change applied to the drawn peptides in the affected
+        batches, e.g. ``-1`` halves their intensities.
+    fraction : float
+        Fraction of eligible proteins, randomly selected, that receive
+        the effect.
+    log_space : bool
+        If False, ``.X`` holds linear intensities and is multiplied by
+        ``2 ** log2fc``. If True, ``.X`` is assumed to be
+        log2-transformed and ``log2fc`` is added. Data in another log
+        base must be converted to log2 first. Missing values stay
+        missing. In linear space zeros stay zero; in log space every
+        finite value is shifted.
+    random_state : int | None
+        Seed for drawing batches, proteins and peptides. If ``None``, a
+        random seed is drawn and stored in
+        ``.uns['batch_effect_simulation']``.
+    inplace : bool
+        If True, modify ``adata`` in place; else return a new AnnData.
+    verbose : bool
+        If True, print the affected batches, the number of affected
+        proteins and peptides, and where the annotations are stored.
+
+    Returns
+    -------
+    AnnData | None
+        ``None`` if ``inplace=True``, else a new AnnData. ``.X`` is
+        dense. Adds:
+
+        - ``.var['batch_effect']``: True for the shifted peptides
+        - ``.uns['batch_effect_simulation']``: ``flavor``,
+          ``batch_key``, the affected ``batches``, ``n_peptides``,
+          ``log2fc``, ``fraction``, ``log_space``, ``random_state``,
+          ``n_proteins_affected`` and ``n_peptides_affected``
+
+        A second call overwrites both; shifts already applied to
+        ``.X`` are kept.
+
+    Raises
+    ------
+    ValueError
+        If ``adata`` is not peptide-level, ``.obs[batch_key]`` has
+        fewer than 2 batches, ``batches`` names unknown batches or all
+        of them, or a numeric argument is out of range.
+    KeyError
+        If ``batch_key`` is not in ``.obs``.
+    TypeError
+        If ``batches`` is neither an int nor a list of labels, or
+        ``log_space`` is not a bool.
+
+    Warns
+    -----
+    UserWarning
+        If ``.X`` is sparse (it is densified), or if ``log_space``
+        does not match whether ``.X`` looks log-transformed.
+
+    Examples
+    --------
+    Shift specific batches of log2-transformed data:
+
+    >>> import proteopy as pr
+    >>> pr.utils.apply_peptide_batch_effects(
+    ...     adata,
+    ...     batch_key="batch",
+    ...     batches=["batch_1"],
+    ...     log2fc=-3,
+    ...     log_space=True,
+    ... )
+    """
+    labels = _validate_args(
+        adata,
+        batch_key,
+        batches,
+        flavor,
+        n_peptides,
+        log2fc,
+        fraction,
+        log_space,
+        random_state,
+    )
+
+    Xsrc = adata.X
+    if sparse.issparse(Xsrc):
+        warnings.warn(
+            "Sparse input is being densified for computation.",
+            UserWarning,
+            stacklevel=2,
+        )
+        X = Xsrc.toarray().astype(float)
+    else:
+        X = np.array(Xsrc, dtype=float)
+
+    looks_log = is_log_transformed(adata)[0]
+    if looks_log != log_space:
+        warnings.warn(
+            f"`log_space={log_space}` but `.X` looks "
+            f"{'log-transformed' if looks_log else 'linear'}, so the "
+            f"effect is applied "
+            f"{'additively' if log_space else 'multiplicatively'}. Check "
+            "that `log_space` matches the scale of `.X`.",
+            UserWarning,
+            stacklevel=2,
+        )
+
+    if random_state is None:
+        random_state = int(
+            np.random.SeedSequence().generate_state(1, np.uint32)[0]
+        )
+    rng = np.random.default_rng(random_state)
+
+    # -- Draw batches, proteins and peptides
+    if _is_int(batches):
+        idx = rng.choice(len(labels), batches, replace=False)
+        affected_batches = [labels[i] for i in sorted(idx)]
+    else:
+        affected_batches = list(batches)
+    # Plain Python scalars keep .uns writable to h5ad
+    affected_batches = [
+        b.item() if isinstance(b, np.generic) else b for b in affected_batches
+    ]
+    var_idx, n_proteins = _draw_peptides(
+        rng,
+        adata.var["protein_id"].astype(str).reset_index(drop=True),
+        n_peptides,
+        fraction,
+    )
+    obs_mask = adata.obs[batch_key].isin(affected_batches).to_numpy()
+
+    # -- Shift
+    block = np.ix_(obs_mask, var_idx)
+    if log_space:
+        X[block] += log2fc
+    else:
+        X[block] *= 2.0**log2fc
+
+    var_flag = np.zeros(adata.n_vars, dtype=bool)
+    var_flag[var_idx] = True
+    params = {
+        "flavor": flavor,
+        "batch_key": batch_key,
+        "batches": affected_batches,
+        "n_peptides": int(n_peptides),
+        "log2fc": float(log2fc),
+        "fraction": float(fraction),
+        "log_space": log_space,
+        "random_state": int(random_state),
+        "n_proteins_affected": int(n_proteins),
+        "n_peptides_affected": int(len(var_idx)),
+    }
+
+    if verbose:
+        print(
+            f"Shifted {len(var_idx)} peptides of {n_proteins} proteins "
+            f"by log2fc={log2fc} in batches {affected_batches} of "
+            f".obs['{batch_key}'] ({int(obs_mask.sum())} samples, "
+            f"random_state={random_state}).\n"
+            f"Flagged shifted peptides in .var['{BATCH_EFFECT_VAR_KEY}'] "
+            f"and stored parameters in .uns['{BATCH_EFFECT_UNS_KEY}']."
+        )
+
+    out = adata if inplace else adata.copy()
+    out.X = X
+    out.var[BATCH_EFFECT_VAR_KEY] = var_flag
+    out.uns[BATCH_EFFECT_UNS_KEY] = params
+    check_proteodata(out)
+    return None if inplace else out
diff --git a/tests/utils/test_batch_effect_simulation.py b/tests/utils/test_batch_effect_simulation.py
new file mode 100644
index 0000000..2f5a6bc
--- /dev/null
+++ b/tests/utils/test_batch_effect_simulation.py
@@ -0,0 +1,317 @@
+"""Tests for ``pr.utils.apply_peptide_batch_effects``."""
+
+import warnings
+
+import anndata as ad
+import numpy as np
+import pandas as pd
+import pytest
+from scipy import sparse
+
+import proteopy as pr
+
+N_PER_BATCH = 6
+BATCHES = ["b1", "b2", "b3"]
+
+
+# --------------------------------------------------------------------
+# Fixtures and helpers
+# --------------------------------------------------------------------
+
+
+def make_peptide_adata(peptides_per_protein, seed=0):
+    """Raw-scale peptide-level data with three batches."""
+    rng = np.random.default_rng(seed)
+    var_rows = []
+    for p, n_peptides in enumerate(peptides_per_protein):
+        for i in range(n_peptides):
+            var_rows.append(
+                {
+                    "peptide_id": f"PROT{p:03d}_PEP{i}",
+                    "protein_id": f"PROT{p:03d}",
+                }
+            )
+    var = pd.DataFrame(var_rows)
+    var.index = var["peptide_id"].to_numpy()
+    n_obs = N_PER_BATCH * len(BATCHES)
+    obs = pd.DataFrame(
+        {
+            "sample_id": [f"s{i}" for i in range(n_obs)],
+            "batch": np.repeat(BATCHES, N_PER_BATCH),
+        }
+    )
+    obs.index = obs["sample_id"].to_numpy()
+    X = 2.0 ** rng.normal(20, 1, (n_obs, len(var)))
+    return ad.AnnData(X=X, obs=obs, var=var)
+
+
+@pytest.fixture
+def adata():
+    # 10 proteins with 3-5 peptides, 2 with only 2 peptides
+    return make_peptide_adata([3, 4, 5, 3, 4, 5, 3, 4, 5, 3, 2, 2])
+
+
+def apply(adata, **kwargs):
+    kwargs = {"batch_key": "batch", "batches": ["b2"], **kwargs}
+    return pr.utils.apply_peptide_batch_effects(adata, inplace=False, **kwargs)
+
+
+def ratio(out, adata):
+    return out.X / adata.X
+
+
+# --------------------------------------------------------------------
+# Shift
+# --------------------------------------------------------------------
+
+
+def test_shifts_only_flagged_peptides_in_affected_batch(adata):
+    out = apply(adata, log2fc=-1.0, random_state=0)
+    flagged = out.var["batch_effect"].to_numpy()
+    in_batch = (out.obs["batch"] == "b2").to_numpy()
+    r = ratio(out, adata)
+
+    np.testing.assert_allclose(r[np.ix_(in_batch, flagged)], 0.5)
+    np.testing.assert_allclose(r[np.ix_(~in_batch, flagged)], 1.0)
+    np.testing.assert_allclose(r[:, ~flagged], 1.0)
+
+
+def test_log2fc_values(adata):
+    expected = {-1.0: 0.5, 1.0: 2.0, 2.0: 4.0, 0.0: 1.0}
+    for log2fc, factor in expected.items():
+        out = apply(adata, log2fc=log2fc, random_state=0)
+        flagged = out.var["batch_effect"].to_numpy()
+        in_batch = (out.obs["batch"] == "b2").to_numpy()
+        np.testing.assert_allclose(
+            ratio(out, adata)[np.ix_(in_batch, flagged)], factor
+        )
+
+
+def test_log_space_adds_log2fc(adata):
+    logged = adata.copy()
+    logged.X = np.log2(adata.X)
+    for log2fc in (-1.0, 0.5, 2.0):
+        out = apply(logged, log2fc=log2fc, log_space=True, random_state=0)
+        flagged = out.var["batch_effect"].to_numpy()
+        in_batch = (out.obs["batch"] == "b2").to_numpy()
+        diff = out.X - logged.X
+        np.testing.assert_allclose(diff[np.ix_(in_batch, flagged)], log2fc)
+        np.testing.assert_allclose(diff[np.ix_(~in_batch, flagged)], 0.0)
+        np.testing.assert_allclose(diff[:, ~flagged], 0.0)
+
+
+def test_log_space_matches_linear_space(adata):
+    linear = apply(adata, log2fc=-1.5, random_state=0)
+    logged = adata.copy()
+    logged.X = np.log2(adata.X)
+    out = apply(logged, log2fc=-1.5, log_space=True, random_state=0)
+    np.testing.assert_allclose(out.X, np.log2(linear.X))
+
+
+def test_n_peptides_per_affected_protein(adata):
+    out = apply(adata, n_peptides=2, random_state=0)
+    per_protein = out.var.groupby("protein_id")["batch_effect"].sum()
+    # Proteins with exactly 2 peptides are not eligible
+    assert (per_protein[["PROT010", "PROT011"]] == 0).all()
+    assert (per_protein.drop(["PROT010", "PROT011"]) == 2).all()
+
+
+def test_fraction_limits_affected_proteins(adata):
+    out = apply(adata, fraction=0.5, random_state=0)
+    per_protein = out.var.groupby("protein_id")["batch_effect"].sum()
+    assert (per_protein > 0).sum() == 5
+    assert out.uns["batch_effect_simulation"]["n_proteins_affected"] == 5
+
+
+def test_preserves_within_batch_correlation(adata):
+    out = apply(adata, log2fc=-2.0, random_state=0)
+    for b in BATCHES:
+        mask = (out.obs["batch"] == b).to_numpy()
+        np.testing.assert_allclose(
+            np.corrcoef(out.X[mask].T),
+            np.corrcoef(adata.X[mask].T),
+        )
+
+
+def test_nan_and_zero_stay(adata):
+    adata.X[0, 0] = np.nan
+    adata.X[N_PER_BATCH, :] = 0.0
+    out = apply(adata, random_state=0)
+    assert np.isnan(out.X[0, 0])
+    assert (out.X[N_PER_BATCH, :] == 0).all()
+
+
+def test_missing_batch_label_not_affected(adata):
+    adata.obs["batch"] = adata.obs["batch"].astype(object)
+    adata.obs.iloc[N_PER_BATCH, adata.obs.columns.get_loc("batch")] = None
+    out = apply(adata, random_state=0)
+    np.testing.assert_allclose(out.X[N_PER_BATCH], adata.X[N_PER_BATCH])
+
+
+# --------------------------------------------------------------------
+# Batch selection
+# --------------------------------------------------------------------
+
+
+def test_int_batches_draws_that_many(adata):
+    out = apply(adata, batches=2, random_state=0)
+    affected = out.uns["batch_effect_simulation"]["batches"]
+    assert len(affected) == 2
+    assert set(affected) <= set(BATCHES)
+    changed = ~np.isclose(ratio(out, adata), 1.0).all(axis=1)
+    assert set(out.obs.loc[changed, "batch"]) == set(affected)
+
+
+def test_numeric_batch_labels(adata):
+    adata.obs["batch"] = np.repeat([1, 2, 3], N_PER_BATCH)
+    out = apply(adata, batches=[3], random_state=0)
+    assert out.uns["batch_effect_simulation"]["batches"] == [3]
+    assert isinstance(out.uns["batch_effect_simulation"]["batches"][0], int)
+
+
+def test_categorical_batch_key(adata):
+    adata.obs["batch"] = pd.Categorical(adata.obs["batch"])
+    out = apply(adata, batches=["b3"], random_state=0)
+    in_batch = (out.obs["batch"] == "b3").to_numpy()
+    flagged = out.var["batch_effect"].to_numpy()
+    np.testing.assert_allclose(
+        ratio(out, adata)[np.ix_(in_batch, flagged)], 0.5
+    )
+
+
+# --------------------------------------------------------------------
+# Reproducibility, annotations and modes
+# --------------------------------------------------------------------
+
+
+def test_same_seed_same_result(adata):
+    a = apply(adata, batches=1, random_state=3)
+    b = apply(adata, batches=1, random_state=3)
+    np.testing.assert_array_equal(a.X, b.X)
+    assert a.uns["batch_effect_simulation"] == b.uns["batch_effect_simulation"]
+
+
+def test_none_seed_is_stored_and_reproduces(adata):
+    a = apply(adata, batches=1)
+    seed = a.uns["batch_effect_simulation"]["random_state"]
+    b = apply(adata, batches=1, random_state=seed)
+    np.testing.assert_array_equal(a.X, b.X)
+
+
+def test_uns_records_parameters(adata):
+    out = apply(adata, n_peptides=1, log2fc=0.5, random_state=1)
+    params = out.uns["batch_effect_simulation"]
+    assert params["flavor"] == "subset_shift"
+    assert params["batch_key"] == "batch"
+    assert params["batches"] == ["b2"]
+    assert params["n_peptides"] == 1
+    assert params["log2fc"] == 0.5
+    assert params["n_peptides_affected"] == out.var["batch_effect"].sum()
+
+
+def test_inplace_modifies_and_returns_none(adata):
+    X0 = adata.X.copy()
+    assert (
+        pr.utils.apply_peptide_batch_effects(
+            adata, batch_key="batch", batches=["b1"], random_state=0
+        )
+        is None
+    )
+    assert "batch_effect" in adata.var
+    assert not np.allclose(adata.X, X0)
+
+
+def test_not_inplace_leaves_input(adata):
+    X0 = adata.X.copy()
+    apply(adata, random_state=0)
+    np.testing.assert_array_equal(adata.X, X0)
+    assert "batch_effect" not in adata.var
+
+
+def test_sparse_input_warns_and_matches_dense(adata):
+    dense = apply(adata, random_state=0)
+    sp = adata.copy()
+    sp.X = sparse.csr_matrix(sp.X)
+    with pytest.warns(UserWarning, match="densified"):
+        out = apply(sp, random_state=0)
+    assert not sparse.issparse(out.X)
+    np.testing.assert_allclose(out.X, dense.X)
+
+
+def test_warns_when_log_data_treated_as_linear(adata):
+    adata.X = np.log2(adata.X)
+    with pytest.warns(UserWarning, match="looks log-transformed"):
+        apply(adata, random_state=0)
+
+
+def test_warns_when_linear_data_treated_as_log(adata):
+    with pytest.warns(UserWarning, match="looks linear"):
+        apply(adata, log_space=True, random_state=0)
+
+
+def test_no_warning_for_raw_data(adata):
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        apply(adata, random_state=0)
+
+
+def test_works_on_simulated_data(adata):
+    sim = next(
+        pr.utils.proteoform_group_simulator(
+            adata, n_mixed=3, n_unmixed=3, random_state=0
+        )
+    )
+    pr.utils.apply_peptide_batch_effects(
+        sim, batch_key="batch", batches=1, random_state=0
+    )
+    assert sim.var["batch_effect"].any()
+
+
+# --------------------------------------------------------------------
+# Validation
+# --------------------------------------------------------------------
+
+
+@pytest.mark.parametrize(
+    "kwargs, error",
+    [
+        ({"flavor": "peptide_mix"}, ValueError),
+        ({"batch_key": "missing"}, KeyError),
+        ({"batches": 0}, ValueError),
+        ({"batches": 3}, ValueError),
+        ({"batches": BATCHES}, ValueError),
+        ({"batches": []}, ValueError),
+        ({"batches": ["b1", "b1"]}, ValueError),
+        ({"batches": ["b9"]}, ValueError),
+        ({"batches": "b1"}, TypeError),
+        ({"n_peptides": 0}, ValueError),
+        ({"log2fc": np.inf}, ValueError),
+        ({"fraction": 0}, ValueError),
+        ({"fraction": 1.5}, ValueError),
+        ({"log_space": 2}, TypeError),
+        ({"log_space": None}, TypeError),
+        ({"random_state": -1}, ValueError),
+    ],
+)
+def test_invalid_arguments(adata, kwargs, error):
+    with pytest.raises(error):
+        apply(adata, **kwargs)
+
+
+def test_single_batch_raises(adata):
+    adata.obs["batch"] = "b1"
+    with pytest.raises(ValueError, match="at least 2 batches"):
+        apply(adata, batches=["b1"])
+
+
+def test_protein_level_raises():
+    var = pd.DataFrame({"protein_id": ["P1", "P2"]}, index=["P1", "P2"])
+    obs = pd.DataFrame(
+        {"sample_id": ["s1", "s2"], "batch": ["b1", "b2"]},
+        index=["s1", "s2"],
+    )
+    prot = ad.AnnData(X=np.ones((2, 2)) * 1e6, obs=obs, var=var)
+    with pytest.raises(ValueError, match="peptide-level"):
+        pr.utils.apply_peptide_batch_effects(
+            prot, batch_key="batch", batches=["b1"]
+        )