In [1]:
import numpy as np
import matplotlib.pyplot as pl

import scvelo as scv
scv.logging.print_version()
Running scvelo 0.1.24 (python 3.7.3) on 2019-10-28 18:26.
In [2]:
scv.settings.set_figure_params('scvelo', dpi_save=200, dpi=80, transparent=True)  # vectorized: pdf or svg
scv.settings.verbosity = 2

Capturing non-stationary populations

In [3]:
n_vars = 2000

mu = np.array([5, .2, .05])

R = np.array([[1., .2, .2], 
              [.2, 1., .8],
              [.2, .8, 1.]])

C = np.array([.4, .4, .4])[:, None]

cov = C.dot(C.T) * R

alpha, beta, gamma = np.exp(np.random.multivariate_normal(mu, cov, size=n_vars).T)  # multivariate log-normal
beta /= 3
gamma /= 3

# remove outliers
idx = (alpha < np.percentile(alpha, 99)) & (beta < np.percentile(beta, 99)) & (gamma < np.percentile(gamma, 99))
alpha = alpha[idx]
beta = beta[idx]
gamma = gamma[idx]
n_vars = np.sum(idx)

switches = np.random.uniform(.1, .5, size=n_vars)

adata = scv.datasets.simulation(n_obs=1000, t_max=20, n_vars=n_vars, switches=switches,
                                alpha=alpha, beta=beta, gamma=gamma, noise_level=.8)
# adata.var['switches'] = switches

scv.pp.neighbors(adata)
scv.tl.velocity(adata, mode='steady_state', use_raw=True)

var_names = adata.var_names[:4]
scv.tl.recover_dynamics(adata, var_names=var_names, max_iter=100)

scv.pl.scatter(adata, basis=var_names)
computing neighbors
    finished (0:00:04)
computing velocities
    finished (0:00:00)
recovering dynamics
    finished (0:00:01)

outputs model fit of gene: 3
In [4]:
# adata.write('data/toy/Fig1_identifiability.h5ad')
# adata = scv.read('data/toy/Fig1_identifiability.h5ad')
In [5]:
dm = scv.tl.recover_dynamics(adata, var_names='all', max_iter=100, use_raw=True)

var_names = adata.var_names[:4]
scv.pl.scatter(adata, basis=var_names, use_raw=True)