theislab/scvelo

scv.tl.recover_dynamics() IndexError: list index out of range

AAA-3 opened this issue · 2 comments

Hello all,

Thanks for this really nice tool - I am having trouble with my analysis and am hoping someone can tell me what is going on here. I am using hte same pipeline as I did a couple of years ago on single nuclei sequencing.

I have processed files from 10X through velocyto, concatenated the multiple loom files, and subset the samples to include only cells that passed my filtering in seurat. I imported the PCA etc. from the seurat metadata into the loom file.

I got an error at the scv.tl.recover_dynamics() steps.

Tried with untouched velocyto output, and got the same error, but scv.datasets.pancreas() works. I can't figure out what I should change since this is simply a velocyto loom...

Can anyone tell me what is going on?

import scanpy as sc
import anndata
import scvelo as scv
import pandas as pd
pd.options.display.max_columns = None
import numpy as np
import matplotlib.pyplot as plt
import rpy2
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
import umap
%matplotlib inline

H9 = scv.read("/home/ali/Dokumente/15d _Organoids_Multiome/Separate/Velocyto/Results/All/H9_raw.h5ad")

scv.pp.filter_and_normalize(H9)
scv.pp.moments(H9)

#here is where I get the error
scv.tl.recover_dynamics(H9, plot_results=True, n_jobs=30, n_top_genes=2000)
Error output
recovering dynamics (using 30/40 cores)

100%
 375/375 [05:21<00:00,  7.00s/gene]

/home/ali/anaconda3/envs/scvelo/lib/python3.12/multiprocessing/popen_fork.py:66: DeprecationWarning: This process (pid=4504) is multi-threaded, use of fork() may lead to deadlocks in the child.
  self.pid = os.fork()
/home/ali/anaconda3/envs/scvelo/lib/python3.12/site-packages/joblib/externals/loky/backend/fork_exec.py:38: DeprecationWarning: This process (pid=4504) is multi-threaded, use of fork() may lead to deadlocks in the child.
  pid = os.fork()

    finished (0:15:07) --> added 
    'fit_pars', fitted parameters for splicing dynamics (adata.var)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[14], line 2
      1 #velocity analysis
----> 2 scv.tl.recover_dynamics(H9, plot_results=True, n_jobs=30, n_top_genes=2000)
      3 scv.tl.velocity(H9, mode='dynamical', groupby="Genotype")
      4 scv.tl.velocity_graph(H9, n_jobs=30)

File ~/anaconda3/envs/scvelo/lib/python3.12/site-packages/scvelo/tools/_em_model_core.py:653, in recover_dynamics(data, var_names, n_top_genes, max_iter, assignment_mode, t_max, fit_time, fit_scaling, fit_steady_states, fit_connected_states, fit_basal_transcription, use_raw, load_pars, return_model, plot_results, steady_state_prior, add_key, copy, n_jobs, backend, show_progress_bar, **kwargs)
    651 if t_max is not False:
    652     mi = dm.m[var_id]
--> 653     P[var_id] *= np.array([1 / mi, 1 / mi, 1 / mi, mi, 1])[:, None]
    654 ax = axes[var_id] if n_rows > 1 else axes
    655 for j, pij in enumerate(P[var_id]):

IndexError: list index out of range

output of the velocyto loom H9

AnnData object with n_obs × n_vars = 71269 × 36601
    obs: 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'gene_count_corr', 'fit_r2', 'fit_alpha', 'fit_beta', 'fit_gamma', 'fit_t_', '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'
    uns: 'log1p', 'pca', 'neighbors', 'recover_dynamics'
    obsm: 'X_pca'
    varm: 'PCs', 'loss'
    layers: 'ambiguous', 'matrix', 'spliced', 'unspliced', 'Ms', 'Mu', 'fit_t', 'fit_tau', 'fit_tau_'
    obsp: 'distances', 'connectivities'

Hello!
Somehow running scv.tl.recover_dynamics(H9) seemed to work....