diff --git a/src/squidpy/_constants/_constants.py b/src/squidpy/_constants/_constants.py index f88e41cb6..fb735130e 100644 --- a/src/squidpy/_constants/_constants.py +++ b/src/squidpy/_constants/_constants.py @@ -94,6 +94,7 @@ class Symbol(ModeEnum): class SpatialAutocorr(ModeEnum): MORAN = "moran" GEARY = "geary" + SPAGFT = "spagft" @unique diff --git a/src/squidpy/gr/_ppatterns.py b/src/squidpy/gr/_ppatterns.py index 292c75994..1e744b7ae 100644 --- a/src/squidpy/gr/_ppatterns.py +++ b/src/squidpy/gr/_ppatterns.py @@ -34,6 +34,7 @@ __all__ = ["spatial_autocorr", "co_occurrence"] +from squidpy.gr._spagft import _spagft it = nt.int32 ft = nt.float32 @@ -188,6 +189,11 @@ def extract_obsm(adata: AnnData, ixs: int | Sequence[int] | None) -> tuple[NDArr params["stat"] = "C" params["expected"] = 1.0 params["ascending"] = True + elif mode.s == "spagft": + params["func"] = _spagft + params["stat"] = "GFT" + params["expected"] = 0.0 + params["ascending"] = False else: raise NotImplementedError(f"Mode `{mode}` is not yet implemented.") diff --git a/src/squidpy/gr/_spagft.py b/src/squidpy/gr/_spagft.py new file mode 100644 index 000000000..b49d100d4 --- /dev/null +++ b/src/squidpy/gr/_spagft.py @@ -0,0 +1,43 @@ +from __future__ import annotations + +import numpy as np +from scipy.sparse import spmatrix + +from squidpy._utils import NDArrayA + + +def _spagft(g: spmatrix, vals: NDArrayA) -> NDArrayA: + """ + SpaGFT: Identify spatially variable genes using graph Fourier transform. + Returns a score per gene indicating spatial variability. + """ + from scipy.sparse import csgraph + from scipy.sparse.linalg import eigsh + + # g: adjacency matrix (n_cells x n_cells) + # vals: (n_cells x n_genes) + if vals.shape[0] != g.shape[0]: + if vals.shape[1] == g.shape[0]: + vals = vals.T + else: + raise ValueError("vals must have shape (n_cells, n_genes), where n_cells == g.shape[0].") + vals_proc = vals + + # Compute normalized Laplacian + lap = csgraph.laplacian(g, normed=True) + # Compute eigenvectors (graph Fourier basis) + n_eig = min(20, lap.shape[0] - 2) + if n_eig <= 0: + from scipy.sparse.linalg import ArpackError + + raise ArpackError("Number of eigenvectors requested must be positive.") + eigvals, eigvecs = eigsh(lap, k=n_eig, which="SM") + + # Project each gene onto Fourier basis, score by energy in low-frequency components + scores = [] + for gene in vals_proc.T: + coeffs = eigvecs.T @ gene + # SVG score: sum squared coeffs for lowest frequencies (spatially smooth signal) + lf_energy = np.sum(coeffs[: n_eig // 2] ** 2) + scores.append(lf_energy) + return np.array(scores) diff --git a/tests/graph/test_spagft.py b/tests/graph/test_spagft.py new file mode 100644 index 000000000..283d73b6e --- /dev/null +++ b/tests/graph/test_spagft.py @@ -0,0 +1,70 @@ +from __future__ import annotations + +import pytest + +from squidpy._constants._constants import SpatialAutocorr +from squidpy.gr import spatial_autocorr + + +def test_spagft_incompatible_shapes(): + import numpy as np + from scipy.sparse import lil_matrix + + from squidpy.gr._spagft import _spagft + + n = 10 + g = lil_matrix((n, n)) + for i in range(n): + g[i, (i + 1) % n] = 1 + g[i, (i - 1) % n] = 1 + g = g.tocsr() + vals = np.random.rand(5, 7) + with pytest.raises(ValueError): + _spagft(g, vals) + + +def test_spagft_svg_identification(): + import numpy as np + from anndata import AnnData + + from squidpy.gr import spatial_autocorr + + n_cells = 50 + np.random.seed(42) + spatial_pattern = np.sin(np.linspace(0, 2 * np.pi, n_cells)) + random_gene = np.random.normal(size=n_cells) + X = np.vstack([spatial_pattern, random_gene]) + adata = AnnData(X=X.T) + from scipy.sparse import lil_matrix + + g = lil_matrix((n_cells, n_cells)) + for i in range(n_cells): + g[i, (i + 1) % n_cells] = 1 + g[i, (i - 1) % n_cells] = 1 + adata.obsp["spatial_connectivities"] = g.tocsr() + df = spatial_autocorr(adata, mode="spagft", copy=True) + assert "GFT" in df.columns + assert df["GFT"].iloc[0] > df["GFT"].iloc[1] + + +def test_spagft_enum_recognition(): + # Check that the enum contains "spagft" + assert hasattr(SpatialAutocorr, "SPAGFT") + # Check that spatial_autocorr accepts the enum member + import numpy as np + from anndata import AnnData + + n_cells = 10 + np.random.seed(0) + X = np.random.normal(size=(n_cells, 2)) + adata = AnnData(X=X) + from scipy.sparse import lil_matrix + + g = lil_matrix((n_cells, n_cells)) + for i in range(n_cells): + g[i, (i + 1) % n_cells] = 1 + g[i, (i - 1) % n_cells] = 1 + adata.obsp["spatial_connectivities"] = g.tocsr() + # Should not raise + df = spatial_autocorr(adata, mode=SpatialAutocorr.SPAGFT, copy=True) + assert "GFT" in df.columns