Closed MoLuLuMo closed 3 years ago
@Lobbii, seems like 'ILIRN'
is either not in adata.var_names
or no model has been fitted.
@Lobbii, seems like
'ILIRN'
is either not inadata.var_names
or no model has been fitted.
I have run the dynamic model, how to perform the fitting ?
I have run the script shown below:
scv.pp.filter_and_normalize(WT_WT, min_shared_counts=20, n_top_genes=2000) scv.pp.moments(WT_WT, n_pcs=30, n_neighbors=30) scv.tl.recover_dynamics(WT_WT) scv.tl.velocity(WT_WT, mode='dynamical') scv.tl.velocity_graph(WT_WT) scv.tl.latent_time(WT_WT) scv.pl.scatter(WT_WT, color='latent_time', color_map='gnuplot', size=80) top_genes = WT_WT.var['fit_likelihood'].sort_values(ascending=False).index[:300] scv.pl.heatmap(WT_WT, var_names=top_genes, sortby='latent_time', col_color='SCT_snn_res.0.2', n_convolve=100)
Sorry, not sure I understand your question correctly. How is this related to the issue? scv.tl.recover_dynamics
fits the models. Please have a look at our docs and the tutorials shown there.
Sorry, not sure I understand your question correctly. How is this related to the issue?
scv.tl.recover_dynamics
fits the models. Please have a look at our docs and the tutorials shown there.
I followed the tutorial and run the dynamical model and calculate latent time. If I need to perform the fitting, how to do it? Which function should I use?
Thanks
I followed the tutorial and run the dynamical model and calculate latent time. If I need to perform the fitting, how to do it? Which function should I use?
As mentioned here: scv.tl.recover_dynamics
is used to fit the model.
I followed the tutorial and run the dynamical model and calculate latent time. If I need to perform the fitting, how to do it? Which function should I use?
As mentioned here:
scv.tl.recover_dynamics
is used to fit the model.
I have performed scv.tl.recover_dynamics, and I found the fitting information in anndata.var. In this case, how to generate the scatter plot?
AnnData object with n_obs × n_vars = 3129 × 2000 obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'Cells', 'Update_Cells', 'nCount_spliced', 'nFeature_spliced', 'nCount_unspliced', 'nFeature_unspliced', 'nCount_ambiguous', 'nFeature_ambiguous', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size', 'n_counts', 'velocity_self_transition', 'root_cells', 'end_points', 'velocity_pseudotime', 'latent_time' var: 'features', 'SCT_features', 'ambiguous_features', 'spliced_features', 'unspliced_features', 'means', 'dispersions', 'dispersions_norm', 'highly_variable', 'fit_r2', 'fit_alpha', 'fit_beta', 'fit_gamma', 'fitt', 'fit_scaling', 'fit_std_u', 'fit_std_s', 'fit_likelihood', 'fit_u0', 'fit_s0', 'fit_pval_steady', 'fit_steady_u', 'fit_steady_s', 'fit_variance', 'fit_alignment_scaling', 'velocity_genes' uns: 'pca', 'neighbors', 'recover_dynamics', 'velocity_params', 'velocity_graph', 'velocity_graph_neg' obsm: 'X_umap', 'X_pca', 'velocity_umap' varm: 'PCs', 'loss' layers: 'SCT', 'ambiguous', 'spliced', 'unspliced', 'Ms', 'Mu', 'fit_t', 'fit_tau', 'fittau', 'velocity', 'velocity_u' obsp: 'distances', 'connectivities'
Please have a look at this tutorial.
I followed the tutorial, and I have run the dynamical model to get top-likelihood genes based on fit_likelihood. Is there any thing should I check? Thanks
Please have a look at this tutorial.
Besides, I can generate scatter for dynamics genes to show the learned dynamics (spliced vs unspliced). I have checked the tutorial multiple times, and followed the colab jupyter notebook. But I don't know why I cannot generate the scatter plot to show transcriptional switches.
Please have a look at this tutorial.
Can I know which adata.var (for var_names ) is used for scv.pl.scatter(WT_WT, x='latent_time', y= var_names , frameon=False)?
@WeilerP
Package version:
scv.logging.print_versions()
scvelo==0.2.3 scanpy==1.7.2 anndata==0.7.5 loompy==3.0.6 numpy==1.20.3 scipy==1.6.3 matplotlib==3.4.2 sklearn==0.24.2 pandas==1.2.4
Here is my script
import scvelo as scv
import scanpy as sc
import numpy as np
scv.logging.print_version()
scv.settings.verbosity = 3
scv.settings.presenter_view = True
scv.settings.set_figure_params('scvelo')
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
Filtered out 6699 genes that are detected 20 counts (shared).
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set enforce=True
.
Normalized count data: spliced, unspliced.
Extracted 2000 highly variable genes.
WARNING: Did not modify X as it looks preprocessed already.
computing neighbors
finished (0:00:23) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:01) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
scv.tl.recover_dynamics(adata)
recovering dynamics (using 1/32 cores)
finished (0:04:26) --> added
'fit_pars', fitted parameters for splicing dynamics (adata.var)
scv.tl.latent_time(adata)
computing velocities
finished (0:00:01) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph
finished (0:00:02) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing terminal states
identified 0 region of root cells and 1 region of end points .
finished (0:00:00) --> added
'root_cells', root cells of Markov diffusion process (adata.obs)
'end_points', end points of Markov diffusion process (adata.obs)
WARNING: No root cells detected. Consider specifying root cells to improve latent time prediction.
computing latent time using root_cells as prior
finished (0:00:01) --> added
'latent_time', shared time (adata.obs)
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[1], frameon=False)
scv.pl.scatter(adata, x='latent_time', y=top_genes[1], frameon=False)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local)
200 if self.has_resolvers:
--> 201 return self.resolvers[key]
202
~/.conda/envs/py_seq/lib/python3.8/collections/__init__.py in __getitem__(self, key)
897 pass
--> 898 return self.__missing__(key) # support subclasses that define __missing__
899
~/.conda/envs/py_seq/lib/python3.8/collections/__init__.py in __missing__(self, key)
889 def __missing__(self, key):
--> 890 raise KeyError(key)
891
KeyError: 'GPR34'
During handling of the above exception, another exception occurred:
KeyError Traceback (most recent call last)
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local)
211 # e.g., df[df > 0]
--> 212 return self.temps[key]
213 except KeyError as err:
KeyError: 'GPR34'
The above exception was the direct cause of the following exception:
UndefinedVariableError Traceback (most recent call last)
<ipython-input-9-0f4129e17c43> in <module>
----> 1 scv.pl.scatter(adata, x='latent_time', y=top_genes[1], frameon=False)
~/.local/lib/python3.8/site-packages/scvelo/plotting/scatter.py in scatter(adata, basis, x, y, vkey, color, use_raw, layer, color_map, colorbar, palette, size, alpha, linewidth, linecolor, perc, groups, sort_order, components, projection, legend_loc, legend_loc_lines, legend_fontsize, legend_fontweight, legend_fontoutline, xlabel, ylabel, title, fontsize, figsize, xlim, ylim, add_density, add_assignments, add_linfit, add_polyfit, add_rug, add_text, add_text_pos, add_outline, outline_width, outline_color, n_convolve, smooth, rescale_color, color_gradients, dpi, frameon, zorder, ncols, nrows, wspace, hspace, show, save, ax, **kwargs)
431 obs = adata.obs[obs_keys]
432 x = obs.astype(np.float32).eval(x)
--> 433 y = obs.astype(np.float32).eval(y)
434 else:
435 raise ValueError(
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/frame.py in eval(self, expr, inplace, **kwargs)
3597 kwargs["resolvers"] = kwargs.get("resolvers", ()) + tuple(resolvers)
3598
-> 3599 return _eval(expr, inplace=inplace, **kwargs)
3600
3601 def select_dtypes(self, include=None, exclude=None) -> DataFrame:
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/eval.py in eval(expr, parser, engine, truediv, local_dict, global_dict, resolvers, level, target, inplace)
340 )
341
--> 342 parsed_expr = Expr(expr, engine=engine, parser=parser, env=env)
343
344 # construct the engine and evaluate the parsed expression
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in __init__(self, expr, engine, parser, env, level)
796 self.parser = parser
797 self._visitor = PARSERS[parser](self.env, self.engine, self.parser)
--> 798 self.terms = self.parse()
799
800 @property
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in parse(self)
815 Parse an expression.
816 """
--> 817 return self._visitor.visit(self.expr)
818
819 @property
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit(self, node, **kwargs)
399 method = "visit_" + type(node).__name__
400 visitor = getattr(self, method)
--> 401 return visitor(node, **kwargs)
402
403 def visit_Module(self, node, **kwargs):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit_Module(self, node, **kwargs)
405 raise SyntaxError("only a single expression is allowed")
406 expr = node.body[0]
--> 407 return self.visit(expr, **kwargs)
408
409 def visit_Expr(self, node, **kwargs):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit(self, node, **kwargs)
399 method = "visit_" + type(node).__name__
400 visitor = getattr(self, method)
--> 401 return visitor(node, **kwargs)
402
403 def visit_Module(self, node, **kwargs):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit_Expr(self, node, **kwargs)
408
409 def visit_Expr(self, node, **kwargs):
--> 410 return self.visit(node.value, **kwargs)
411
412 def _rewrite_membership_op(self, node, left, right):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit(self, node, **kwargs)
399 method = "visit_" + type(node).__name__
400 visitor = getattr(self, method)
--> 401 return visitor(node, **kwargs)
402
403 def visit_Module(self, node, **kwargs):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/expr.py in visit_Name(self, node, **kwargs)
533
534 def visit_Name(self, node, **kwargs):
--> 535 return self.term_type(node.id, self.env, **kwargs)
536
537 def visit_NameConstant(self, node, **kwargs):
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/ops.py in __init__(self, name, env, side, encoding)
84 tname = str(name)
85 self.is_local = tname.startswith(LOCAL_TAG) or tname in DEFAULT_GLOBALS
---> 86 self._value = self._resolve_name()
87 self.encoding = encoding
88
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/ops.py in _resolve_name(self)
101
102 def _resolve_name(self):
--> 103 res = self.env.resolve(self.local_name, is_local=self.is_local)
104 self.update(res)
105
~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local)
215 from pandas.core.computation.ops import UndefinedVariableError
216
--> 217 raise UndefinedVariableError(key, is_local) from err
218
219 def swapkey(self, old_key: str, new_key: str, new_value=None):
UndefinedVariableError: name 'GPR34' is not defined
Does the issue persist when you calculate velocities and the velocity graph explicitly, i.e.
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
scv.tl.latent_time(adata)
?
Just ran your pipeline on the Pancreas dataset used in the tutorials and everything works fine. Do you encounter the same problem there?
Does the issue persist when you calculate velocities and the velocity graph explicitly, i.e.
scv.tl.recover_dynamics(adata) scv.tl.velocity(adata, mode='dynamical') scv.tl.velocity_graph(adata) scv.tl.latent_time(adata)
?
No, it works well. I can generate stream plot of velocity
Just ran your pipeline on the Pancreas dataset used in the tutorials and everything works fine. Do you encounter the same problem there?
I didn't encounter the same problem. I converted my dataset from Seurat Object to Anndata.
No, it works well. I can generate stream plot of velocity
So after running the explicit velocity calculation you can plot the gene expression vs. latent time?
I didn't encounter the same problem. I converted my dataset from Seurat Object to Anndata.
Just to confirm:
cells x genes
form, right?Would you please check immediately prior to the line causing the problem that
'latent_time' in list(adata.obs.keys()) + list(adata.layers.keys())
'GPR34' in adata.var_names
both return True
?Also,
scv.pl.scatter(adata, x='latent_time', y=top_genes[:4], frameon=False)
does not work as well, I presume?
No, it works well. I can generate stream plot of velocity
So after running the explicit velocity calculation you can plot the gene expression vs. latent time?
I can generate the plot (scatter & heatmap) the gene expression vs. latent time for scv.datasets.pancreas()
I didn't encounter the same problem. I converted my dataset from Seurat Object to Anndata.
Just to confirm:
- You ran the tutorial notebook in the same environment?
- Your data is in
cells x genes
form, right?
I ran the tutorial notebook in the same environment. My data is in 'cells * genes' form
AnnData object with n_obs × n_vars = 3129 × 2000
(3129 cells * 2000 highly_variable genes
Would you please check immediately prior to the line causing the problem that
'latent_time' in list(adata.obs.keys()) + list(adata.layers.keys())
'GPR34' in adata.var_names
both returnTrue
?
yes, Both return True
Also,
scv.pl.scatter(adata, x='latent_time', y=top_genes[:4], frameon=False)
does not work as well, I presume?
It doesn't work as well. If I added layer='spliced' argument, I can generate a scatterplot which is different from plots shown in tutorial.
Okay, if both return True
, the correct if clause here should be entered since is_timeseries
evaluates to True
and not this case. So it seems
As a sanity check: Could you please check that scv.pl.scatter(adata, x='latent_time', y='GPR34', frameon=False)
fails as well? Also that top_genes[1] in adata.var_names
also returns True
?
Okay, if both return
True
, the correct if clause here should be entered sinceis_timeseries
evaluates toTrue
and not this case. So it seemsAs a sanity check: Could you please check that
scv.pl.scatter(adata, x='latent_time', y='GPR34', frameon=False)
fails as well? Also thattop_genes[1] in adata.var_names
also returnsTrue
?
top_genes[1]
GPR34
top_genes[1] in adata.var_names
True
scv.pl.scatter(adata, x='latent_time', y='GPR34', frameon=False)
KeyError Traceback (most recent call last) ~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local) 200 if self.has_resolvers: --> 201 return self.resolvers[key] 202
~/.conda/envs/py_seq/lib/python3.8/collections/init.py in getitem(self, key) 897 pass --> 898 return self.missing(key) # support subclasses that define missing 899
~/.conda/envs/py_seq/lib/python3.8/collections/init.py in missing(self, key) 889 def missing(self, key): --> 890 raise KeyError(key) 891
KeyError: 'GPR34'
Hm, okay, sorry I have no idea why this is happening then. Would have to debug this in the code base. Unless you can share your h5ad file, you'll have to clone the scvelo repo and debug why this evaluates to False
.
@WeilerP
I used the example dataset from http://htmlpreview.github.io/?https://github.com/satijalab/seurat-wrappers/blob/master/docs/scvelo.html. I performed seurat pipeline for UMAP and clustering and converted seurat object to Anndata. Then I had the same issue for scatter plot.
@Lobbii, couldn't run the R code out of the box and do not have the time to debug it. Can you provide the h5ad file you are working with?
Before you do so: There seem to have been some problems with pandas==1.3.0
(see e.g. here). Could you check if rolling back your version of pandas fixes the issue?
@Lobbii, any update on this?
https://drive.google.com/drive/folders/10DYf7Sc9tqyfElIMXaHKxchRjGvfx-hM?usp=sharing
I followed the tutorial provided by Seurat to generate H5ad file
library(Seurat)
library(SeuratDisk)
library(SeuratWrappers)
ldat <- ReadVelocity(file = "./SCG71.loom")
bm <- as.Seurat(x = ldat)
bm[["RNA"]] <- bm[["spliced"]]
bm <- SCTransform(bm)
bm <- RunPCA(bm)
bm <- RunUMAP(bm, dims = 1:20)
bm <- FindNeighbors(bm, dims = 1:20)
bm <- FindClusters(bm)
DefaultAssay(bm) <- "RNA"
SaveH5Seurat(bm, filename = "mouseBM.h5Seurat")
Convert("mouseBM.h5Seurat", dest = "h5ad")
@Lobbii, the problem are the variable names in adata.raw
:
>>> adata.raw.var_names
Index(['0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
...
'24411', '24412', '24413', '24414', '24415', '24416', '24417', '24418',
'24419', '24420'],
dtype='object', length=24421)
Everything works with use_raw=False
, i.e.,
scv.pl.scatter(adata, x='latent_time', y=top_genes[1], frameon=False, use_raw=False)
scv.pl.scatter cannot generate plot to show transcriptional switches. I convert seurat object to anndata, and then perform rna velocity (dynamical model) through scVelo. I can generate the heatmap of top likelihood genes.
Error
```pytb KeyError Traceback (most recent call last) ~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local) 200 if self.has_resolvers: --> 201 return self.resolvers[key] 202 ~/.conda/envs/py_seq/lib/python3.8/collections/__init__.py in __getitem__(self, key) 897 pass --> 898 return self.__missing__(key) # support subclasses that define __missing__ 899 ~/.conda/envs/py_seq/lib/python3.8/collections/__init__.py in __missing__(self, key) 889 def __missing__(self, key): --> 890 raise KeyError(key) 891 KeyError: 'IL1RN' During handling of the above exception, another exception occurred: KeyError Traceback (most recent call last) ~/.conda/envs/py_seq/lib/python3.8/site-packages/pandas/core/computation/scope.py in resolve(self, key, is_local) 211 # e.g., df[df > 0] --> 212 return self.temps[key] 213 except KeyError as err: KeyError: 'IL1RN' The above exception was the direct cause of the following exception: UndefinedVariableError Traceback (most recent call last)Versions:
scvelo==0.2.3 scanpy==1.7.2 anndata==0.7.5 > ![download](https://user-images.githubusercontent.com/25945712/123481309-8b7c6a80-d5d1-11eb-9a19-7e9e3bf9ac9b.png)