statbiomed / unitvelo Goto Github PK

View Code? Open in Web Editor NEW

23.0 3.0 9.0 20.72 MB

UniTVelo, Temporally Unified RNA Velocity for single cell trajectory inference

Home Page: https://unitvelo.readthedocs.io/en/latest/

License: BSD 3-Clause "New" or "Revised" License

Jupyter Notebook 98.38% Python 1.62%

rna velocity transcriptomics trajectory inference

unitvelo's Introduction

UniTVelo for RNA Velocity Analysis

Temporally unified RNA velocity for single cell trajectory inference (UniTVelo) is implementated on Python 3 and TensorFlow 2. The model estimates velocity of each gene and updates cell time based on phase portraits concurrently.

The major features of UniTVelo are,

Using spliced RNA oriented design to model RNA velocity and transcription rates
Introducing a unified latent time (Unified-time mode) across whole transcriptome to incorporate stably and monotonically changed genes
Retaining gene-spcific time matrics (Independent mode) for complex datasets

UniTVelo has proved its robustness in 10 different datasets. Details can be found via our manuscript in bioRxiv which is currently under review (UniTVelo).

Installation

GPU Acceleration

UniTVelo is designed based on TensorFlow's automatic differentiation architecture. Please make sure TensorFlow 2 and relative CUDA dependencies are correctly installed.

Use the following scripts to confirm TensorFlow is using the GPU.

import tensorflow as tf
print ("Num GPUs Available: ", len(tf.config.list_physical_devices('GPU')))

If GPU is not available, UniTVelo will automatically switch to CPU for model fitting or it can be spcified in config.py (see Getting Started below).

Main Module

(Optional) Create a separate conda environment for version control and to avoid potential conflicts.

conda create -n unitvelo python=3.7
conda activate unitvelo

UniTVelo package can be conveniently installed via PyPI or directly from GitHub repository.

pip install unitvelo

pip install git+https://github.com/StatBiomed/UniTVelo

Getting Started

Public Datasets

Examples of UniTVelo and steps for reproducible results are provided in Jupyter Notebook under notebooks folder. Specifically, please refer to records analyzing Mouse Erythroid and Human Bone Marrow datasets.

UniTVelo has proved its performance through 10 different datasets and 4 of them have been incorporated within scVelo package, see datasets. Others can be obtained via link.

RNA Velocity on New Dataset

UniTVelo provides an integrated function for velocity analysis by default whilst specific configurations might need to be adjusted accordingly.

Import package

import unitvelo as utv

Sub-class and override base configuration file (here lists a few frequently used), please refer config.py for detailed arguments.

velo = utv.config.Configuration()
velo.R2_ADJUST = True 
velo.IROOT = None
velo.FIT_OPTION = '1'
velo.GPU = 0

Arguments:
- -- velo.R2_ADJUST (bool), linear regression R-squared on extreme quantile (default) or full data (adjusted)
- -- velo.IROOT (str), specify root cell cluster would enable diffusion map based time initialization, default None
- -- velo.FIT_OPTION (str), '1' Unified-time mode (default), '2' Independent mode
- -- velo.GPU (int), specify the GPU card used for fitting, -1 will switch to CPU mode, default 0.

Run model (label refers to column name in adata.obs specifying celltypes)

adata = utv.run_model(path_to_adata, label, config_file=velo)
scv.pl.velocity_embedding_stream(adata, color=label, dpi=100, title='')

Evaluation metrics (Optional)

# Cross Boundary Direction Correctness
# Ground truth should be given via `cluster_edges`
metrics = {}
metrics = utv.evaluate(adata, cluster_edges, label, 'velocity')

# Latent time estimation
scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=20)

# Phase portraits for individual genes (experimental)
utv.pl.plot_range(gene_name, adata, velo, show_ax=True, time_metric='latent_time')

unitvelo's People

Contributors

Stargazers

Watchers

Forkers

michaelgmz shobistassen milescsmith altairwei auniany idriskb molecules hoimanc dp-sun

unitvelo's Issues

Error in low UMI count dataset

Hi!

