3. scVelo: Differential Kinetics

https://scvelo.readthedocs.io/DifferentialKinetics.html

最終更新日: 2023/10/1

重要な懸念事項の一つは、複数の系統やプロセスを表す系を扱うことで、遺伝子が亜集団間で異なる状態を示す可能性があることです。 異なる細胞状態や系統は通常異なる遺伝子制御ネットワークによって支配されており、それゆえ異なるスプライシング動態を示す可能性があります。これにより、位相空間で複数の軌跡を示す遺伝子が発生することになります。

この問題に対処するために、力学モデルを使用して、キネティクスの差分の尤度比検定を実行することができます。こうすることで、全体的なダイナミクスの単一モデルではうまく説明できないような挙動を示すクラスタを検出することができます。細胞種ごとに異なるキネティクスレジームにクラスタリングすることで、各レジームを個別に適合させることができます。

説明のために、歯状回神経新生に微分動力学解析を適用してみる。 例として、複数の異質な部分集団から構成されている歯状回神経新生 (dentate gyrus neurogenesis) に差分キネティクス解析を適用してみましょう。

import scvelo as scv
scv.logging.print_version()
scv.settings.verbosity = 3  # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True  # set max width size for presenter view
scv.settings.set_figure_params('scvelo')  # for beautified visualization

Running scvelo 0.2.5 (python 3.8.16) on 2023-10-01 19:08.

ERROR: XMLRPC request failed [code: -32500]
RuntimeError: PyPI no longer supports 'pip search' (or XML-RPC search). Please use https://pypi.org/search (via a browser) instead. See https://warehouse.pypa.io/api-reference/xml-rpc.html#deprecated-methods for more information.

3.1. データの準備

前準備として、HVGの選択、総サイズによる遺伝子発現量の正規化・対数化、速度推定のためのモーメント計算を実行します。 詳しい説明は前のチュートリアルを参照。

adata = scv.datasets.dentategyrus()
scv.pp.filter_and_normalize(adata, min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

Filtered out 11019 genes that are detected 30 counts (shared).
Normalized count data: X, spliced, unspliced.
Extracted 2000 highly variable genes.
Logarithmized X.
computing neighbors
    finished (0:00:05) --> added
    'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
    finished (0:00:00) --> added
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)

3.2. Basic Velocity Estimation

scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis='umap')

computing velocities
    finished (0:00:00) --> added
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/40 cores)

    finished (0:00:05) --> added
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)

3.3. Differential Kinetic Test

異なる細胞種や系統は異なる遺伝子ネットワークによって支配され得るので、異なるキネティクスレジームを示すことがあります。細胞種や系統が関連していても、alternative splicing、alternative polyadenylation、分解における調節により、キネティクスが異なる場合があります。

力学モデルによって、全体的なキネティクスを表す単一のモデルでは十分に説明できないような動態を示すクラスターを検出するために、キネティクスの差異に関する尤度比検定でこの問題に対処することができる。各細胞種は独立したフィットが有意に改善された尤度をもたらすかどうかについてテストされる。

尤度比は漸近的カイ二乗分布に従い、有意性を検定することができる。効率上の理由から、デフォルトでは、クラスタが全体的な力学によって十分に説明されているか、または異なる力学を示すかをテストするために、完全な位相軌跡の代わりに直交回帰が使用されていることに注意。

var_names = ['Tmsb10', 'Fam155a', 'Hn1', 'Rpl6']
scv.tl.differential_kinetic_test(adata, var_names=var_names, groupby='clusters')

recovering dynamics (using 1/40 cores)
    finished (0:00:02) --> added
    'fit_pars', fitted parameters for splicing dynamics (adata.var)

outputs model fit of gene: Rpl6
testing for differential kinetics
    finished (0:00:00) --> added
    'fit_diff_kinetics', clusters displaying differential kinetics (adata.var)
    'fit_pvals_kinetics', p-values of differential kinetics (adata.var)

outputs model fit of gene: Rpl6

<scvelo.tools.dynamical_model.DynamicsRecovery at 0x7f4ab5d0b160>

scv.get_df(adata[:, var_names], ['fit_diff_kinetics', 'fit_pval_kinetics'], precision=2)

	fit_diff_kinetics	fit_pval_kinetics
index
Tmsb10	Endothelial	6.02e-16
Fam155a	Cajal Retzius	8.35e-161
Hn1	Microglia	3.02e-03
Rpl6	Microglia	5.69e-16

kwargs = dict(linewidth=2, add_linfit=True, frameon=False)
scv.pl.scatter(adata, basis=var_names, add_outline='fit_diff_kinetics', **kwargs)

例えばTmsb10では、内皮は全体的なダイナミクス（紫色の曲線）ではうまく説明できないような動力学的挙動を示す（縁取られた◯、黒線がフィットされた線）。

