theislab / cellrank

CellRank: dynamics from multi-view single-cell data
https://cellrank.org
BSD 3-Clause "New" or "Revised" License
350 stars 46 forks source link

A question about ‘g_fwd.compute_fate_probabilities(preconditioner='ilu', tol=1e-7, solver='direct')’ #1194

Open liuyang218622 opened 7 months ago

liuyang218622 commented 7 months ago

Hi, when I run this code below:

import sys import cellrank as cr import scanpy as sc import pandas as pd sc.settings.set_figure_params(frameon=False, dpi=100) import numpy as np cr.settings.verbosity = 2 import warnings warnings.simplefilter("ignore", category=UserWarning)

Load data

adata=sc.read_h5ad('data.h5ad') sc.pl.embedding(adata, basis='velocity_umap', color=['celltype','leiden_velocity'], legend_loc='on data')

Prepare Cellrank

sc.tl.pca(adata, n_comps=50, svd_solver = 'arpack') sc.pp.neighbors(adata, n_neighbors=15, n_pcs=40) sc.tl.diffmap(adata) root = adata.obs.index.get_loc(adata[adata.obs['celltype']=='C13'].obs.index[np.argmin(adata[adata.obs['celltype']=='C13'].obsm['X_umap'][:, 1])]) adata.uns['iroot'] = root sc.tl.dpt(adata) adata.obs['root'] = np.nan adata.obs['root'][root] = 1 sc.pl.embedding(adata, basis='velocity_umap', color=['dpt_pseudotime','celltype'],vmax=0.5)

Computing fate probability

from cellrank.kernels import PseudotimeKernel pk = PseudotimeKernel(adata,time_key="dpt_pseudotime") pk.compute_transition_matrix(threshold_scheme='soft',nu=0.01, b=20.0) ck = cr.kernels.ConnectivityKernel(adata) ck.compute_transition_matrix() kernel = 0.3ck + 0.7pk from cellrank.estimators import GPCCA g_fwd = GPCCA(kernel) print(g_fwd) g_fwd.fit(n_states=3, cluster_key="celltype")

predict_terminal_status

g_fwd.plot_macrostates(which="all") g_fwd.predict_terminal_states(method="top_n", n_states=5) g_fwd.plot_macrostates(which="terminal") CD4_T_cells = (adata[adata.obs['celltype'] == 'CD4 T cells']).obs_names[:30] CD8_T_cells = (adata[adata.obs['celltype'] == 'CD8 T cells']).obs_names[:30] Tgd_cells = (adata[adata.obs['celltype'] == 'Tgd cells']).obs_names[:30] tLable = pd.DataFrame({ 'CD4 T cells': CD4_T_cells, 'CD8 T cells': CD8_T_cells, 'Tgd cells': Tgd_cells }) tLable = tLable.to_dict(orient='list') g_fwd.set_terminal_states(states=tLable, cluster_key='celltype') g_fwd.plot_macrostates(which="terminal")

Compute fate probabilities and driver genes

g_fwd.compute_fate_probabilities(preconditioner='ilu', tol=1e-7, solver='direct')

I got this error

Computing fate probabilities 100% 3/3 [19:59<00:00, 398.50s/] WARNING: 3 solution(s) did not converge

ValueError Traceback (most recent call last) Cell In[69], line 1 ----> 1 g_fwd.compute_fate_probabilities(preconditioner='ilu', tol=1e-7, solver='direct')

File ~/anaconda3/envs/cellrank/lib/python3.10/site-packages/cellrank/estimators/mixins/_fate_probabilities.py:220, in FateProbsMixin.compute_fate_probabilities(self, keys, solver, use_petsc, n_jobs, backend, show_progress_bar, tol, preconditioner) 218 start = logg.info("Computing fate probabilities") 219 data = self._rec_trans_states(keys, ctx="fate_probs") --> 220 fate_probs = self._compute_fate_probabilities( 221 data.q, 222 data.s, 223 trans_indices=data.trans_indices, 224 term_states=data.term_states, 225 solver=solver, 226 use_petsc=use_petsc, 227 n_jobs=n_jobs, 228 backend=backend, 229 tol=tol, 230 show_progress_bar=show_progress_bar, 231 preconditioner=preconditioner, 232 ) 233 fate_probs = Lineage( 234 fate_probs, 235 names=list(data.term_states.cat.categories), 236 colors=data.term_states_colors, 237 ) 239 params = self._create_params(remove=["use_petsc", "n_jobs", "backend", "show_progress_bar"])

File ~/anaconda3/envs/cellrank/lib/python3.10/site-packages/cellrank/estimators/mixins/_fate_probabilities.py:494, in FateProbsMixin._compute_fate_probabilities(self, q, s, trans_indices, term_states, solver, use_petsc, n_jobs, backend, tol, show_progress_bar, preconditioner) 492 mask = np.isclose(abs_classes.sum(1), 1.0, rtol=1e-3) 493 if not np.all(mask): --> 494 raise ValueError( 495 f"{np.sum(~mask)} value(s) do not sum to 1 (rtol=1e-3). Try decreasing the tolerance as tol=..., " 496 f"specifying a preconditioner as preconditioner=... or " 497 f"use a direct solver as solver='direct' if the matrix is small." 498 ) 500 return abs_classes

