Fix adata.raw handling after slicing (Python ≥3.11 / 3.13)

.github/workflows/test_and_deploy.yml

@@ -21,7 +21,7 @@ jobs:
     strategy:
       matrix:
         platform: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: [3.8, 3.9, '3.10']
+        python-version: ['3.11', '3.12', '3.13']
 
     steps:
       - uses: actions/checkout@v2

pyproject.toml

@@ -1,11 +1,59 @@
 [build-system]
-requires = ["setuptools", "wheel"]
+requires = ["setuptools>=42.0.0", "wheel"]
 build-backend = "setuptools.build_meta"
 
-[tool.briefcase]
-project_name = "sc-3D"
-author = "Leo Guignard"
+[tool.setuptools.packages.find]
+where = ["src"]
+
+[project]
+authors = [
+    { name = "Léo Guignard", email = "leo.guignard@univ-amu.fr"},
+    { name = "Miquel Sendra"},
+]
+maintainers = [
+    {name = "Léo Guignard", email = "leo.guignard@univ-amu.fr"}
+]
+name = "sc-3D"
+description = "Array alignment and 3D differential expression for 3D spatial omics"
+version = "2.0.0"
 license = "MIT"
+license-files = [ "LICENSE" ]
+readme = {file = "README.md", content-type = "text/markdown"}
+requires-python = ">= 3.11"
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Intended Audience :: Developers",
+    "Operating System :: OS Independent",
+    "Programming Language :: Python",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3 :: Only",
+    "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+]
+
+dependencies = [
+    "scipy",
+    "numpy",
+    "matplotlib",
+    "pandas",
+    "seaborn",
+    "scikit-learn",
+    "anndata",
+]
+
+[project.optional-dependencies]
+testing = [
+    "tox",
+    "pytest",
+    "pytest-cov",
+]
+
+[project.urls]
+"Bug Tracker" = "https://github.com/GuignardLab/sc3D/issues"
+"Documentation" = "https://github.com/GuignardLab/sc3D#README.md"
+"Source Code" = "https://github.com/GuignardLab/sc3D"
+"User Support" = "https://github.com/GuignardLab/sc3D/issues"
 
 [tool.black]
 line-length = 79
@@ -25,13 +73,11 @@ push = false
 [tool.bumpver.file_patterns]
 "pyproject.toml" = [
     'current_version = "{version}"',
+    'version = "{version}"',
 ]
 "src/sc3D/__init__.py" = [
     '__version__ = "{version}"',
 ]
-"setup.cfg" = [
-    'version = {version}',
-]
 "CITATION.cff" = [
     "version: {version}",
 ]

src/sc3D/_tests/test_sc3D.py

@@ -1,33 +1,67 @@
 from sc3D import SpatialOmicArray
 import numpy as np
 
+em = SpatialOmicArray("data/data_test.h5ad", store_anndata=True)
+
+em2 = SpatialOmicArray(
+    "data/DLPFC.h5ad",
+    tissue_id="layer_guess",
+    pos_id="spatial",
+    array_id="z",
+    z_space=30,
+    store_anndata=True,
+)
+
 
 def test_sc3D():
-    em = SpatialOmicArray("data/data_test.h5ad", store_anndata=True)
     assert len(em.all_cells) == 120
+
+
+def test_smooth():
     em.smooth_data()
+
+
+def test_plot_coverslip():
     em.plot_coverslip(7)
+
+
+def test_3D_diff():
     em.get_3D_differential_expression([21])
+
+
+def test_plot_diff():
     em.plot_top_3D_diff_expr_genes([21, 23], repetition_allowed=True)
     em.plot_top_3D_diff_expr_genes([21, 23], repetition_allowed=False)
+
+
+def test_vol_neighbs():
     em.plot_volume_vs_neighbs(21)
+
+
+def test_spatial_outliers():
     em.removing_spatial_outliers()
+
+
+def test_volumes():
     em.compute_volumes()
+
+
+def test_z_pos():
     em.set_zpos()
 
-    em = SpatialOmicArray(
-        "data/DLPFC.h5ad",
-        tissue_id="layer_guess",
-        pos_id="spatial",
-        array_id="z",
-        z_space=30,
-        store_anndata=True,
-    )
-    em.produce_em()
-    em.registration_3d()
-    origin = np.mean([em.final[c] for c in em.all_cells], axis=0)
+
+def test_produce():
+    em2.produce_em()
+
+
+def test_registration():
+    em2.registration_3d()
+
+
+def test_plot_slices():
+    origin = np.mean([em2.final[c] for c in em2.all_cells], axis=0)
     origin = np.hstack([origin, 80])
     angles = np.array([-5.0, 5.0, 0.0])