diff_clusters=list(adata[:, var_names].var['fit_diff_kinetics'])
scv.pl.scatter(adata, legend_loc='right', size=60, title='diff kinetics',
               add_outline=diff_clusters, outline_width=(.8, .2))

3.4. Testing top-likelihood genes

尤度の高い遺伝子をスクリーニングすると、複数のキネティックレジームを示す遺伝子単位のダイナミクスがいくつか見つかります。

scv.tl.recover_dynamics(adata)
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:100]
scv.tl.differential_kinetic_test(adata, var_names=top_genes, groupby='clusters')

recovering dynamics (using 1/40 cores)

    finished (0:03:33) --> added
    'fit_pars', fitted parameters for splicing dynamics (adata.var)
testing for differential kinetics
    finished (0:00:20) --> added
    'fit_diff_kinetics', clusters displaying differential kinetics (adata.var)
    'fit_pvals_kinetics', p-values of differential kinetics (adata.var)

特に、Cck/Tox、GABA、内皮、ミクログリアなど、主顆粒とは異なる細胞タイプが頻繁に出現する。

scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, add_outline='fit_diff_kinetics', **kwargs)

scv.pl.scatter(adata, basis=top_genes[15:30], ncols=5, add_outline='fit_diff_kinetics', **kwargs)

3.5. Recompute velocities

最終的に、複数の競合するキネティクスレジームの情報を利用して速度を再計算することができます。

scv.tl.velocity(adata, diff_kinetics=True)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding(adata, dpi=120, arrow_size=2, arrow_length=2)

computing velocities
    finished (0:00:02) --> added
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/40 cores)

    finished (0:00:05) --> added
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)

import session_info
session_info.show()

Click to view session information

-----
anndata             0.9.1
pandas              1.5.3
scvelo              0.2.5
session_info        1.0.0
-----

Click to view modules imported as dependencies

PIL                         9.4.0
asciitree                   NA
asttokens                   NA
backcall                    0.2.0
beta_ufunc                  NA
binom_ufunc                 NA
bottleneck                  1.3.5
cffi                        1.15.1
cloudpickle                 2.2.1
colorama                    0.4.6
comm                        0.1.2
cycler                      0.10.0
cython_runtime              NA
dask                        2023.5.0
dateutil                    2.8.2
debugpy                     1.5.1
decorator                   5.1.1
defusedxml                  0.7.1
entrypoints                 0.4
executing                   0.8.3
fasteners                   0.18
h5py                        3.8.0
hypergeom_ufunc             NA
igraph                      0.10.4
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.19.2
ipython_genutils            0.2.0
ipywidgets                  8.0.4
jedi                        0.18.1
jinja2                      3.1.2
joblib                      1.2.0
jupyter_server              1.23.4
kiwisolver                  1.4.4
leidenalg                   0.9.1
llvmlite                    0.39.1
louvain                     0.8.0
markupsafe                  2.1.1
matplotlib                  3.7.1
matplotlib_inline           0.1.6
mpl_toolkits                NA
msgpack                     1.0.5
natsort                     8.3.1
nbinom_ufunc                NA
ncf_ufunc                   NA
numba                       0.56.4
numcodecs                   0.11.0
numexpr                     2.8.4
numpy                       1.23.5
packaging                   23.0
parso                       0.8.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
platformdirs                2.5.2
plotly                      5.14.1
prompt_toolkit              3.0.36
psutil                      5.9.0
ptyprocess                  0.7.0
pure_eval                   0.2.2
pycparser                   2.21
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.6.0
pydevd_concurrency_analyser NA
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.15.1
pynndescent                 0.5.10
pyparsing                   3.0.9
pytz                        2022.7
scanpy                      1.9.3
scipy                       1.9.1
setuptools                  66.0.0
six                         1.16.0
sklearn                     1.2.2
stack_data                  0.2.0
texttable                   1.6.7
threadpoolctl               3.1.0
tlz                         0.12.0
toolz                       0.12.0
tornado                     6.2
tqdm                        4.65.0
traitlets                   5.7.1
typing_extensions           NA
umap                        0.5.3
wcwidth                     0.2.5
yaml                        6.0
zarr                        2.14.2
zipp                        NA
zmq                         25.0.2
zope                        NA

-----
IPython             8.12.0
jupyter_client      8.1.0
jupyter_core        5.3.0
jupyterlab          3.5.3
notebook            6.5.4
-----
Python 3.8.16 (default, Mar  2 2023, 03:21:46) [GCC 11.2.0]
Linux-5.15.0-79-generic-x86_64-with-glibc2.17
-----
Session information updated at 2023-10-01 19:13