diff --git a/benchmarks/benchmarks/preprocessing_log.py b/benchmarks/benchmarks/preprocessing_log.py
index 3243ac1ec1..c9f30c51e4 100644
--- a/benchmarks/benchmarks/preprocessing_log.py
+++ b/benchmarks/benchmarks/preprocessing_log.py
@@ -11,6 +11,7 @@
 import anndata as ad
 import numpy as np
 import zarr
+from anndata.acc import A
 
 import scanpy as sc
 
@@ -97,10 +98,10 @@ def peakmem_neighbors(self, *_) -> None:
         sc.pp.neighbors(self.adata)
 
     def time_bbknn(self, *_) -> None:
-        sc.pp.bbknn(self.adata, batch_key=self.adata.uns["batch_key"])
+        sc.pp.bbknn(self.adata, batches=A.obs[self.adata.uns["batch_key"]])
 
     def peakmem_bbknn(self, *_) -> None:
-        sc.pp.bbknn(self.adata, batch_key=self.adata.uns["batch_key"])
+        sc.pp.bbknn(self.adata, batches=A.obs[self.adata.uns["batch_key"]])
 
 
 class HVGSuite:  # noqa: D101
diff --git a/docs/conf.py b/docs/conf.py
index 5448ef2e61..f7769e2635 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -161,6 +161,7 @@
     matplotlib=("https://matplotlib.org/stable/", None),
     networkx=("https://networkx.org/documentation/stable/", None),
     numpy=("https://numpy.org/doc/stable/", None),
+    palantir=("https://palantir.readthedocs.io/en/stable/", None),
     pandas=("https://pandas.pydata.org/pandas-docs/stable/", None),
     pydeseq2=("https://pydeseq2.readthedocs.io/en/stable/", None),
     pynndescent=("https://pynndescent.readthedocs.io/en/latest/", None),
diff --git a/docs/external/index.md b/docs/external/index.md
index ba17865182..618514564d 100644
--- a/docs/external/index.md
+++ b/docs/external/index.md
@@ -5,8 +5,8 @@
 .. module:: scanpy.external
 ```
 
-```{warning}
-We are no longer accepting new tools into `scanpy.external`.
+```{deprecated} 1.13.0
+`scanpy.external` is deprecated and will be removed in a future release.
 For tools that integrate well with scanpy and anndata, or to submit your own, see the [scverse ecosystem](https://scverse.org/packages/#ecosystem).
 ```
 
diff --git a/docs/release-notes/3645.chore.md b/docs/release-notes/3645.chore.md
new file mode 100644
index 0000000000..28e6862948
--- /dev/null
+++ b/docs/release-notes/3645.chore.md
@@ -0,0 +1 @@
+Deprecate {mod}`scanpy.external` and each of its functions. Each function’s deprecation message names its replacement: an upstream package, a successor project, or the scanpy function it will move to. See also the [scverse ecosystem](https://scverse.org/packages/#ecosystem) {smaller}`P Angerer`
diff --git a/src/scanpy/__init__.py b/src/scanpy/__init__.py
index a7d06628c3..f1bb3d7e69 100644
--- a/src/scanpy/__init__.py
+++ b/src/scanpy/__init__.py
@@ -10,12 +10,13 @@
 
 from anndata import AnnData, concat
 
-from . import datasets, experimental, external, get, io, logging, metrics, queries
+from . import datasets, experimental, get, io, logging, metrics, queries
 from . import plotting as pl
 from . import preprocessing as pp
 from . import tools as tl
 from ._utils import annotate_doc_types
 from .neighbors import Neighbors
+from .plotting.legacy.mpl_settings import set_figure_params
 
 if TYPE_CHECKING:
     from typing import Any
@@ -28,7 +29,6 @@
     "concat",
     "datasets",
     "experimental",
-    "external",
     "get",
     "io",
     "logging",
@@ -42,11 +42,8 @@
 ]
 
 
-from .plotting.legacy.mpl_settings import set_figure_params
-
 annotate_doc_types(sys.modules[__name__], "scanpy")
 
-# has to be done at the end, after everything has been imported
 sys.modules.update({f"{__name__}.{m}": globals()[m] for m in ["tl", "pp", "pl"]})
 
 
@@ -62,6 +59,11 @@ def __getattr__(name: str) -> Any:
         warn(msg, FutureWarning)
         return version("scanpy")
 
+    if name == "external":
+        import scanpy.external
+
+        return scanpy.external
+
     if name in {"read", "read_10x_h5", "read_10x_mtx", "write"} or (
         name in io.__all__ and name in anndata.io.__all__
     ):