I tried to run UniTVelo on the MouseBoneMarrow dataset with very low UMI counts that you gave to me at some point. I found that that while the "unified" model works the "independent model" gives the below error. Do you know what to do about this?

Thanks a lot!
Best wishes,

Alexander

Preprint link not working

Hi all, this looks interesting but the link to your preprint is broken - can you please point us to it? Thanks!

Ping @michalk8

Why does UniTVelo v0.2.4.2 depend on a specified pandas version (pandas==1.4.2)

Hi,

Thank you for this fantastic tool! I found a dependency problem when attempted to upgrade UniTVelo from 0.2.3 to 0.2.4. My project depends on the scprep package which requires a certain version (<1.4,>=0.25) of pandas. However, in the commit 927bc2a , a specified pandas version (v1.4.2) was added to setup.py , which conflicts with scprep.

Is this a bug? Will I miss critical features if I stay on version 0.2.3?

Best,
Altair

peak time fit_t0 not found

I'm trying to use the var of "fit_t0" to identify repression/induction genes, but this var is not found. I found there is a var of "fit_t", does it equals to "fit_t0"?

Thank you for your awesome tool and answering my questions!

TypeError: 'NoneType' object is not iterable

Hi, I'm running UniTVelo and obtained some strange results with the default settings. Specifically, the velocity direction was inverted compared to scvelo and biologically feasible results. However, when I modify the default settings, the following error appears. Any ideas?

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/optimize_utils.py:187, in Model_Utils.gene_prior_perc(self, dis)
    184 perc, perc_idx = [], []
    185 aggregrate_weight = np.ones((dis.shape[1], ), dtype=np.float32)
--> 187 for prior in self.config.GENE_PRIOR:
    188     expr = np.array(self.adata[:, prior[0]].layers['Ms'])
    189     perc.append(np.max(expr) * 0.75) # modify 0.75 for parameter tuning

TypeError: 'NoneType' object is not iterable

Parameter choices are as follows:

velo_config = utv.config.Configuration()
velo_config.R2_ADJUST = True
velo_config.IROOT = None
velo_config.FIT_OPTION = '2' ## Default is 1
velo_config.AGENES_R2 = 1

System and version information:

unitvelo==0.2.5
DISTRIB_DESCRIPTION="Ubuntu 20.04.5 LTS"
Linux version 5.10.16.3-microsoft-standard-WSL2
Python 3.9.15

The full error message:

------> Manully Specified Parameters <------
FIT_OPTION:	2
DENSITY:	Raw
REORDER_CELL:	Hard
AGGREGATE_T:	False
------> Model Configuration Settings <------
N_TOP_GENES:	2000
LEARNING_RATE:	0.01
R2_ADJUST:	True
GENE_PRIOR:	None
VGENES:	basic
IROOT:	None
--------------------------------------------

Current working dir is /mnt/c/workdir.
Results will be stored in res folder
Extracted 2000 highly variable genes.
Computing moments for 2000 genes with n_neighbors: 30 and n_pcs: 30

# of velocity genes 561 (Criterion: positive regression coefficient between un/spliced counts)
# of velocity genes 535 (Criterion: std of un/spliced reads should be moderate, w/o extreme values)
# of velocity genes 535 (Criterion: genes have reads in more than 5% of total cells)
No GPU device has been detected. Switch to CPU mode.

Loss (Total): 593.525, (Spliced): 290.466, (Unspliced): 303.059: 100%|██████████████████████████████████████████████████████████████████████████████████████████▉| 11999/12000 [1:25:32<00:00,  2.25it/s]

Loss (Total): 593.525, (Spliced): 290.466, (Unspliced): 303.059: 100%|██████████████████████████████████████████████████████████████████████████████████████████▉| 11999/12000 [1:26:54<00:00,  2.30it/s]

Total loss 593.525, vgene loss 593.525


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 adata = utv.run_model(adata, 
      2                       label, config_file=velo_config)

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/main.py:64, in run_model(adata, label, config_file, normalize)
     61 model = Velocity(adata, config=config)
     62 model.get_velo_genes()