-    _ = em.plot_slice(
+    _ = em2.plot_slice(
         angles, color_map="viridis", origin=origin, thickness=30, nb_interp=5
     )

src/sc3D/sc3D.py

@@ -17,7 +17,9 @@
 from scipy.spatial.distance import cdist
 from scipy.optimize import linear_sum_assignment
 from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
-from scipy.stats import zscore, linregress
+from scipy.stats import zscore
+from scipy.stats import linregress
+from scipy.sparse import issparse
 from seaborn import scatterplot
 import json
 from pathlib import Path
@@ -148,10 +150,10 @@ def read_anndata(
                 + orig[-self.nb_CS_end_ignore :]
             )
             data = data[~(data.obs[array_id].isin(cs_to_remove))]
-        if data.raw is not None:
-            data.raw = data.raw.to_adata()
-        else:
-            data.raw = data.copy()
+        # if data.raw is not None:
+        #     data.raw = data.raw.to_adata()
+        # else:
+        #     data.raw = data.copy()
         ids = range(len(data))
         self.all_cells = list(ids)
         self.cell_names = dict(
@@ -196,7 +198,7 @@ def read_anndata(
             self.gene_expression = {id_: [] for id_ in ids}
 
         self.array_id_num_pos = array_id_num_pos
-        if array_id in data.obs_keys():
+        if array_id in data.obs:
             if data.obs[array_id].dtype != int:
                 exp = re.compile("[0-9]+")
                 cs = list(
@@ -724,12 +726,10 @@ def smooth_data(self, inplace=True):
         dist_sum = GG.sum(axis=1)
         product_n = product / dist_sum.reshape(-1, 1)
         product_sparse = sp.sparse.csr_array(product_n)
-        tmp_raw = self.anndata.raw.to_adata()
-        tmp_raw.X = product_sparse.toarray()
-        if inplace:
-            self.anndata.raw = tmp_raw
-        else:
-            return tmp_raw
+        # tmp_raw = self.anndata.raw.to_adata()
+        print(f"{self.anndata.X.shape=}\n{product_sparse.toarray().shape=}")
+        self.anndata.raw._X = product_sparse.toarray()
+        return self.anndata
 
     def downsample(self, spacing=10, pos_id="pos_3D"):
         """
@@ -991,7 +991,6 @@ def removing_spatial_outliers(self, th=0.2, n_components=3):
         l_all = list(self.all_cells)
         if hasattr(self, "anndata"):
             self.anndata = self.anndata[l_all]
-            self.anndata.raw = self.anndata.raw.to_adata()
         for t, c in self.cells_from_cover_slip.items():
             c.intersection_update(self.filtered_cells)
         for t, c in self.cells_from_tissue.items():
@@ -1096,9 +1095,9 @@ def registration_3d(
                     "no filtering will be applied"
                 )
             if work_with_raw:
-                raw_data = self.anndata.raw.to_adata()
+                raw_data = self.anndata.raw
             else:
-                raw_data = self.anndata.copy()
+                raw_data = self.anndata
             if sc_imp:
                 if min_counts_genes is not None:
                     filter_1 = sc.pp.filter_genes(
@@ -1810,7 +1809,56 @@ def save_anndata(self, output_path):
         Args:
             output_path (str): path to the output anndata file ('.h5ad' file)
         """
-        data_tmp = self.anndata.copy()
+        a = self.anndata
+
+        if a.is_view:
+            ref = a._adata_ref
+
+            rows = ref.obs_names.get_indexer(a.obs_names)
+            cols = ref.var_names.get_indexer(a.var_names)
+            if (rows < 0).any() or (cols < 0).any():
+                raise ValueError(
+                    "Could not map view obs/var names back to parent AnnData."
+                )
+
+            X = ref.X[rows, :]
+            X = X[:, cols]
+            if issparse(X):
+                X = X.tocsr()
+
+            data_tmp = anndata.AnnData(X=X, obs=a.obs.copy(), var=a.var.copy())
+            data_tmp.uns = dict(a.uns)
+
+            # --- IMPORTANT: keep all genes accessible for the napari viewer ---
+            # The viewer expects embryo.anndata.raw to contain the full gene set.
+            if ref.raw is not None:
+                Xraw = ref.raw.X[rows, :]  # keep ALL genes (no col subset)
+                if issparse(Xraw):
+                    Xraw = Xraw.tocsr()
+                raw_tmp = anndata.AnnData(
+                    X=Xraw, obs=a.obs.copy(), var=ref.raw.var.copy()
+                )
+                data_tmp.raw = raw_tmp
+
+        else:
+            data_tmp = a  # already materialized
+            # Ensure raw exists for viewer if possible
+            if (
+                data_tmp.raw is None
+                and getattr(a, "_adata_ref", None) is not None
+                and a._adata_ref.raw is not None
+            ):
+                ref = a._adata_ref
+                rows = ref.obs_names.get_indexer(a.obs_names)
+                Xraw = ref.raw.X[rows, :]
+                if issparse(Xraw):
+                    Xraw = Xraw.tocsr()
+                raw_tmp = anndata.AnnData(
+                    X=Xraw, obs=a.obs.copy(), var=ref.raw.var.copy()
+                )
+                data_tmp.raw = raw_tmp
+
+        # add registered coords
         all_c_sorted = sorted(self.all_cells)
         pos_final = np.array([self.pos_3D[c] for c in all_c_sorted])
         data_tmp.obsm["X_spatial_registered"] = pos_final

tox.ini

@@ -1,13 +1,13 @@
 # For more information about tox, see https://tox.readthedocs.io/en/latest/
 [tox]
-envlist = py{38,39,310}-{linux,macos,windows}
+envlist = py{311,312,313}-{linux,macos,windows}
 isolated_build=true
 
 [gh-actions]
 python =
-    3.8: py38
-    3.9: py39
-    3.10: py310
+    3.11: py311
+    3.12: py312
+    3.13: py313
 
 [gh-actions:env]
 PLATFORM =