Gene trend plotting with pygam with scipy>1.14.0
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
@michalk8, can we offer an easy solution here without pinning scipy and waiting for pyGAM release?
I am getting the same errors for gene trends and then for the drivers heatmap; I have tried Python 3.11 and 3.12 starting from a new environment installing cellrank with conda install -c conda-forge cellrank.
What's the solution at the moment? I would really like to get the heatmaps out in particular
@georginaalbadri, the issue is not the Python version but scipy. After installing CellRank, you can reinstall scipy<1.14.0 via pip install "scipy<1.14", for example.
Thanks for your swift response, I'd tried 1.14.0 instead of 1.14.1 but not below, it's working now 👍