ValueError: 4791 value(s) do not sum to 1 (rtol=1e-3). Try decreasing the tolerance as tol=..., specifying a preconditioner as preconditioner=... or use a direct solver as solver='direct' if the matrix is small.

I try to set the tol=1e-7, tol=1e-20,tol=1e-25. I also used g_fwd = GPCCA(kernel) instead of g_fwd = GPCCA(pk). But I still get the error. Can you give me some suggestions?

WeilerP commented 7 months ago

@liuyang218622, can you check values larger than 1e-7? There was a counterintuitive way to decrease the threshold; I'm just not sure if it was here or if we updated it. Please also make sure you are running the latest version of CellRank.

liuyang218622 commented 7 months ago

The cellrank version I used is v2.0.4

image

I used the tol=1e-2,tol=1e-3,tol=1e-5 but still get the same error

Computing fate probabilities WARNING: Unable to import petsc4py. For installation, please refer to: https://petsc4py.readthedocs.io/en/stable/install.html. Defaulting to 'gmres' solver. 100%  3/3 [00:00<00:00, 15.53/s]

ValueError Traceback (most recent call last) Cell In[29], line 1 ----> 1 g_fwd.compute_fate_probabilities(preconditioner='ilu', tol=1e-2, solver='direct')

File ~/anaconda3/envs/Cellrank2/lib/python3.10/site-packages/cellrank/estimators/mixins/_fate_probabilities.py:220, in FateProbsMixin.compute_fate_probabilities(self, keys, solver, use_petsc, n_jobs, backend, show_progress_bar, tol, preconditioner) 218 start = logg.info("Computing fate probabilities") 219 data = self._rec_trans_states(keys, ctx="fate_probs") --> 220 fate_probs = self._compute_fate_probabilities( 221 data.q, 222 data.s, 223 trans_indices=data.trans_indices, 224 term_states=data.term_states, 225 solver=solver, 226 use_petsc=use_petsc, 227 n_jobs=n_jobs, 228 backend=backend, 229 tol=tol, 230 show_progress_bar=show_progress_bar, 231 preconditioner=preconditioner, 232 ) 233 fate_probs = Lineage( 234 fate_probs, 235 names=list(data.term_states.cat.categories), 236 colors=data.term_states_colors, 237 ) 239 params = self._create_params(remove=["use_petsc", "n_jobs", "backend", "show_progress_bar"])

File ~/anaconda3/envs/Cellrank2/lib/python3.10/site-packages/cellrank/estimators/mixins/_fate_probabilities.py:494, in FateProbsMixin._compute_fate_probabilities(self, q, s, trans_indices, term_states, solver, use_petsc, n_jobs, backend, tol, show_progress_bar, preconditioner) 492 mask = np.isclose(abs_classes.sum(1), 1.0, rtol=1e-3) 493 if not np.all(mask): --> 494 raise ValueError( 495 f"{np.sum(~mask)} value(s) do not sum to 1 (rtol=1e-3). Try decreasing the tolerance as tol=..., " 496 f"specifying a preconditioner as preconditioner=... or " 497 f"use a direct solver as solver='direct' if the matrix is small." 498 ) 500 return abs_classes

ValueError: 6411 value(s) do not sum to 1 (rtol=1e-3). Try decreasing the tolerance as tol=..., specifying a preconditioner as preconditioner=... or use a direct solver as solver='direct' if the matrix is small.

WeilerP commented 7 months ago

@liuyang218622, could you please use CellRank with petsc and slepsc installed, and check again?

liuyang218622 commented 7 months ago

I am installing them now, but I have another question. Is this condition associated with the fact that I set the terminal status manually?

Predict_terminal_status

g_fwd.plot_macrostates(which="all") g_fwd.predict_terminal_states(method="top_n", n_states=5) g_fwd.plot_macrostates(which="terminal") CD4_T_cells = (adata[adata.obs['celltype'] == 'CD4 T cells']).obs_names[:30] CD8_T_cells = (adata[adata.obs['celltype'] == 'CD8 T cells']).obs_names[:30] Tgd_cells = (adata[adata.obs['celltype'] == 'Tgd cells']).obs_names[:30] tLable = pd.DataFrame({ 'CD4 T cells': CD4_T_cells, 'CD8 T cells': CD8_T_cells, 'Tgd cells': Tgd_cells }) tLable = tLable.to_dict(orient='list')

g_fwd.set_terminal_states(states=tLable, cluster_key='celltype')

g_fwd.plot_macrostates(which="terminal")

liuyang218622 commented 7 months ago

And I found If I don't set terminal_status manually. I can run g_fwd.compute_fate_probabilities(preconditioner='ilu', tol=1e-7, solver='direct') successfully