diff --git a/pyproject.toml b/pyproject.toml
index d2d7858ff1..c08a34a4f0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -107,10 +107,7 @@ test-full = [
     "zappy",
     # additional tested algorithms
     "scanpy[louvain]",
-    "scanpy[magic]",
     "scanpy[skmisc]",
-    "scanpy[harmony]",
-    "scanpy[scanorama]",
     "scanpy[dask-ml]",
 ]
 doc = [
@@ -142,11 +139,7 @@ dev = [
 paga = [ "igraph" ]
 louvain = [ "igraph", "louvain>=0.6.0,!=0.6.2" ] # Louvain community detection
 leiden = [ "igraph>=0.10", "leidenalg>=0.9.0" ]  # Leiden community detection
-bbknn = [ "bbknn" ]                              # Batch balanced KNN (batch correction)
-magic = [ "magic-impute>=2.0" ]                  # MAGIC imputation method
 skmisc = [ "scikit-misc>=0.1.3" ]                # highly_variable_genes method 'seurat_v3'
-harmony = [ "harmonypy" ]                        # Harmony dataset integration
-scanorama = [ "scanorama" ]                      # Scanorama dataset integration
 scrublet = [ "scikit-image>=0.20" ]              # Doublet detection with automatic thresholds
 # Acceleration
 rapids = [ "cudf>=0.9", "cuml>=0.9", "cugraph>=0.9" ] # GPU accelerated calculation of neighbors
diff --git a/src/scanpy/__init__.py b/src/scanpy/__init__.py
index 9212318f6a..92dcbd29f7 100644
--- a/src/scanpy/__init__.py
+++ b/src/scanpy/__init__.py
@@ -45,7 +45,7 @@
     )
 from anndata import AnnData, concat
 
-from . import datasets, experimental, external, get, logging, metrics, queries
+from . import datasets, experimental, get, logging, metrics, queries
 from . import plotting as pl
 from . import preprocessing as pp
 from . import tools as tl
@@ -67,7 +67,6 @@
     "concat",
     "datasets",
     "experimental",
-    "external",
     "get",
     "logging",
     "metrics",
@@ -91,3 +90,12 @@
     "tl",
     "write",
 ]
+
+
+def __getattr__(name: str) -> object:
+    if name == "external":
+        from . import external
+
+        return external
+
+    raise AttributeError(name)
diff --git a/src/scanpy/external/__init__.py b/src/scanpy/external/__init__.py
index 2087fc2671..c841b328f0 100644
--- a/src/scanpy/external/__init__.py
+++ b/src/scanpy/external/__init__.py
@@ -3,11 +3,16 @@
 from __future__ import annotations
 
 import sys
+import warnings
 
 from .. import _utils
 from . import exporting, pl, pp, tl
 
+__all__: list[str] = ["exporting", "pl", "pp", "tl"]
+
 _utils.annotate_doc_types(sys.modules[__name__], "scanpy")
-del sys, _utils
 