diff --git a/src/scanpy/external/pl.py b/src/scanpy/external/pl.py
index 368fa68564..9bf8e60fa2 100644
--- a/src/scanpy/external/pl.py
+++ b/src/scanpy/external/pl.py
@@ -38,6 +38,11 @@
 ]
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0", "Use :func:`scanpy.pl.embedding` with ``basis='phate'`` instead."
+    )
+)
 @doctest_needs("phate")
 @_wraps_plot_scatter
 @_doc_params(
@@ -77,6 +82,8 @@ def phate(adata: AnnData, **kwargs) -> list[Axes] | None:
     >>> adata = AnnData(data)
     >>> adata.obs["branches"] = branches
     >>> sce.tl.phate(adata, k=5, a=20, t=150)
+    FutureWarning: The function phate is deprecated and will be removed in the future. Use the `phate <https://github.com/KrishnaswamyLab/PHATE>`_ package directly.
+        sce.tl.phate(adata, k=5, a=20, t=150)
     >>> adata.obsm["X_phate"].shape
     (2000, 2)
     >>> sce.pl.phate(
@@ -84,11 +91,18 @@ def phate(adata: AnnData, **kwargs) -> list[Axes] | None:
     ...     color="branches",
     ...     color_map="tab20",
     ... )
+    FutureWarning: The function phate is deprecated and will be removed in the future. Use :func:`scanpy.pl.embedding` with ``basis='phate'`` instead.
+        sce.pl.phate(
 
     """
     return embedding(adata, "phate", **kwargs)
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0", "Use :func:`scanpy.pl.embedding` with ``basis='trimap'`` instead."
+    )
+)
 @_wraps_plot_scatter
 @_doc_params(
     adata_color_etc=doc_adata_color_etc,
@@ -114,6 +128,11 @@ def trimap(adata: AnnData, **kwargs) -> Axes | list[Axes] | None:
     return embedding(adata, "trimap", **kwargs)
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0", "Use :func:`scanpy.pl.embedding` with ``basis='harmony'`` instead."
+    )
+)
 @_wraps_plot_scatter
 @_doc_params(
     adata_color_etc=doc_adata_color_etc,
@@ -162,6 +181,7 @@ def harmony_timeseries(
     return axes
 
 
+@deprecated(Deprecation("1.13.0", "Use :func:`scanpy.pl.embedding` instead."))
 def sam(
     adata: AnnData,
     projection: str | np.ndarray = "X_umap",
@@ -256,6 +276,13 @@ def sam(
     return axes
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`palantir.plot.plot_gene_trends` from Wishbone’s successor "
+        "`palantir <https://github.com/dpeerlab/Palantir>`_ instead.",
+    )
+)
 @_doc_params(show_save_ax=doc_show_save_ax)
 def wishbone_marker_trajectory(  # noqa: PLR0913
     adata: AnnData,
diff --git a/src/scanpy/external/pp/_magic.py b/src/scanpy/external/pp/_magic.py
index e9b6130d55..ab1f2db934 100644
--- a/src/scanpy/external/pp/_magic.py
+++ b/src/scanpy/external/pp/_magic.py
@@ -6,6 +6,7 @@
 from typing import TYPE_CHECKING
 
 from packaging.version import Version
+from scverse_misc import Deprecation, deprecated
 
 from ... import logging as logg
 from ..._compat import pkg_version
@@ -23,6 +24,13 @@
 MIN_VERSION = "2.0"
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`palantir.utils.run_magic_imputation` from MAGIC’s successor "
+        "`palantir <https://github.com/dpeerlab/Palantir>`_ instead.",
+    )
+)
 @doctest_needs("magic")
 def magic(  # noqa: PLR0913
     adata: AnnData,
@@ -127,12 +135,18 @@ def magic(  # noqa: PLR0913
     >>> sc.pp.normalize_total(adata)
     >>> sc.pp.sqrt(adata)  # or sc.pp.log1p(adata)
     >>> adata_magic = sce.pp.magic(adata, name_list=["Mpo", "Klf1", "Ifitm1"], knn=5)
+    FutureWarning: The function magic is deprecated and will be removed in the future. Use :func:`palantir.utils.run_magic_imputation` from MAGIC’s successor `palantir <https://github.com/dpeerlab/Palantir>`_ instead.
+        adata_magic = sce.pp.magic(adata, name_list=["Mpo", "Klf1", "Ifitm1"], knn=5)
     >>> adata_magic.shape
     (2730, 3)
     >>> sce.pp.magic(adata, name_list="pca_only", knn=5)
+    FutureWarning: The function magic is deprecated and will be removed in the future. Use :func:`palantir.utils.run_magic_imputation` from MAGIC’s successor `palantir <https://github.com/dpeerlab/Palantir>`_ instead.
+        sce.pp.magic(adata, name_list="pca_only", knn=5)
     >>> adata.obsm["X_magic"].shape
     (2730, 100)
     >>> sce.pp.magic(adata, name_list="all_genes", knn=5)
+    FutureWarning: The function magic is deprecated and will be removed in the future. Use :func:`palantir.utils.run_magic_imputation` from MAGIC’s successor `palantir <https://github.com/dpeerlab/Palantir>`_ instead.
+        sce.pp.magic(adata, name_list="all_genes", knn=5)
     >>> adata.X.shape
     (2730, 3451)
 
diff --git a/src/scanpy/external/pp/_mnn_correct.py b/src/scanpy/external/pp/_mnn_correct.py
index 4400aded12..26278e536a 100644
--- a/src/scanpy/external/pp/_mnn_correct.py
+++ b/src/scanpy/external/pp/_mnn_correct.py
@@ -2,6 +2,8 @@
 
 from typing import TYPE_CHECKING
 
+from scverse_misc import Deprecation, deprecated
+
 from ..._settings import settings
 
 if TYPE_CHECKING:
@@ -13,6 +15,12 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "For batch correction, use :func:`scanpy.pp.bbknn` or :func:`scanpy.pp.combat` instead.",
+    )
+)
 def mnn_correct(  # noqa: PLR0913
     *datas: AnnData | np.ndarray,
     var_index: Collection[str] | None = None,
diff --git a/src/scanpy/external/pp/_scanorama_integrate.py b/src/scanpy/external/pp/_scanorama_integrate.py
index 24f35ef954..d6d37e2748 100644
--- a/src/scanpy/external/pp/_scanorama_integrate.py
+++ b/src/scanpy/external/pp/_scanorama_integrate.py
@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+from scverse_misc import Deprecation, deprecated
 
 from ..._utils._doctests import doctest_needs
 
@@ -12,6 +13,12 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use `snapatac2.pp.scanorama_integrate <https://scverse.org/SnapATAC2/api/_autosummary/snapatac2.pp.scanorama_integrate.html>`_ instead.",
+    )
+)
 @doctest_needs("scanorama")
 def scanorama_integrate(
     adata: AnnData,
@@ -97,6 +104,8 @@ def scanorama_integrate(
     ``adata.obsm`` containing the Scanorama embeddings.
 
     >>> sce.pp.scanorama_integrate(adata, "batch", verbose=1)
+    FutureWarning: The function scanorama_integrate is deprecated and will be removed in the future. Use `snapatac2.pp.scanorama_integrate <https://scverse.org/SnapATAC2/api/_autosummary/snapatac2.pp.scanorama_integrate.html>`_ instead.
+        sce.pp.scanorama_integrate(adata, "batch", verbose=1)
     Processing datasets a <=> b
     >>> "X_scanorama" in adata.obsm
     True
diff --git a/src/scanpy/external/tl/_harmony_timeseries.py b/src/scanpy/external/tl/_harmony_timeseries.py
index e59fd3503d..5a53045782 100644
--- a/src/scanpy/external/tl/_harmony_timeseries.py
+++ b/src/scanpy/external/tl/_harmony_timeseries.py
@@ -6,6 +6,7 @@
 
 import numpy as np
 import pandas as pd
+from scverse_misc import Deprecation, deprecated
 
 from ... import logging as logg
 from ..._utils._doctests import doctest_needs
@@ -14,6 +15,12 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use the `harmonyTS <https://github.com/dpeerlab/Harmony>`_ package directly.",
+    )
+)
 @doctest_needs("harmony")
 def harmony_timeseries(
     adata: AnnData,
@@ -118,10 +125,14 @@ def harmony_timeseries(
     Run harmony_timeseries
 
     >>> sce.tl.harmony_timeseries(adata, tp="time_points", n_components=500)
+    FutureWarning: The function harmony_timeseries is deprecated and will be removed in the future. Use the `harmonyTS <https://github.com/dpeerlab/Harmony>`_ package directly.
+        sce.tl.harmony_timeseries(adata, tp="time_points", n_components=500)
 
     Plot time points:
 
     >>> sce.pl.harmony_timeseries(adata)
+    FutureWarning: The function harmony_timeseries is deprecated and will be removed in the future. Use :func:`scanpy.pl.embedding` with ``basis='harmony'`` instead.
+        sce.pl.harmony_timeseries(adata)
 
     For further demonstration of Harmony visualizations please follow the notebook
     `Harmony_sample_notebook.ipynb
diff --git a/src/scanpy/external/tl/_palantir.py b/src/scanpy/external/tl/_palantir.py
index ad24229ff7..4fbf8bfeab 100644
--- a/src/scanpy/external/tl/_palantir.py
+++ b/src/scanpy/external/tl/_palantir.py
@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING
 
 import pandas as pd
+from scverse_misc import Deprecation, deprecated
 
 from ... import logging as logg
 from ..._utils._doctests import doctest_needs
@@ -13,6 +14,13 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`palantir.utils.run_diffusion_maps` and "
+        ":func:`palantir.utils.determine_multiscale_space` directly.",
+    )
+)
 @doctest_needs("palantir")
 def palantir(
     adata: AnnData,
@@ -135,6 +143,8 @@ def palantir(
     dimensional phenotypic manifold of the data.
 
     >>> sce.tl.palantir(adata, n_components=5, knn=30)
+    FutureWarning: The function palantir is deprecated and will be removed in the future. Use :func:`palantir.utils.run_diffusion_maps` and :func:`palantir.utils.determine_multiscale_space` directly.
+        sce.tl.palantir(adata, n_components=5, knn=30)
 
     if pre-computed distances are to be used,
 
@@ -145,6 +155,8 @@ def palantir(
     ...     use_adjacency_matrix=True,
     ...     distances_key="distances",
     ... )
+    FutureWarning: The function palantir is deprecated and will be removed in the future. Use :func:`palantir.utils.run_diffusion_maps` and :func:`palantir.utils.determine_multiscale_space` directly.
+        sce.tl.palantir(
 
     **Visualizing Palantir results**
 
@@ -183,6 +195,8 @@ def palantir(
     ...     ms_data="X_palantir_multiscale",
     ...     num_waypoints=500,
     ... )
+    FutureWarning: The function palantir_results is deprecated and will be removed in the future. Use :func:`palantir.core.run_palantir` directly.
+        pr_res = sce.tl.palantir_results(
 
     .. note::
        A `start_cell` must be defined for every data set. The start cell for
@@ -244,6 +258,12 @@ def palantir(
     return adata if copy else None
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`palantir.core.run_palantir` directly.",
+    )
+)
 def palantir_results(
     adata: AnnData,
     early_cell: str,
diff --git a/src/scanpy/external/tl/_phate.py b/src/scanpy/external/tl/_phate.py
index 54d12643cf..d50000a49e 100644
--- a/src/scanpy/external/tl/_phate.py
+++ b/src/scanpy/external/tl/_phate.py
@@ -4,6 +4,8 @@
 
 from typing import TYPE_CHECKING
 
+from scverse_misc import Deprecation, deprecated
+
 from ... import logging as logg
 from ..._settings import settings
 from ..._utils._doctests import doctest_needs
@@ -16,6 +18,12 @@
     from ..._utils.random import _LegacyRandom
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use the `phate <https://github.com/KrishnaswamyLab/PHATE>`_ package directly.",
+    )
+)
 @doctest_needs("phate")
 def phate(  # noqa: PLR0913
     adata: AnnData,
@@ -123,9 +131,13 @@ def phate(  # noqa: PLR0913
     (2000, 100)
     >>> adata = AnnData(tree_data)
     >>> sce.tl.phate(adata, k=5, a=20, t=150)
+    FutureWarning: The function phate is deprecated and will be removed in the future. Use the `phate <https://github.com/KrishnaswamyLab/PHATE>`_ package directly.
+        sce.tl.phate(adata, k=5, a=20, t=150)
     >>> adata.obsm["X_phate"].shape
     (2000, 2)
     >>> sce.pl.phate(adata)
+    FutureWarning: The function phate is deprecated and will be removed in the future. Use :func:`scanpy.pl.embedding` with ``basis='phate'`` instead.
+        sce.pl.phate(adata)
 
     """
     start = logg.info("computing PHATE")
diff --git a/src/scanpy/external/tl/_phenograph.py b/src/scanpy/external/tl/_phenograph.py
index 6a2b59f1d2..8e0bae7694 100644
--- a/src/scanpy/external/tl/_phenograph.py
+++ b/src/scanpy/external/tl/_phenograph.py
@@ -7,6 +7,7 @@
 import pandas as pd
 from anndata import AnnData
 from packaging.version import Version
+from scverse_misc import Deprecation, deprecated
 
 from ... import logging as logg
 from ..._compat import pkg_version
@@ -21,6 +22,14 @@
     from ...tools._leiden import MutableVertexPartition
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`scanpy.pp.neighbors` with ``method='jaccard'``, followed by "
+        ":func:`scanpy.tl.leiden` for clustering and "
+        "`cellmapper <https://github.com/quadbio/cellmapper>`_ for label transfer.",
+    )
+)
 @doctest_needs("phenograph")
 def phenograph(  # noqa: PLR0913
     data: AnnData | np.ndarray | SpBase,
@@ -161,14 +170,20 @@ def phenograph(  # noqa: PLR0913
     **Louvain** community detection
 
     >>> sce.tl.phenograph(adata, clustering_algo="louvain", k=30)
+    FutureWarning: The function phenograph is deprecated and will be removed in the future. Use :func:`scanpy.pp.neighbors` with ``method='jaccard'``, followed by :func:`scanpy.tl.leiden` for clustering and `cellmapper <https://github.com/quadbio/cellmapper>`_ for label transfer.
+        sce.tl.phenograph(adata, clustering_algo="louvain", k=30)
 
     **Leiden** community detection
 
     >>> sce.tl.phenograph(adata, clustering_algo="leiden", k=30)
+    FutureWarning: The function phenograph is deprecated and will be removed in the future. Use :func:`scanpy.pp.neighbors` with ``method='jaccard'``, followed by :func:`scanpy.tl.leiden` for clustering and `cellmapper <https://github.com/quadbio/cellmapper>`_ for label transfer.
+        sce.tl.phenograph(adata, clustering_algo="leiden", k=30)
 
     Return only `Graph` object
 
     >>> sce.tl.phenograph(adata, clustering_algo=None, k=30)
+    FutureWarning: The function phenograph is deprecated and will be removed in the future. Use :func:`scanpy.pp.neighbors` with ``method='jaccard'``, followed by :func:`scanpy.tl.leiden` for clustering and `cellmapper <https://github.com/quadbio/cellmapper>`_ for label transfer.
+        sce.tl.phenograph(adata, clustering_algo=None, k=30)
 
     Now to show phenograph on tSNE (for example):
 
@@ -197,6 +212,8 @@ def phenograph(  # noqa: PLR0913
     >>> adata = AnnData(df)
     >>> sc.pp.pca(adata, n_comps=20)
     >>> sce.tl.phenograph(adata, clustering_algo="leiden", k=50)
+    FutureWarning: The function phenograph is deprecated and will be removed in the future. Use :func:`scanpy.pp.neighbors` with ``method='jaccard'``, followed by :func:`scanpy.tl.leiden` for clustering and `cellmapper <https://github.com/quadbio/cellmapper>`_ for label transfer.
+        sce.tl.phenograph(adata, clustering_algo="leiden", k=50)
     >>> sc.tl.tsne(adata, random_state=1)
     >>> sc.pl.tsne(
     ...     adata,
diff --git a/src/scanpy/external/tl/_pypairs.py b/src/scanpy/external/tl/_pypairs.py
index 6b7eabee8b..9c1032fa66 100644
--- a/src/scanpy/external/tl/_pypairs.py
+++ b/src/scanpy/external/tl/_pypairs.py
@@ -6,6 +6,7 @@
 from typing import TYPE_CHECKING
 
 from packaging.version import Version
+from scverse_misc import Deprecation, deprecated
 
 from ..._compat import pkg_version
 from ..._settings import settings
@@ -20,6 +21,12 @@
     type Genes = Collection[str | int | bool]
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "For cell cycle scoring, use :func:`scanpy.tl.score_genes_cell_cycle` instead.",
+    )
+)
 @doctest_needs("pypairs")
 def sandbag(
     adata: AnnData,
@@ -65,6 +72,8 @@ def sandbag(
     >>> from pypairs import datasets
     >>> adata = datasets.leng15()
     >>> marker_pairs = sandbag(adata, fraction=0.5)
+    FutureWarning: The function sandbag is deprecated and will be removed in the future. For cell cycle scoring, use :func:`scanpy.tl.score_genes_cell_cycle` instead.
+        marker_pairs = sandbag(adata, fraction=0.5)
 
     """
     _check_available()
@@ -86,6 +95,12 @@ def sandbag(
     )
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "For cell cycle scoring, use :func:`scanpy.tl.score_genes_cell_cycle` instead.",
+    )
+)
 def cyclone(
     adata: AnnData,
     marker_pairs: Mapping[str, Collection[tuple[str, str]]] | None = None,
diff --git a/src/scanpy/external/tl/_sam.py b/src/scanpy/external/tl/_sam.py
index ea2b8fc94e..19dc13e206 100644
--- a/src/scanpy/external/tl/_sam.py
+++ b/src/scanpy/external/tl/_sam.py
@@ -4,6 +4,8 @@
 
 from typing import TYPE_CHECKING
 
+from scverse_misc import Deprecation, deprecated
+
 from ... import logging as logg
 from ..._utils._doctests import doctest_needs
 
@@ -18,6 +20,12 @@
         SAM = type("SAM", (), dict(__module__="samalg"))
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use the `sc-sam <https://github.com/atarashansky/self-assembling-manifold>`_ package directly.",
+    )
+)
 @doctest_needs("samalg")
 def sam(  # noqa: PLR0913
     adata: AnnData,
@@ -153,6 +161,8 @@ def sam(  # noqa: PLR0913
     algorithm as follows:
 
     >>> sam_obj = sce.tl.sam(adata, inplace=True)
+    FutureWarning: The function sam is deprecated and will be removed in the future. Use the `sc-sam <https://github.com/atarashansky/self-assembling-manifold>`_ package directly.
+        sam_obj = sce.tl.sam(adata, inplace=True)
 
     The input AnnData object should contain unstandardized, non-negative
     expression values. Preferably, the data should be log-normalized and no
@@ -168,6 +178,8 @@ def sam(  # noqa: PLR0913
     To visualize the output, we can use:
 
     >>> sce.pl.sam(adata, projection="X_umap")
+    FutureWarning: The function sam is deprecated and will be removed in the future. Use :func:`scanpy.pl.embedding` instead.
+        sce.pl.sam(adata, projection="X_umap")
 
     `sce.pl.sam` accepts all keyword arguments used in the
     `matplotlib.pyplot.scatter` function.
diff --git a/src/scanpy/external/tl/_trimap.py b/src/scanpy/external/tl/_trimap.py
index 1239014672..d079f81223 100644
--- a/src/scanpy/external/tl/_trimap.py
+++ b/src/scanpy/external/tl/_trimap.py
@@ -4,6 +4,8 @@
 
 from typing import TYPE_CHECKING
 
+from scverse_misc import Deprecation, deprecated
+
 from ... import logging as logg
 from ..._compat import CSBase
 from ..._keys import _existing_preset_keys
@@ -16,6 +18,13 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use the `trimap <https://github.com/eamid/trimap>`_ package directly, or its "
+        "`JAX implementation <https://github.com/google-research/google-research/tree/master/trimap>`_.",
+    )
+)
 @doctest_needs("trimap")
 def trimap(  # noqa: PLR0913
     adata: AnnData,
@@ -87,7 +96,11 @@ def trimap(  # noqa: PLR0913
     >>> import scanpy.external as sce
     >>> pbmc = sc.datasets.pbmc68k_reduced()
     >>> pbmc = sce.tl.trimap(pbmc, copy=True)
+    FutureWarning: The function trimap is deprecated and will be removed in the future. Use the `trimap <https://github.com/eamid/trimap>`_ package directly, or its `JAX implementation <https://github.com/google-research/google-research/tree/master/trimap>`_.
+        pbmc = sce.tl.trimap(pbmc, copy=True)
     >>> sce.pl.trimap(pbmc, color=["bulk_labels"], s=10)
+    FutureWarning: The function trimap is deprecated and will be removed in the future. Use :func:`scanpy.pl.embedding` with ``basis='trimap'`` instead.
+        sce.pl.trimap(pbmc, color=["bulk_labels"], s=10)
 
     """
     try:
diff --git a/src/scanpy/external/tl/_wishbone.py b/src/scanpy/external/tl/_wishbone.py
index cad8f490e4..6d720397ca 100644
--- a/src/scanpy/external/tl/_wishbone.py
+++ b/src/scanpy/external/tl/_wishbone.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 import pandas as pd
+from scverse_misc import Deprecation, deprecated
 
 from ... import logging
 from ..._utils._doctests import doctest_needs
@@ -15,6 +16,13 @@
     from anndata import AnnData
 
 
+@deprecated(
+    Deprecation(
+        "1.13.0",
+        "Use :func:`palantir.core.run_palantir` from Wishbone’s successor "
+        "`palantir <https://github.com/dpeerlab/Palantir>`_ instead.",
+    )
+)
 @doctest_needs("wishbone")
 def wishbone(
     adata: AnnData,
@@ -86,12 +94,16 @@ def wishbone(
     ...     adata=adata, start_cell='ACAAGAGACTTATC-1',
     ...     components=[2, 3], num_waypoints=150,
     ... )
+    FutureWarning: The function wishbone is deprecated and will be removed in the future. Use :func:`palantir.core.run_palantir` from Wishbone’s successor `palantir <https://github.com/dpeerlab/Palantir>`_ instead.
+        sce.tl.wishbone(
 
     **Visualizing Wishbone results**
 
     >>> sc.pl.tsne(adata, color=['trajectory_wishbone', 'branch_wishbone'])
     >>> markers = ['C1QA', 'PSAP', 'CD79A', 'CD79B', 'CST3', 'LYZ', 'MALAT1']
     >>> sce.pl.wishbone_marker_trajectory(adata, markers, show=True)
+    FutureWarning: The function wishbone_marker_trajectory is deprecated and will be removed in the future. Use :func:`palantir.plot.plot_gene_trends` from Wishbone’s successor `palantir <https://github.com/dpeerlab/Palantir>`_ instead.
+        sce.pl.wishbone_marker_trajectory(adata, markers, show=True)
 
     For further demonstration of Wishbone methods and visualization please follow the
     notebooks in the package `Wishbone_for_single_cell_RNAseq.ipynb
diff --git a/tests/external/test_harmony_timeseries.py b/tests/external/test_harmony_timeseries.py
deleted file mode 100644
index 66e37795e4..0000000000
--- a/tests/external/test_harmony_timeseries.py
+++ /dev/null
@@ -1,35 +0,0 @@
-from __future__ import annotations
-
-from itertools import product
-
-from anndata import concat
-
-import scanpy as sc
-import scanpy.external as sce
-from testing.scanpy._helpers.data import pbmc3k
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.harmony]
-
-
-def test_load_timepoints_from_anndata_list():
-    adata_ref = pbmc3k()
-    start = [596, 615, 1682, 1663, 1409, 1432]
-    adata = concat(
-        [adata_ref[i : i + 1000] for i in start],
-        join="outer",
-        label="sample",
-        keys=[f"sa{i}_Rep{j}" for i, j in product((1, 2, 3), (1, 2))],
-        index_unique="-",
-    )
-    adata.obs["time_points"] = adata.obs["sample"].str.split("_", expand=True)[0]
-    adata.obs["time_points"] = adata.obs["time_points"].astype("category")
-    sc.pp.normalize_total(adata, target_sum=10000)
-    sc.pp.log1p(adata)
-    sc.pp.highly_variable_genes(adata, n_top_genes=1000, subset=True)
-
-    sce.tl.harmony_timeseries(adata=adata, tp="time_points", n_components=None)
-    assert all([
-        adata.obsp["harmony_aff"].shape[0],
-        adata.obsp["harmony_aff_aug"].shape[0],
-    ]), "harmony_timeseries augmented affinity matrix Error!"
diff --git a/tests/external/test_magic.py b/tests/external/test_magic.py
deleted file mode 100644
index 243ce3b968..0000000000
--- a/tests/external/test_magic.py
+++ /dev/null
@@ -1,56 +0,0 @@
-from __future__ import annotations
-
-import numpy as np
-from anndata import AnnData
-
-import scanpy as sc
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.magic]
-
-A_list = [
-    [0, 0, 7, 0, 0],
-    [8, 5, 0, 2, 0],
-    [6, 0, 0, 2, 5],
-    [0, 0, 0, 1, 0],
-    [8, 8, 2, 1, 0],
-    [0, 0, 0, 4, 5],
-]
-
-
-def test_magic_default():
-    a = np.array(A_list, dtype="float32")
-    adata = AnnData(a)
-    sc.external.pp.magic(adata, knn=1)
-    # check raw unchanged
-    np.testing.assert_array_equal(adata.raw.X, a)
-    # check .X changed
-    assert not np.all(a == adata.X)
-    # check .X shape unchanged
-    assert adata.X.shape == a.shape
-
-
-def test_magic_pca_only():
-    a = np.array(A_list, dtype="float32")
-    # pca only
-    adata = AnnData(a)
-    n_pca = 3
-    sc.external.pp.magic(adata, knn=1, name_list="pca_only", n_pca=n_pca)
-    # check raw unchanged
-    np.testing.assert_array_equal(adata.X, a)
-    # check .X shape consistent with n_pca
-    assert adata.obsm["X_magic"].shape == (a.shape[0], n_pca)
-
-
-def test_magic_copy():
-    a = np.array(A_list, dtype="float32")
-    adata = AnnData(a)
-    adata_copy = sc.external.pp.magic(adata, knn=1, copy=True)
-    # check adata unchanged
-    np.testing.assert_array_equal(adata.X, a)
-    # check copy raw unchanged
-    np.testing.assert_array_equal(adata_copy.raw.X, a)
-    # check .X changed
-    assert not np.all(a == adata_copy.X)
-    # check .X shape unchanged
-    assert adata_copy.X.shape == a.shape
diff --git a/tests/external/test_palantir.py b/tests/external/test_palantir.py
deleted file mode 100644
index b6b084be3a..0000000000
--- a/tests/external/test_palantir.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from __future__ import annotations
-
-import scanpy.external as sce
-from testing.scanpy._helpers.data import pbmc3k_processed
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.palantir]
-
-
-def test_palantir_core():
-    adata = pbmc3k_processed()
-
-    sce.tl.palantir(adata=adata, n_components=5, knn=30)
-    assert adata.layers["palantir_imp"].shape[0], "palantir_imp matrix Error!"
diff --git a/tests/external/test_phenograph.py b/tests/external/test_phenograph.py
deleted file mode 100644
index 6a6e17abec..0000000000
--- a/tests/external/test_phenograph.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from __future__ import annotations
-
-import numpy as np
-import pandas as pd
-from anndata import AnnData
-
-import scanpy as sc
-import scanpy.external as sce
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.phenograph]
-
-
-def test_phenograph():
-    rng = np.random.default_rng()
-    df = rng.random((1000, 40))
-    dframe = pd.DataFrame(df)
-    dframe.index, dframe.columns = (map(str, dframe.index), map(str, dframe.columns))
-    adata = AnnData(dframe)
-    sc.pp.pca(adata, n_comps=20)
-    sce.tl.phenograph(adata, clustering_algo="leiden", k=50)
-    assert adata.obs["pheno_leiden"].shape[0], "phenograph_Community Detection Error!"
diff --git a/tests/external/test_sam.py b/tests/external/test_sam.py
deleted file mode 100644
index 8f6b79b93f..0000000000
--- a/tests/external/test_sam.py
+++ /dev/null
@@ -1,23 +0,0 @@
-from __future__ import annotations
-
-import numpy as np
-
-import scanpy as sc
-import scanpy.external as sce
-from testing.scanpy._helpers.data import pbmc3k
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.samalg]
-
-
-def test_sam():
-    rng = np.random.default_rng()
-    adata_ref = pbmc3k()
-    ix = rng.choice(adata_ref.shape[0], size=200, replace=False)
-    adata = adata_ref[ix, :].copy()
-    sc.pp.normalize_total(adata, target_sum=10000)
-    sc.pp.log1p(adata)
-    sce.tl.sam(adata, inplace=True)
-    uns_keys = list(adata.uns.keys())
-    obsm_keys = list(adata.obsm.keys())
-    assert all(["sam" in uns_keys, "X_umap" in obsm_keys, "neighbors" in uns_keys])
diff --git a/tests/external/test_scanorama_integrate.py b/tests/external/test_scanorama_integrate.py
deleted file mode 100644
index df90368861..0000000000
--- a/tests/external/test_scanorama_integrate.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from __future__ import annotations
-
-import scanpy as sc
-import scanpy.external as sce
-from testing.scanpy._helpers.data import pbmc68k_reduced
-from testing.scanpy._pytest.marks import needs
-
-pytestmark = [needs.scanorama]
-
-
-def test_scanorama_integrate():
-    """Test that Scanorama integration works.
-
-    This is a very simple test that just checks to see if the Scanorama
-    integrate wrapper succesfully added a new field to ``adata.obsm``
-    and makes sure it has the same dimensions as the original PCA table.
-    """
-    adata = pbmc68k_reduced()
-    sc.pp.pca(adata)
-    adata.obs["batch"] = 350 * ["a"] + 350 * ["b"]
-    sce.pp.scanorama_integrate(adata, "batch", approx=False)
-    assert adata.obsm["X_scanorama"].shape == adata.obsm["X_pca"].shape
diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py
index 9b8505e6fa..6b2b3a69b8 100644
--- a/tests/test_package_structure.py
+++ b/tests/test_package_structure.py
@@ -59,6 +59,14 @@ def test_descend_classes_and_funcs():
     assert {p.values[0] for p in api_functions} == funcs
 
 
+@pytest.mark.parametrize(
+    ("f", "qualname"),
+    [p for p in api_functions if p.id.startswith("sc.external.")],
+)
+def test_external_funcs_are_deprecated(f, qualname) -> None:
+    assert getattr(f, "__deprecated__", None), f"{qualname} is not deprecated"
+
+
 @pytest.mark.filterwarnings("error::FutureWarning:.*Import anndata.*")
 def test_import_future_anndata_import_warning():
     import scanpy