12.1 前言

單細(xì)胞數(shù)據(jù)集允許以高分辨率研究生物過程，例如早期發(fā)育鹅心。例如，雖然分析的是單個(gè)細(xì)胞而不是整個(gè)組織旭愧，但細(xì)胞表型的變化無法隨著時(shí)間的推移進(jìn)行跟蹤虐秋。這一事實(shí)源于單細(xì)胞測序方案的局限性。對細(xì)胞進(jìn)行測序后客给，它會(huì)被破壞，因此無法在以后的時(shí)間點(diǎn)再次測量其定義特征靶剑。值得注意的是，實(shí)驗(yàn)技術(shù)不僅無法測量不同時(shí)間的一般細(xì)胞概況缎讼，而且無法測量這些變化發(fā)生的速度坑匠。可以使用軌跡推理(trajectory inference厘灼，TI)領(lǐng)域的工具來及時(shí)恢復(fù)沿發(fā)展景觀的位置。然而设凹，經(jīng)典的TI方法不提供任何定向的動(dòng)態(tài)信息。此外月匣，這些算法傳統(tǒng)上不考慮轉(zhuǎn)錄組讀數(shù)和相似性之外的信息钻洒。

12.2 RNA velocity模型

細(xì)胞轉(zhuǎn)錄組譜的變化是由一系列事件觸發(fā)的：廣義上講素标，DNA 被轉(zhuǎn)錄產(chǎn)生所謂的未剪接前體信使 RNA (pre-mRNA)院刁。未剪接的前mRNA包含與翻譯相關(guān)的區(qū)域（外顯子）以及非編碼區(qū)域（內(nèi)含子）。這些非編碼區(qū)被剪接退腥，即去除，以形成剪接的成熟mRNA狡刘。雖然單細(xì)胞RNA測序 (scRNA-seq) 方案無法在多個(gè)時(shí)間點(diǎn)捕獲轉(zhuǎn)錄組，但它們確實(shí)包含了分離未剪接和剪接mRNA讀數(shù)所需的信息剑按。

12.3 胰腺內(nèi)分泌發(fā)生pancreatic endocrinogenesis中的RNA velocity推斷

作為如何推斷RNA velocity的實(shí)際示例，我們分析了胰腺的內(nèi)分泌發(fā)育[Bastidas-Ponce et al., 2019]艺蝴。在該系統(tǒng)中鸟废，前內(nèi)分泌細(xì)胞（Ductal、Ngn3 low EP盒延、Ngn3 high EP、Pre-endocrine）發(fā)育成四種內(nèi)分泌細(xì)胞類型（Alpha胯盯、Beta、Delta博脑、Epsilon）票罐。在這里，我們使用scVelo[Bergen et al., 2020]來推斷RNA velocity胶坠。

12.3.1 配置環(huán)境

from pathlib import Path

import scanpy as sc
import scvelo as scv

12.3.2 常規(guī)設(shè)置

scv.settings.set_figure_params("scvelo")
DATA_DIR = Path("../../data/")
DATA_DIR.mkdir(parents=True, exist_ok=True)

FILE_PATH = DATA_DIR / "pancreas.h5ad"

12.3.3 加載數(shù)據(jù)

為了使用scVelo估計(jì)RNA velocity繁堡，需要將未剪接和剪接計(jì)數(shù)存儲(chǔ)在AnnData的層槽中乡数。我們建議將整個(gè)計(jì)數(shù)（即未處理的數(shù)據(jù)）傳遞到scVelo管道闻牡。

adata = scv.datasets.pancreas(file_path=FILE_PATH)
adata

AnnData object with n_obs × n_vars = 3696 × 27998
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
    var: 'highly_variable_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'distances', 'connectivities'

12.4 數(shù)據(jù)預(yù)處理

由于scRNA-seq數(shù)據(jù)充滿噪聲且稀疏，因此必須對數(shù)據(jù)進(jìn)行預(yù)處理玖翅，以便使用穩(wěn)態(tài)或EM模型推斷RNA velocity。第一步金度，我們過濾掉未剪接和剪接RNA均未充分表達(dá)的基因（此處至少20個(gè)）严沥。接下來，對未剪接和剪接RNA的細(xì)胞大小進(jìn)行歸一化消玄，并在adata.X log1p中進(jìn)行計(jì)數(shù)轉(zhuǎn)換，以減少異常值的影響受扳。接下來，我們還識(shí)別和過濾高度可變的基因（此處為2000）勘高。

scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)

Filtered out 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Extracted 2000 highly variable genes.
Logarithmized X.

到目前為止浮定，數(shù)據(jù)預(yù)處理與經(jīng)典的scRNA-seq工作流程類似。對于RNA velocity桦卒，我們還通過其附近的平均表達(dá)來觀察結(jié)果。 這可以使用scVelo的moments函數(shù)來完成方灾。

sc.tl.pca(adata)
sc.pp.neighbors(adata)
scv.pp.moments(adata, n_pcs=None, n_neighbors=None)

computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)

在經(jīng)典的工作流程中，我們會(huì)對數(shù)據(jù)進(jìn)行聚類洞慎，推斷細(xì)胞類型嘿棘，并以二維嵌入方式可視化數(shù)據(jù)。幸運(yùn)的是鸟妙，對于胰腺數(shù)據(jù)挥吵，這些信息已經(jīng)被先驗(yàn)計(jì)算并可以直接使用花椭。

scv.pl.scatter(adata, basis="umap", color="clusters")

12.4.1 RNA velocity推斷-穩(wěn)態(tài)模型

第一步，我們計(jì)算穩(wěn)態(tài)模型下的RNA velocity矿辽。在這種情況下，我們用mode="deterministic"來調(diào)用scVelo的velocity函數(shù)雕蔽。

scv.tl.velocity(adata, mode="deterministic")

computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)

雖然我們不鼓勵(lì)過度解釋高維速度向量到數(shù)據(jù)低維表示的投影，但scVelo提供了一種簡單的方法萎羔。

scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

computing velocity graph (using 8/96 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)

12.4.2 RNA velocity推斷-EM模型

為了用EM模型計(jì)算RNA velocity贾陷，需要首先推斷剪接動(dòng)力學(xué)參數(shù)。scVelo的recover_dynamics函數(shù)負(fù)責(zé)完成髓废。

scv.tl.recover_dynamics(adata, n_jobs=8)

recovering dynamics (using 8/96 cores)

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

拼接模型的參數(shù)是通過最大化給定的可能性來推斷的该抒。為了研究哪些基因最適合scVelo，我們可以研究相應(yīng)的相圖以及推斷的軌跡（以紫色繪制）和穩(wěn)態(tài)比率（紫色虛線）凑保。此處，顯示的五個(gè)基因中的三個(gè)（Pcsk2欧引、Top2a、Ppp1r1a）呈現(xiàn)出（部分）杏仁形狀的相圖芝此。我們觀察到單一細(xì)胞類型（Top2a、Ppp1r1a）內(nèi)或跨多種細(xì)胞類型（Pcsk2岸更，從前內(nèi)分泌細(xì)胞到 Alpha 和 Beta）的明顯轉(zhuǎn)變。就Nfib而言怎炊，我們觀察到兩個(gè)處于穩(wěn)態(tài)的細(xì)胞群。這很可能是對Ngn3低/高EP細(xì)胞周圍表型流形進(jìn)行欠采樣的人為因素评肆。類似地，Ghrl在Epsilon細(xì)胞中高度表達(dá)，但由于簇大小較小院仿，表達(dá)量很少。雖然當(dāng)前的最佳實(shí)踐僅限于手動(dòng)分析模型擬合度和置信度歹垫，但最近提出的方法可以幫助自動(dòng)化該過程（新方向）。 在這里排惨，Nfib和Ghrl 將被分配較低的置信度分?jǐn)?shù)。

top_genes = adata.var["fit_likelihood"].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[:5], color="clusters", frameon=False)

scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap")

computing velocities
    finished (0:00:02) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 8/96 cores)
    finished (0:00:02) --> 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)

基于2D投影凄鼻，EM模型更穩(wěn)健地捕獲導(dǎo)管細(xì)胞中的細(xì)胞周期聚假。此外，穩(wěn)態(tài)模型的投影表現(xiàn)出從阿爾法細(xì)胞到前內(nèi)分泌細(xì)胞的“回流”膘格。然而，為了進(jìn)行嚴(yán)格的定量分析瘪贱，我們建議使用CellRank等下游工具來評估模型差異并得出結(jié)論。