Gene trend plotting with pygam with scipy>1.14.0

[WeilerP]

Plotting gene trends fails with pyGAM, ATM.

import cellrank as cr
import scvelo as scv

adata = cr.datasets.bone_marrow()

scv.pp.filter_and_normalize(
    adata, min_shared_counts=20, n_top_genes=2000, subset_highly_variable=False, log=False
)
sc.pp.log1p(adata)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=30, random_state=0)

ptk = cr.kernels.PseudotimeKernel(adata, time_key="palantir_pseudotime")
ptk.compute_transition_matrix()

estimator = cr.estimators.GPCCA(ptk)
estimator.compute_schur()

estimator.compute_macrostates(5, cluster_key="clusters")
estimator.predict_terminal_states()
estimator.compute_fate_probabilities()

model = cr.models.GAM(adata)
cr.pl.gene_trends(
    adata,
    model=model,
    genes="ITGA2B",
    time_key="palantir_pseudotime",
    hide_cells=True,
    same_plot=True,
)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/models/_pygam_model.py:203, in GAM.fit(self, x, y, w, **kwargs)
    202 try:
--> 203     self.model.fit(self.x, self.y, weights=self.w, **kwargs)
    204     return self

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/pygam/pygam.py:913, in GAM.fit(self, X, y, weights)
    912 # optimize
--> 913 self._pirls(X, y, weights)
    914 # if self._opt == 0:
    915 #     self._pirls(X, y, weights)
    916 # if self._opt == 1:
    917 #     self._pirls_naive(X, y)

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/pygam/pygam.py:739, in GAM._pirls(self, X, Y, weights)
    733 if (
    734     not self._is_fitted
    735     or len(self.coef_) != self.terms.n_coefs
    736     or not np.isfinite(self.coef_).all()
    737 ):
    738     # initialize the model
--> 739     self.coef_ = self._initial_estimate(Y, modelmat)
    741 assert np.isfinite(
    742     self.coef_
    743 ).all(), "coefficients should be well-behaved, but found: {}".format(self.coef_)

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/pygam/pygam.py:706, in GAM._initial_estimate(self, y, modelmat)
    704 # solve the linear problem
    705 return np.linalg.solve(
--> 706     load_diagonal(modelmat.T.dot(modelmat).A), modelmat.T.dot(y_)
    707 )

AttributeError: 'csr_matrix' object has no attribute 'A'

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/models/_base_model.py:65, in _handle_exception.<locals>.handle.<locals>.wrapper(wrapped, instance, args, kwargs)
     64         instance.reraise()
---> 65     return wrapped(*args, **kwargs)
     66 except Exception as e:  # noqa: BLE001

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/models/_pygam_model.py:206, in GAM.fit(self, x, y, w, **kwargs)
    205 except Exception as e:  # noqa: BLE001
--> 206     raise RuntimeError(
    207         f"Unable to fit `{type(self).__name__}` for gene "
    208         f"`{self._gene!r}` in lineage `{self._lineage!r}`. Reason: `{e}`"
    209     ) from e

RuntimeError: Unable to fit `GAM` for gene `'ITGA2B'` in lineage `'Mega'`. Reason: `'csr_matrix' object has no attribute 'A'`

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
Cell In[17], line 4
      1 model = cr.models.GAM(adata)
      3 with mplscience.style_context():
----> 4     cr.pl.gene_trends(
      5         adata,
      6         model=model,
      7         genes="ITGA2B",
      8         time_key="palantir_pseudotime",
      9         hide_cells=True,
     10         same_plot=True,
     11     )
     12     plt.show()

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/_utils/_utils.py:1574, in _genesymbols(wrapped, instance, args, kwargs)
   1566 @wrapt.decorator
   1567 def _genesymbols(wrapped: Callable[..., Any], instance: Any, args: Any, kwargs: Any) -> Any:
   1568     with _gene_symbols_ctx(
   1569         args[0],
   1570         key=kwargs.pop("gene_symbols", None),
   1571         make_unique=True,
   1572         use_raw=kwargs.get("use_raw", False),
   1573     ):
-> 1574         return wrapped(*args, **kwargs)

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/pl/_gene_trend.py:207, in gene_trends(adata, model, genes, time_key, lineages, backward, data_key, time_range, transpose, callback, conf_int, same_plot, hide_cells, perc, lineage_cmap, fate_prob_cmap, cell_color, cell_alpha, lineage_alpha, size, lw, cbar, margins, sharex, sharey, gene_as_title, legend_loc, obs_legend_loc, ncols, suptitle, return_models, n_jobs, backend, show_progress_bar, figsize, dpi, save, return_figure, plot_kwargs, **kwargs)
    204 kwargs["conf_int"] = conf_int  # prepare doesnt take or need this
    205 models = _create_models(model, genes, lineages)
--> 207 all_models, models, genes, lineages = _fit_bulk(
    208     models,
    209     _create_callbacks(adata, callback, genes, lineages, **kwargs),
    210     genes,
    211     lineages,
    212     time_range,
    213     return_models=True,
    214     filter_all_failed=False,
    215     parallel_kwargs={
    216         "show_progress_bar": show_progress_bar,
    217         "n_jobs": _get_n_cores(n_jobs, len(genes)),
    218         "backend": _get_backend(models, backend),
    219     },
    220     **kwargs,
    221 )
    223 lineages = sorted(lineages)
    224 probs = probs[lineages]

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/pl/_utils.py:320, in _fit_bulk(models, callbacks, genes, lineages, time_range, parallel_kwargs, return_models, filter_all_failed, **kwargs)
    304 models = parallelize(
    305     _fit_bulk_helper,
    306     genes,
   (...)
    316     **kwargs,
    317 )
    318 logg.info("    Finish", time=start)
--> 320 return _filter_models(models, return_models=return_models, filter_all_failed=filter_all_failed)

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/pl/_utils.py:356, in _filter_models(models, return_models, filter_all_failed)
    354         for model in ms.values():
    355             assert isinstance(model, FailedModel), f"Expected `FailedModel`, found `{type(model).__name__!r}`."
--> 356             model.reraise()
    358 if not np.all(modelmask.values):
    359     failed_models = modelmat.values[~modelmask.values]

File ~/miniforge3/envs/crp-py310/lib/python3.10/site-packages/cellrank/models/_base_model.py:1234, in FailedModel.reraise(self)
   1232 """Raise the original exception with additional model information."""
   1233 # retain the exception type and also the original exception
-> 1234 raise type(self._exc)(f"Fatal model failure `{self}`.") from self._exc

RuntimeError: Fatal model failure `<FailedModel[origin=GAM[gene='ITGA2B', lineage='Mega', model=GammaGAM(callbacks=[Deviance(), Diffs()], fit_intercept=True, max_iter=2000, scale=None, terms=s(0) + intercept, tol=0.0001, verbose=False)]]>`.

Versions:

cellrank==2.0.6 scanpy==1.10.4 anndata==0.11.1 numpy==1.26.4 numba==0.60.0 scipy==1.14.1 pandas==2.2.3 pygpcca==1.0.4 scikit-learn==1.6.0 statsmodels==0.14.4 python-igraph==0.11.8 scvelo==0.3.3 pygam==0.9.1 matplotlib==3.10.0 seaborn==0.13.2

[WeilerP]

@michalk8, can we offer an easy solution here without pinning scipy and waiting for pyGAM release?