-__all__ = ["exporting", "pl", "pp", "tl"]
+msg = "The `scanpy.external` module is deprecated and will be removed in a future version."
+warnings.warn(msg, DeprecationWarning, stacklevel=2)
+
+del msg, sys, _utils, warnings
diff --git a/tests/external/test_harmony_integrate.py b/tests/external/test_harmony_integrate.py
deleted file mode 100644
index 2844354a2f..0000000000
--- a/tests/external/test_harmony_integrate.py
+++ /dev/null
@@ -1,23 +0,0 @@
-from __future__ import annotations
-
-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.harmonypy]
-
-
-def test_harmony_integrate():
-    """Test that Harmony integrate works.
-
-    This is a very simple test that just checks to see if the Harmony
-    integrate wrapper succesfully added a new field to ``adata.obsm``
-    and makes sure it has the same dimensions as the original PCA table.
-    """
-    adata = pbmc3k()
-    sc.pp.recipe_zheng17(adata)
-    sc.pp.pca(adata)
-    adata.obs["batch"] = 1350 * ["a"] + 1350 * ["b"]
-    sce.pp.harmony_integrate(adata, "batch")
-    assert adata.obsm["X_pca_harmony"].shape == adata.obsm["X_pca"].shape
diff --git a/tests/external/test_harmony_timeseries.py b/tests/external/test_harmony_timeseries.py
deleted file mode 100644
index 3b51c1ecfd..0000000000
--- a/tests/external/test_harmony_timeseries.py
+++ /dev/null
@@ -1,33 +0,0 @@
-from __future__ import annotations
-
-from itertools import product
-
-from anndata import AnnData
-
-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 = AnnData.concatenate(
-        *(adata_ref[i : i + 1000] for i in start),
-        join="outer",
-        batch_key="sample",
-        batch_categories=[f"sa{i}_Rep{j}" for i, j in product((1, 2, 3), (1, 2))],
-    )
-    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_hashsolo.py b/tests/external/test_hashsolo.py
deleted file mode 100644
index e481a25e2a..0000000000
--- a/tests/external/test_hashsolo.py
+++ /dev/null
@@ -1,38 +0,0 @@
-from __future__ import annotations
-
-import numpy as np
-import pandas as pd
-from anndata import AnnData
-
-import scanpy.external as sce
-
-
-def test_cell_demultiplexing():
-    import random
-
-    from scipy import stats
-
-    random.seed(52)
-    signal = stats.poisson.rvs(1000, 1, 990)
-    doublet_signal = stats.poisson.rvs(1000, 1, 10)
-    x = np.reshape(stats.poisson.rvs(500, 1, 10000), (1000, 10))
-    for idx, signal_count in enumerate(signal):
-        col_pos = idx % 10
-        x[idx, col_pos] = signal_count
-
-    for idx, signal_count in enumerate(doublet_signal):
-        col_pos = (idx % 10) - 1
-        x[idx, col_pos] = signal_count
-
-    test_data = AnnData(np.random.randint(0, 100, size=x.shape), obs=pd.DataFrame(x))
-    sce.pp.hashsolo(test_data, test_data.obs.columns)
-
-    doublets = ["Doublet"] * 10
-    classes = np.repeat(np.arange(10), 98).reshape(98, 10, order="F").ravel().tolist()
-    negatives = ["Negative"] * 10
-    expected = pd.array(doublets + classes + negatives, dtype="string")
-    classification = test_data.obs["Classification"].array.astype("string")
-    # This is a bit flaky, so allow some mismatches:
-    if (expected != classification).sum() > 3:
-        # Compare lists for better error message
-        assert classification.tolist() == expected.tolist()
diff --git a/tests/external/test_magic.py b/tests/external/test_magic.py
deleted file mode 100644
index abf8c46fa0..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(adata.X == A)
-    # 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(adata_copy.X == A)
-    # 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 cee9211eb4..0000000000
--- a/tests/external/test_phenograph.py
+++ /dev/null
@@ -1,21 +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():
-    df = np.random.rand(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 b1b5b56f00..0000000000
--- a/tests/external/test_sam.py
+++ /dev/null
@@ -1,22 +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():
-    adata_ref = pbmc3k()
-    ix = np.random.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/external/test_wishbone.py b/tests/external/test_wishbone.py
deleted file mode 100644
index 001b7d7b4f..0000000000
--- a/tests/external/test_wishbone.py
+++ /dev/null
@@ -1,27 +0,0 @@
-from __future__ import annotations
-
-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.wishbone]
-
-
-def test_run_wishbone():
-    adata = pbmc3k()
-    sc.pp.normalize_per_cell(adata)
-    sc.pp.neighbors(adata, n_pcs=15, n_neighbors=10)
-    sc.pp.pca(adata)
-    sc.tl.tsne(adata=adata, n_pcs=5, perplexity=30)
-    sc.tl.diffmap(adata, n_comps=10)
-
-    sce.tl.wishbone(
-        adata=adata,
-        start_cell="ACAAGAGACTTATC-1",
-        components=[2, 3],
-        num_waypoints=150,
-    )
-    assert all(k in adata.obs for k in ["trajectory_wishbone", "branch_wishbone"]), (
-        "Run Wishbone Error!"
-    )
diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py
index 067c7fa1a9..f9007e7934 100644
--- a/tests/test_package_structure.py
+++ b/tests/test_package_structure.py
@@ -27,10 +27,6 @@
     "sc.tl",
     "sc.pl",
     "sc.experimental.pp",
-    "sc.external.pp",
-    "sc.external.tl",
-    "sc.external.pl",
-    "sc.external.exporting",
     "sc.get",
     "sc.logging",
     # "sc.neighbors",  # Not documented
@@ -54,11 +50,19 @@
 
 
 @pytest.mark.xfail(reason="TODO: unclear if we want this to totally match, let’s see")
+@pytest.mark.filterwarnings("ignore::DeprecationWarning:scanpy.external")
 def test_descend_classes_and_funcs():
     funcs = set(descend_classes_and_funcs(scanpy, "scanpy"))
     assert {p.values[0] for p in api_functions} == funcs
 
 
+def test_external_is_deprecated() -> None:
+    with pytest.warns(DeprecationWarning, match=r"scanpy.external"):
+        import scanpy.external
+
+        importlib.reload(scanpy.external)
+
+
 @pytest.mark.filterwarnings("error::FutureWarning:.*Import anndata.*")
 def test_import_future_anndata_import_warning():
     import scanpy