---> 64 adata = model.fit_velo_genes(basis=basis, rep=rep)
     65 pre_time_gm = adata.obs['latent_time'].values
     67 if config.GENERAL != 'Linear':

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/velocity.py:268, in Velocity.fit_velo_genes(self, basis, rep)
    264 assert self.config.GENERAL in ['Deterministic', 'Curve', 'Linear'], \
    265     'Please specify the correct self.GENERAL in configuration file.'
    267 if self.config.GENERAL == 'Curve':
--> 268     residual, adata = self.fit_curve(self.adata, idx, Ms_scale, Mu_scale, rep=rep)
    270 if self.config.GENERAL == 'Deterministic':
    271     residual = self.fit_deterministic(idx, self.Ms, self.Mu, Ms_scale, Mu_scale)

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/velocity.py:246, in Velocity.fit_curve(self, adata, idx, Ms_scale, Mu_scale, rep)
    243     print (f'Using GPU card: {self.config.GPU}')
    245 with tf.device(device):
--> 246     residual, adata = lagrange(
    247         adata, idx=idx,
    248         Ms=Ms_scale, Mu=Mu_scale, 
    249         rep=rep, config=self.config
    250     )
    252 return residual, adata

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/model.py:416, in lagrange(adata, idx, Ms, Mu, var_names, rep, config)
    404 var_names = make_unique_list(var_names, allow_array=True)
    406 model = Recover_Paras(
    407     adata,
    408     Ms,
   (...)
    413     config=config
    414 )
--> 416 latent_time_gm, s_derivative, adata = model.fit_likelihood()
    418 if 'latent_time' in adata.obs.columns:
    419     del adata.obs['latent_time']

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/model.py:340, in Recover_Paras.fit_likelihood(self)
    337 self.adata.var['velocity_genes'] = self.total_genes if not self.flag else self.idx
    338 self.adata.layers['fit_t'][:, ~self.adata.var['velocity_genes'].values] = np.nan
--> 340 return self.get_interim_t(t_cell, self.adata.var['velocity_genes'].values), s_derivative.numpy(), self.adata

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/optimize_utils.py:351, in Model_Utils.get_interim_t(self, t_cell, idx)
    347         temp = np.reshape(t_cell[:, i], (-1, 1))
    348         t_interim[:, i] = \
    349             np.squeeze(col_minmax(temp, self.adata.var.index[i]))
--> 351 t_interim = self.max_density(t_interim)
    352 t_interim = tf.reshape(t_interim, (-1, 1))
    353 t_interim = tf.broadcast_to(t_interim, self.adata.shape)

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/optimize_utils.py:180, in Model_Utils.max_density(self, dis)
    177     return tf.cast(mean(dis_approx, axis=1), tf.float32)
    179 if self.config.DENSITY == 'Raw':
--> 180     weight = self.gene_prior_perc(dis)
    181     return tf.cast(sum(dis * weight, axis=1), tf.float32)

File ~/miniconda3/lib/python3.9/site-packages/unitvelo/optimize_utils.py:187, in Model_Utils.gene_prior_perc(self, dis)
    184 perc, perc_idx = [], []
    185 aggregrate_weight = np.ones((dis.shape[1], ), dtype=np.float32)
--> 187 for prior in self.config.GENE_PRIOR:
    188     expr = np.array(self.adata[:, prior[0]].layers['Ms'])
    189     perc.append(np.max(expr) * 0.75) # modify 0.75 for parameter tuning

TypeError: 'NoneType' object is not iterable

preprocessing script for input to `SuppFig6_Oligodendrocyte.ipynb`

Would you be able to provide the pre-processing script used to produce

11 "dataset = '../data/Oligo/oligo_lite.h5ad'

presumably from OligoPons_metafromSten.loom or oligos_info.rds that is
referenced at https://github.com/StatBiomed/UniTVelo/blame/8370e50e6cec261b6791f0e714736cb54f3c4088/notebooks/SuppFig6_Oligodendrocyte.ipynb#L11 ?

Many thanks!

Error: C stack usage is too close to the limit

Hello,

I recently upgraded UniTVelo from v0.2.3 to v0.2.5 and found that it didn't work. UniTVelo v0.2.3 was successfully run with the same data and codes.

The following is the log:

(Running UniTVelo 0.2.5)
2022-12-08 15:16:21
------> Manully Specified Parameters <------
NUM_REP:        2
NUM_REP_TIME:   re_init
BASIS:  pca
------> Model Configuration Settings <------
N_TOP_GENES:    2000
LEARNING_RATE:  0.01
FIT_OPTION:     1
DENSITY:        SVD
REORDER_CELL:   Soft_Reorder
AGGREGATE_T:    True
R2_ADJUST:      True
GENE_PRIOR:     None
VGENES: basic
IROOT:  None
--------------------------------------------

Current working dir is /home/XXX.
Results will be stored in res folder
Filtered out 54980 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.
Extracted 2000 highly variable genes.
Computing moments for 2000 genes with n_neighbors: 30 and n_pcs: 30
computing neighbors
    finished (0:00:08) --> 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)

# of velocity genes 982 (Criterion: positive regression coefficient between un/spliced counts)
# of velocity genes 953 (Criterion: std of un/spliced reads should be moderate, w/o extreme values)
# of velocity genes 949 (Criterion: genes have reads in more than 5% of total cells)
Using GPU card: 0
Loss (Total): 58018.953, (Spliced): 30779.621, (Unspliced): 27239.332:  38%|██████████████████▌                              | 4545/12000 [11:18<15:35,  7.97it/s]Build failed. See log files for full details.
Error: C stack usage  7970020 is too close to the limit

cluster specific velocity genes

KeyError: 'highly_variable'

Hello dev team!
Thank you so much for developing such an amazing tool!

I'm trying to run unitvelo and got an error about "highly_variable". I checked the adata.var and there are highly variable genes there. Could you help figure out how to solve this issue?

What I did a step before is to integrate adata obtained from Seurat convertion with loom files with unspliced and spliced:

These are the variables in the adata

adata
AnnData object with n_obs × n_vars = 329 × 1807
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.9', 'seurat_clusters', 'RNA_snn_res.1', 'RNA_snn_res.1.2', 'integrated_snn_res.0.8', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size'
    var: 'features', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    uns: 'neighbors'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
    obsp: 'distances'

And I added the highly variable through the command:
sc.pp.highly_variable_genes(adata)

adata
AnnData object with n_obs × n_vars = 329 × 1442
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.9', 'seurat_clusters', 'RNA_snn_res.1', 'RNA_snn_res.1.2', 'integrated_snn_res.0.8', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'features', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'neighbors', 'temp', 'hvg'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
    obsp: 'distances'

When I run the UniTVelo command, I get the following error

adata = utv.run_model("results/velocity_integrated.h5ad", label="seurat_clusters")
Model configuration file not specified. Default settings with unified-time mode will be used.
------> Manully Specified Parameters <------
------> Model Configuration Settings <------
N_TOP_GENES:	2000
LEARNING_RATE:	0.01
FIT_OPTION:	1
DENSITY:	SVD
REORDER_CELL:	Soft_Reorder
AGGREGATE_T:	True
R2_ADJUST:	True
GENE_PRIOR:	None
VGENES:	basic
IROOT:	None
--------------------------------------------

Filtered out 365 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.
Skip filtering by dispersion since number of variables are less than `n_top_genes`.
WARNING: Did not modify X as it looks preprocessed already.
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   3360             try:
-> 3361                 return self._engine.get_loc(casted_key)
   3362             except KeyError as err:

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'highly_variable'

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

KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_1863459/1623242212.py in <module>
----> 1 adata = utv.run_model("results/velocity_integrated.h5ad", label="seurat_clusters")
      2 #scv.pl.velocity_embedding_stream(adata)

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/unitvelo/main.py in run_model(adata, label, config_file, normalize)
     27     from .utils import init_config_summary, init_adata_and_logs
     28     config, _ = init_config_summary(config=config_file)
---> 29     adata, data_path = init_adata_and_logs(adata, config, normalize=normalize)
     30 
     31     scv.settings.presenter_view = True

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/unitvelo/utils.py in init_adata_and_logs(adata, config, normalize)
    605                                     min_shared_counts=config.MIN_SHARED_COUNTS,
    606                                     n_top_genes=config.N_TOP_GENES)
