diff --git a/HISTORY.md b/HISTORY.md
index 8909f00..6309b3e 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -28,6 +28,11 @@ and this project adheres to
   that overlap in the protein sequence, keeping the most abundant
   member of each group. Reimplements CCprofiler's
   `summarizeAlternativePeptideSequences(topN = 1)`
+- **Tools** (`pr.tl`): `infer_proteoform_groups()` wraps the
+  currently recommended proteoform inference approach behind one call.
+  `method` defaults to `"copf"` (correlations, dendrograms, clusters,
+  scores), but the default may change in future releases; pass it
+  explicitly for reproducibility
 - **Tools** (`pr.tl`): `peptide_proximity()` tests whether the peptides
   of a proteoform cluster sit closer together in the protein sequence
   than a random grouping of the same size. A reimplementation of
diff --git a/docs/sphinx/source/api/tl.rst b/docs/sphinx/source/api/tl.rst
index d3168aa..5fe9400 100644
--- a/docs/sphinx/source/api/tl.rst
+++ b/docs/sphinx/source/api/tl.rst
@@ -35,6 +35,7 @@ detection of functional proteoform groups from peptide-level quantitative data
    :toctree: generated/
    :nosignatures:
 
+   proteopy.tl.infer_proteoform_groups
    proteopy.tl.pairwise_peptide_correlations
    proteopy.tl.peptide_dendograms_by_correlation
    proteopy.tl.peptide_clusters_from_dendograms
diff --git a/proteopy/tl/__init__.py b/proteopy/tl/__init__.py
index 1e245f0..6a65e0c 100644
--- a/proteopy/tl/__init__.py
+++ b/proteopy/tl/__init__.py
@@ -11,3 +11,4 @@
     hclustv_cluster_ann,
     hclustv_profiles,
 )