--> 607         print (f"Extracted {adata.var[adata.var['highly_variable'] == True].shape[0]} highly variable genes.")
    608 
    609         print (f'Computing moments for {len(adata.var)} genes with n_neighbors: {config.N_NEIGHBORS} and n_pcs: {config.N_PCS}')

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/pandas/core/frame.py in __getitem__(self, key)
   3456             if self.columns.nlevels > 1:
   3457                 return self._getitem_multilevel(key)
-> 3458             indexer = self.columns.get_loc(key)
   3459             if is_integer(indexer):
   3460                 indexer = [indexer]

~/miniconda3/envs/unitvelo/lib/python3.7/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   3361                 return self._engine.get_loc(casted_key)
   3362             except KeyError as err:
-> 3363                 raise KeyError(key) from err
   3364 
   3365         if is_scalar(key) and isna(key) and not self.hasnans:

KeyError: 'highly_variable'

I appreciate your support in advance.

TypeError: Tensor is unhashable.

Hi, I'm trying to run the example as https://unitvelo.readthedocs.io/en/latest/Figure2_ErythroidMouse.html shows. When I run:

adata = utv.run_model('./data/Gastrulation/erythroid_lineage.h5ad', label, config_file=velo_config)

I received the error:

Traceback (most recent call last):
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/arrays/categorical.py", line 431, in init
codes, categories = factorize(values, sort=True)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/algorithms.py", line 760, in factorize
codes, uniques = factorize_array(
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/algorithms.py", line 562, in factorize_array
uniques, codes = table.factorize(
File "pandas/_libs/hashtable_class_helper.pxi", line 5396, in pandas._libs.hashtable.PyObjectHashTable.factorize
File "pandas/_libs/hashtable_class_helper.pxi", line 5310, in pandas._libs.hashtable.PyObjectHashTable._unique
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/tensorflow/python/framework/ops.py", line 897, in hash
raise TypeError("Tensor is unhashable. "
TypeError: Tensor is unhashable. Instead, use tensor.ref() as the key.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "", line 1, in
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/unitvelo/main.py", line 64, in run_model
adata = model.fit_velo_genes(basis=basis, rep=rep)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/unitvelo/velocity.py", line 284, in fit_velo_genes
scv.tl.latent_time(adata, min_likelihood=None)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/scvelo/tools/dynamical_model.py", line 804, in latent_time
terminal_states(adata, vkey=vkey)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/scvelo/tools/terminal_states.py", line 275, in terminal_states
strings_to_categoricals(adata)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/scvelo/tools/utils.py", line 304, in strings_to_categoricals
c = Categorical(c)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/arrays/categorical.py", line 433, in init
codes, categories = factorize(values, sort=False)
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/algorithms.py", line 760, in factorize
codes, uniques = factorize_array(
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/pandas/core/algorithms.py", line 562, in factorize_array
uniques, codes = table.factorize(
File "pandas/_libs/hashtable_class_helper.pxi", line 5396, in pandas._libs.hashtable.PyObjectHashTable.factorize
File "pandas/_libs/hashtable_class_helper.pxi", line 5310, in pandas._libs.hashtable.PyObjectHashTable._unique
File "/home/linqianyu/anaconda3/envs/UniTVelo/lib/python3.9/site-packages/tensorflow/python/framework/ops.py", line 897, in hash
raise TypeError("Tensor is unhashable. "
TypeError: Tensor is unhashable. Instead, use tensor.ref() as the key.

All the parameters are same with the example, same error also occur when I tried to run my own data.

'KeyError: 'The optimizer cannot recognize variable log_gamma:0.

Hi,

Thank you a lot for developing this great tool! I got the following error when I run the Mouse Erythroid dataset under UniTVelo mode. No problem under scVelo stochastic or scVelo dynamic mode.

'KeyError: 'The optimizer cannot recognize variable log_gamma:0. This usually means you are trying to call the optimizer to update different parts of the model separately. Please call optimizer.build(variables) with the full list of trainable variables before the training loop or use legacy optimizer `tf.keras.optimizers.legacy.{self.class.name}.''

Best