+from .proteoform_group_inference import infer_proteoform_groups
diff --git a/proteopy/tl/proteoform_group_inference.py b/proteopy/tl/proteoform_group_inference.py
new file mode 100644
index 0000000..b0e92ee
--- /dev/null
+++ b/proteopy/tl/proteoform_group_inference.py
@@ -0,0 +1,249 @@
+"""Proteoform group inference: a stable entry point over the methods.
+
+`infer_proteoform_groups` dispatches to the currently recommended
+proteoform inference approach. The building blocks remain available
+individually (e.g. the COPF steps in `proteopy.tl.copf`) for users who
+need finer control over each step.
+"""
+
+from numbers import Real
+
+import anndata as ad
+
+from proteopy.tl.copf import (
+    pairwise_peptide_correlations,
+    peptide_clusters_from_dendograms,
+    peptide_dendograms_by_correlation,
+    proteoform_scores,
+)
+from proteopy.utils.anndata import check_proteodata
+
+DEFAULT_METHOD = "copf"
+
+
+def _infer_copf(
+    adata: ad.AnnData,
+    *,
+    batch_key: str | None,
+    min_score: float,
+    min_adj_pval: float,
+    verbose: bool,
+) -> None:
+    """Run the COPF pipeline in place on ``adata``."""
+    if verbose:
+        batches = (
+            f"batches from .obs['{batch_key}']" if batch_key else "no batches"
+        )
+        print(f"Computing pairwise peptide correlations ({batches}).")
+    pairwise_peptide_correlations(adata, batch_key=batch_key)
+
+    if verbose:
+        print("Building peptide dendrograms from correlations.")
+    peptide_dendograms_by_correlation(
+        adata,
+        method="agglomerative-hierarchical-clustering",
+    )
+
+    if verbose:
+        print("Cutting dendrograms into two peptide clusters.")
+    peptide_clusters_from_dendograms(
+        adata,
+        n_clusters=2,
+        min_peptides_per_cluster=2,
+    )
+
+    if verbose:
+        print(
+            f"Scoring proteoforms (min_score={min_score}, "
+            f"min_adj_pval={min_adj_pval})."
+        )
+    proteoform_scores(
+        adata,
+        min_score=min_score,
+        min_pval_adj=min_adj_pval,
+    )
+
+    if verbose:
+        is_pf = adata.var.groupby("protein_id", observed=True)[
+            "is_proteoform"
+        ].first()
+        print(
+            f"{int((is_pf == 1).sum())} of {len(is_pf)} proteins called "
+            "as having proteoforms. Results stored in .var and .uns."
+        )
+
+
+METHODS = {
+    "copf": _infer_copf,
+}
+
+
+def _validate_args(
+    adata: ad.AnnData,
+    method: str,
+    batch_key: str | None,
+    min_score: float,
+    min_adj_pval: float,
+) -> None:
+    _, level = check_proteodata(adata)
+    if level != "peptide":
+        raise ValueError(
+            "infer_proteoform_groups requires peptide-level data, "
+            f"got {level}-level data."
+        )
+
+    if method not in METHODS:
+        raise ValueError(
+            f"Unknown method '{method}'. "
+            f"Available methods: {sorted(METHODS)}."
+        )
+
+    if batch_key is not None and batch_key not in adata.obs.columns:
+        raise KeyError(f"batch_key '{batch_key}' not found in .obs.")
+
+    for name, value in (
+        ("min_score", min_score),
+        ("min_adj_pval", min_adj_pval),
+    ):
+        if isinstance(value, bool) or not isinstance(value, Real):
+            raise TypeError(
+                f"{name} must be a number, got {type(value).__name__}."
+            )
+
+    if not 0 < min_adj_pval <= 1:
+        raise ValueError(
+            f"min_adj_pval must be in (0, 1], got {min_adj_pval}."
+        )
+
+
+def infer_proteoform_groups(
+    adata: ad.AnnData,
+    method: str = DEFAULT_METHOD,
+    *,
+    batch_key: str | None = None,
+    min_score: float = 0.1,
+    min_adj_pval: float = 0.1,
+    inplace: bool = True,
+    verbose: bool = False,
+) -> ad.AnnData | None:
+    """
+    Infer proteoform groups from peptide-level data.
+
+    A wrapper around the currently recommended proteoform inference
+    approach in ProteoPy.
+
+    .. warning::
+
+       This is a **dynamic** function. The default ``method`` points
+       to the approach the ProteoPy developers currently recommend and
+       may change in future releases as better approaches become
+       available, which would change the results of calls that rely on
+       the default. For reproducible analyses, always pass ``method``
+       explicitly, e.g. ``method="copf"``.
+
+    Available methods:
+
+    ``"copf"`` (current default)
+        COPF, COrrelation based functional ProteoForm assessment
+        [1]_. Equivalent to running::
+
+            pr.tl.pairwise_peptide_correlations(
+                adata, batch_key=batch_key,
+            )
+            pr.tl.peptide_dendograms_by_correlation(
+                adata, method="agglomerative-hierarchical-clustering",
+            )
+            pr.tl.peptide_clusters_from_dendograms(
+                adata, n_clusters=2, min_peptides_per_cluster=2,
+            )
+            pr.tl.proteoform_scores(
+                adata, min_score=min_score, min_pval_adj=min_adj_pval,
+            )
+
+        Call these functions individually for finer control over each
+        step.
+
+    Parameters
+    ----------
+    adata : AnnData
+        Peptide-level :class:`~anndata.AnnData` with ``.var``
+        columns ``peptide_id`` and ``protein_id``.
+    method : str
+        Proteoform inference method. One of ``"copf"``. The default
+        tracks the currently recommended method and may change; pass
+        it explicitly for reproducibility.
+    batch_key : str | None
+        ``.obs`` column defining batches. For ``"copf"``, peptide
+        correlations are computed per batch and pooled across batches
+        (see :func:`~proteopy.tl.pairwise_peptide_correlations`). If
+        ``None``, all observations are treated as a single batch.
+    min_score : float
+        Minimum proteoform score for a protein to be called as having
+        proteoforms (see :func:`~proteopy.tl.proteoform_scores`).
+    min_adj_pval : float
+        Maximum Benjamini-Hochberg adjusted p-value for a protein to
+        be called as having proteoforms.
+    inplace : bool
+        If True, modify ``adata`` in place; else return a modified
+        copy.
+    verbose : bool
+        If True, print status messages about each step.
+
+    Returns
+    -------
+    AnnData | None
+        ``None`` if ``inplace=True``, else a modified copy of
+        ``adata``. For ``"copf"``, the following are written:
+
+        - ``.uns['pairwise_peptide_correlations']`` (and
+          ``.uns['pairwise_peptide_correlations_by_batch']`` if
+          ``batch_key`` is given)
+        - ``.uns['dendograms']`` and ``.uns['clusters']``
+        - ``.var`` columns ``cluster_id``, ``proteoform_id``,
+          ``proteoform_score``, ``proteoform_score_z``,
+          ``proteoform_score_dz``, ``proteoform_score_pval``,
+          ``proteoform_score_pval_adj`` and ``is_proteoform``
+
+    Examples
+    --------
+    >>> import proteopy as pr
+    >>> pr.tl.infer_proteoform_groups(
+    ...     adata,
+    ...     method="copf",
+    ...     batch_key="multiplex_batch",
+    ...     min_score=0.1,
+    ...     min_adj_pval=0.1,
+    ... )
+    >>> adata.var.loc[adata.var["is_proteoform"] == 1, "protein_id"]
+
+    References
+    ----------
+    .. [1] Bludau I, Frank M, Dörig C, Cai Y, Heusel M, Rosenberger G,
+       Picotti P, Collins BC, Röst H, and Aebersold R. "Systematic
+       Detection of Functional Proteoform Groups from Bottom-Up
+       Proteomic Datasets." *Nature Communications*, 12:3810, 2021.
+       DOI: 10.1038/s41467-021-24030-x.
+    """
+    _validate_args(
+        adata,
+        method=method,
+        batch_key=batch_key,
+        min_score=min_score,
+        min_adj_pval=min_adj_pval,
+    )
+
+    if verbose:
+        print(f"Inferring proteoform groups with method:\n=== {method} ===")
+    target = adata if inplace else adata.copy()
+    METHODS[method](
+        target,
+        batch_key=batch_key,
+        min_score=min_score,
+        min_adj_pval=min_adj_pval,
+        verbose=verbose,
+    )
+    check_proteodata(target)
+
+    if inplace:
+        return None
+    return target
diff --git a/tests/tl/test_proteoform_group_inference.py b/tests/tl/test_proteoform_group_inference.py
new file mode 100644
index 0000000..a28a5c6
--- /dev/null
+++ b/tests/tl/test_proteoform_group_inference.py
@@ -0,0 +1,201 @@
+"""Tests for ``pr.tl.infer_proteoform_groups``.
+
+The function is a thin dispatcher over the currently recommended
+proteoform inference method. These tests check that it:
+
+- validates its input
+- selects the right method
+- gives exactly the same result as running that method's individual
+  ``pr.tl`` steps by hand (for ``"copf"``: correlations, dendrograms,
+  clusters, scores)
+- honours the shared ``inplace`` / ``verbose`` conventions
+"""
+
+import anndata as ad
+import numpy as np
+import pandas as pd
+import pytest
+
+import proteopy as pr
+
+N_OBS = 12
+
+
+# --------------------------------------------------------------------
+# Fixtures and helpers
+# --------------------------------------------------------------------
+
+
+@pytest.fixture
+def peptide_adata():
+    """Peptide-level data with one proteoform and one uniform protein.
+
+    PROT_A carries two groups of three peptides with opposite profiles;
+    PROT_B has six peptides sharing one profile. Peptide IDs are
+    sequence-like because `pairwise_peptide_correlations` splits peptide
+    pairs on underscores.
+    """
+    rng = np.random.default_rng(0)
+    base = np.linspace(1.0, 10.0, N_OBS)
+    profiles = {
+        "PROT_A": [base] * 3 + [base[::-1]] * 3,
+        "PROT_B": [base] * 6,
+    }
+
+    columns, var_rows = [], []
+    for protein, protein_profiles in profiles.items():
+        for i, profile in enumerate(protein_profiles):
+            columns.append(profile + rng.normal(0, 0.3, N_OBS) + 20)
+            var_rows.append(
+                {
+                    "peptide_id": f"{protein[-1]}PEPTIDE{'ACDEFG'[i]}K",
+                    "protein_id": protein,
+                }
+            )
+
+    var = pd.DataFrame(var_rows)
+    var.index = var["peptide_id"].to_numpy()
+    obs = pd.DataFrame(
+        {
+            "sample_id": [f"s{i}" for i in range(N_OBS)],
+            "batch": ["b1"] * (N_OBS // 2) + ["b2"] * (N_OBS // 2),
+        }
+    )
+    obs.index = obs["sample_id"].to_numpy()
+    return ad.AnnData(X=np.column_stack(columns), obs=obs, var=var)
+
+
+def run_copf_steps(adata, batch_key=None, min_score=0.1, min_adj_pval=0.1):
+    pr.tl.pairwise_peptide_correlations(adata, batch_key=batch_key)
+    pr.tl.peptide_dendograms_by_correlation(
+        adata, method="agglomerative-hierarchical-clustering"
+    )
+    pr.tl.peptide_clusters_from_dendograms(
+        adata, n_clusters=2, min_peptides_per_cluster=2
+    )
+    pr.tl.proteoform_scores(
+        adata, min_score=min_score, min_pval_adj=min_adj_pval
+    )
+
+
+# --------------------------------------------------------------------
+# Input validation
+# --------------------------------------------------------------------
+
+
+def test_protein_level_raises(peptide_adata):
+    protein = peptide_adata[:, :1].copy()
+    protein.var = pd.DataFrame({"protein_id": ["PROT_A"]}, index=["PROT_A"])
+    with pytest.raises(ValueError, match="peptide-level"):
+        pr.tl.infer_proteoform_groups(protein)
+
+
+def test_unknown_method_raises(peptide_adata):
+    with pytest.raises(ValueError, match="Unknown method"):
+        pr.tl.infer_proteoform_groups(peptide_adata, method="foo")
+
+
+def test_missing_batch_key_raises(peptide_adata):
+    with pytest.raises(KeyError, match="batch_key"):
+        pr.tl.infer_proteoform_groups(peptide_adata, batch_key="missing")
+
+
+@pytest.mark.parametrize("arg", ["min_score", "min_adj_pval"])
+def test_non_numeric_threshold_raises(peptide_adata, arg):
+    with pytest.raises(TypeError, match=arg):
+        pr.tl.infer_proteoform_groups(peptide_adata, **{arg: "0.1"})
+
+
+@pytest.mark.parametrize("min_adj_pval", [0, -0.1, 1.5])
+def test_min_adj_pval_out_of_range_raises(peptide_adata, min_adj_pval):
+    with pytest.raises(ValueError, match="min_adj_pval"):
+        pr.tl.infer_proteoform_groups(peptide_adata, min_adj_pval=min_adj_pval)
+
+
+# --------------------------------------------------------------------
+# Method dispatch
+# --------------------------------------------------------------------
+
+
+def test_default_method_is_copf(peptide_adata):
+    expected = peptide_adata.copy()
+    pr.tl.infer_proteoform_groups(expected, method="copf")
+
+    pr.tl.infer_proteoform_groups(peptide_adata)
+
+    pd.testing.assert_frame_equal(peptide_adata.var, expected.var)
+
+
+# --------------------------------------------------------------------
+# COPF integration: parity with the individual pr.tl steps
+# --------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("batch_key", [None, "batch"])
+def test_copf_matches_individual_steps(peptide_adata, batch_key):
+    expected = peptide_adata.copy()
+    run_copf_steps(expected, batch_key=batch_key)
+
+    pr.tl.infer_proteoform_groups(
+        peptide_adata, method="copf", batch_key=batch_key
+    )
+
+    pd.testing.assert_frame_equal(peptide_adata.var, expected.var)
+    pd.testing.assert_frame_equal(
+        peptide_adata.uns["pairwise_peptide_correlations"],
+        expected.uns["pairwise_peptide_correlations"],
+    )
+    assert ("pairwise_peptide_correlations_by_batch" in peptide_adata.uns) == (
+        batch_key is not None
+    )
+
+
+def test_thresholds_are_forwarded(peptide_adata):
+    expected = peptide_adata.copy()
+    run_copf_steps(expected, min_score=0.5, min_adj_pval=0.01)
+
+    pr.tl.infer_proteoform_groups(
+        peptide_adata, method="copf", min_score=0.5, min_adj_pval=0.01
+    )
+
+    pd.testing.assert_series_equal(
+        peptide_adata.var["is_proteoform"], expected.var["is_proteoform"]
+    )
+
+
+def test_detects_proteoform(peptide_adata):
+    pr.tl.infer_proteoform_groups(peptide_adata, method="copf")
+
+    is_pf = peptide_adata.var.groupby("protein_id")["is_proteoform"].first()
+    assert is_pf["PROT_A"] == 1
+    assert is_pf["PROT_B"] != 1
+
+
+# --------------------------------------------------------------------
+# Output: return value, in-place modification and verbosity
+# --------------------------------------------------------------------
+
+
+def test_inplace_true_returns_none(peptide_adata):
+    assert pr.tl.infer_proteoform_groups(peptide_adata) is None
+    assert "is_proteoform" in peptide_adata.var.columns
+
+
+def test_inplace_false_returns_copy(peptide_adata):
+    original_var = peptide_adata.var.copy()
+
+    out = pr.tl.infer_proteoform_groups(
+        peptide_adata, method="copf", inplace=False
+    )
+
+    assert isinstance(out, ad.AnnData)
+    assert "is_proteoform" in out.var.columns
+    pd.testing.assert_frame_equal(peptide_adata.var, original_var)
+    assert "pairwise_peptide_correlations" not in peptide_adata.uns
+
+
+def test_verbose_prints(peptide_adata, capsys):
+    pr.tl.infer_proteoform_groups(peptide_adata, verbose=True)
+    out = capsys.readouterr().out
+    assert "copf" in out
+    assert "called as having proteoforms" in out