RNAvelocity6：scVelo用于RNA 速率基本流程
分享是一种态度
在这里，您将学习RNA速率分析的基本流程。
示例数据采用胰腺发育数据集，其谱系包含四个主要细胞命运：α、β、δ和ε细胞。有关详细信息，请参阅此处[1]。可以按照相同的思路应用于您自己的数据。
[ ]:
# update to the latest version, if not done yet.
!pip install scvelo --upgrade --quiet
[1]:
import scvelo as scv
scv.logging.print_version()
Running scvelo 0.2.0 (python 3.8.2) on 2020-05-16 12:55.
[2]:
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.set_figure_params('scvelo') # for beautified visualization
加载数据
分析基于内置胰腺数据[2]。
要对自己的数据进行速率分析，请通过adata = scv.read('path/file.loom', cache=True)将您的文件（loom, h5ad, csv …）读取到AnnData数据对象。
如果您想将loom文件合并到已存在的 AnnData 对象中，请使用scv.utils.merge(adata, adata_loom)
[3]:
adata = scv.datasets.pancreas()
adata
[3]:
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'
scVelo 是基于adata，存储了数据矩阵adata.X、观测注释adata.obs、变量adata.var和非结构化注释的对象adata.uns。观测和变量的名称可分别通过adata.obs_names和adata.var_names访问。例如，adata对象可以像 adata_subset = adata[:,list_of_gene_names] 数据帧一样切片。
有关详细信息，请参阅anndata docs[3] 。
[4]:
scv.pl.proportions(adata)
这里，展示了剪切/未剪切count的比例。根据使用的协议（Drop-Seq，Smart-Seq），我们通常有10%-25%的未剪切分子含有内含子序列。我们还建议您检查cluster水平的变化，以验证剪切效率的一致性。在这里，我们发现变化如预期的那样，在循环导管细胞中未剪切的比例略低，在许多基因开始转录的Ngn3高表达的和内分泌前细胞中的比例更高。
预处理数据
预处理包括通过检测（count数最少）和高变异性（分散）进行基因选择，按其总大小使每个细胞标准化。
所有这些都归入一个函数scv.pp.filter_and_normalize，主要运行如下：
scv.pp.filter_genes(adata, min_shared_counts=20)
scv.pp.normalize_per_cell(adata)
scv.pp.filter_genes_dispersion(adata, n_top_genes=2000)
scv.pp.log1p(adata)
此外，我们需要在PCA空间中最近的邻居之间计算的第一和第二顺序时刻，汇总在scv.pp.moments，它计算了scv.pp.pca和scv.pp.neighbors。确定速率估计需要第一个顺序，而随机估计也需要第二个顺序时刻。
[5]:
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
Filtered out 20801 genes that are detected in less than 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Logarithmized X.
computing neighbors
finished (0:00:03) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:00) --> added
'Ms' and 'Mu', moments of spliced/unspliced abundances (adata.layers)
进一步的预处理（如批次校正）可用于消除不必要的变异源。有关更多详细信息，请参阅最佳实践[4]。请注意，任何额外的预处理步骤仅影响count data X，不应用于剪切/未剪切计数。
估计RNA速率
速率是基因表达空间中的载体，代表单个细胞运动的方向和速率。速率通过对拼接动力学的转录动力学进行建模获得，无论是随机模型（默认）还是确定性模型（通过设置mode='deterministic'）。对于每个基因，预成熟（未剪切）和成熟（剪切）mRNA计数的稳定状态比，构成一个恒定的转录状态。然后，从此比率中获取速率作为残差。正速率表示基因被向上调节，这发生在细胞显示该基因的未剪切mRNA的丰度高于预期的稳定状态。相反，负速表示基因被降低调节。
完整的动力学模型的解决方案是通过设置mode='dynamical'获得的，这需要事先运行scv.tl.recover_dynamics(adata)。我们将在下一个教程中详细阐述动态模型。
[6]:
scv.tl.velocity(adata)
computing velocities
finished (0:00:01) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
计算的速率存储在计数矩阵adata.layers中。
然后，基因间速率的组合可用于估计单个细胞的未来状态。为了将速率投射到低维嵌入中，估计了细胞到细胞过渡的概率。即，对于每个速率矢量，我们发现可能的细胞过渡方向。过渡概率是使用潜在细胞到细胞的过渡和速率矢量之间的共生相关性计算的，并存储在表示速率图的矩阵中。由此产生的速率图具有维度n obs×n obs，并总结了通过速率矢量很好地解释的可能的细胞状态变化（对于运行时速率，也可以通过设置approx=True在减少的 PCA 空间上计算）。
[7]:
scv.tl.velocity_graph(adata)
computing velocity graph
finished (0:00:10) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
对于各种应用，速率图可以转换为过渡矩阵，通过应用高斯函数将余弦值相关性转换为实际的过渡概率。这需要设置scv.utils.get_transition_matrix
如前所述，它内部用于将速率投射到低维嵌入中，通过对与 scv.tl.velocity_embedding获得的平均概率过渡概率进行平均转换。此外，我们可以通过scv.tl.terminal_states沿着马尔代夫转移矩阵追踪细胞的起源和潜在命运，从而获得根细胞和终点的轨迹。
预测速率
最后，速率被投影到任何嵌入，通过basis指定，并以以下方式之一可视化：
细胞水平 scv.pl.velocity_embedding， 网格线scv.pl.velocity_embedding_grid 流线型scv.pl.velocity_embedding_stream
请注意，数据已预计算 UMAP 嵌入，并注释了群集。在应用到您自己的数据时，可以使用scv.tl.umap和scv.tl.louvain。有关详细信息，请参阅[scanpy tutorial][5]。此外，所有绘图功能都默认使用，并且，您可以相应地设置basis='umap'和color='clusters'
[8]:
scv.pl.velocity_embedding_stream(adata, basis='umap')
computing velocity embedding
finished (0:00:00) --> added
'velocity_umap', embedded velocity vectors (adata.obsm)
流线显示的速率矢量可对发育过程进行详细分析。它准确地勾画了导管细胞和内分泌祖细胞的循环。此外，它阐明了细胞的谱系命运、细胞周期回归和内分泌细胞分化状态。
我们在单细胞水平上获得的速率矢量的最详细分辨率，每个箭头显示单个细胞运动的方向和速率。这揭示了 Ngn3 细胞（黄色）的早期内分泌命运，以及近端α细胞（蓝色）和瞬态β细胞（绿色）之间的明显差异。
[9]:
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
解释速率
这可能是最重要的部分，因为我们建议用户不要将生物结论限制在预测速率上，而是通过图像来检查个体基因动力学，以了解特定基因如何支持推断方向。
请参阅此处的GIF[6]，了解如何解释剪切与未剪切的图像。基因活动是由转录调节的。特定基因的转录诱导导致（新转录的）前体未剪切 mRNA 增加，而相反，抑制或没有转录会导致未转录 mRNA 的减少。拼接的 mRNA 由未剪切的 mRNA 生成，并遵循相同的趋势，并具有时滞。时间是一个隐藏/潜在的变量。因此，需要从实际测量中推断出动态：相像中显示的剪切和未剪切的 mRNA。
现在，让我们来检查一些标记基因的图像，通过scv.pl.velocity(adata, gene_names)或scv.pl.scatter(adata, gene_names)可视化
[10]:
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
黑线对应于估计的"稳定状态"比率，即未剪切与剪切的 mRNA 丰度之比，后者处于恒定的转录状态。特定基因的RNA速率被确定为残留的，即观察与稳定状态线的偏差程度。正速率表示基因被向上调节，这发生在细胞显示该基因的未剪切mRNA的丰度高于预期的稳定状态。相反，负速表示基因被降低调节。
例如，Cpe解释了上调节 Ngn3（黄色）到前内分泌（橙色）到β细胞（绿色）的方向，而Adk则解释了向下调节的Ductal（深绿色）到 Ngn3（黄色）到剩余内分泌细胞的方向。
[11]:
scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'],
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
识别重要基因
我们需要一种系统的方法来识别基因，这可能有助于解释由此产生的向量场和推断的谱系。为此，我们可以测试哪些基因具有群体特异性微分速率表达，与剩余群相比，其比例要高得多/更低。模块scv.tl.rank_velocity_genes运行差异速率 t 测试，并生成每个集群的基因排名。可以设置阈值（例如min_corr），以限制对基因候选者选择的测试。
[12]:
scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
df.head()
#ranking velocity genes
# finished (0:00:02) --> added
# 'rank_velocity_genes', sorted scores by group ids (adata.uns)
# 'spearmans_score', spearmans correlation scores (adata.var)
[12]:
[13]:
kwargs = dict(frameon=False, size=10, linewidth=1.5,
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
scv.pl.scatter(adata, df['Ngn3 high EP'][:5], ylabel='Ngn3 high EP', **kwargs)
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs)
例如， Ptprs, Pclo, Pam, Abcc8, Gnas等基因支持从Ngn3 高表达的 EP（黄色）到内分泌前（橙色）到Beta（绿色）的方向性。
循环祖细胞的速率
RNA速率检测到的细胞周期，通过细胞周期分数（相标记基因平均表达水平的标准化分数）在生物学上得到了确认。
[14]:
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
calculating cell cycle phase
--> 'S_score' and 'G2M_score', scores of cell cycle phases (adata.obs)
对于循环导管细胞，我们可以筛选通过S和G2M相标记。前一个模块还计算了一个spearmans 相关分数，我们可以用它来对相标记基因进行排名/排序，然后显示其相像。
[15]:
s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)
s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).index
g2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).index
kwargs = dict(frameon=False, ylabel='cell cycle genes')
scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)
特别是Hells 和Top2a非常适合解释循环祖细胞的矢量场。Top2a 在 G2M阶段实际达到峰值前不久被分配了高速。在那里，负速然后完美匹配了接下来的下调。
[16]:
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
速率和连贯性
还有两个有用的统计数据：
分化的速率或速率由速率矢量的长度给出。 矢量场的一致性（即速率矢量如何与其邻近速率相关）提供了置信度的衡量标准。
[17]:
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
--> added 'velocity_length' (adata.obs)
--> added 'velocity_confidence' (adata.obs)
--> added 'velocity_confidence_transition' (adata.obs)
这些提供了细胞在哪里以较慢/更快的速率分化以及方向未确定的见解。
在集群层面，我们发现，细胞周期回归（Ngn3低表达的EP）后分化速率显著加快，在β细胞生产过程中保持步调，同时在阿尔法细胞生产过程中减速。
[18]:
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
[18]:
速率图和伪时间
我们可以可视化速率图，以描绘所有速率推断的细胞到细胞连接/过渡。它可以通过设置 threshold.限制在高概率转换。例如，该图指示来自早期和晚期内分泌细胞的 Epsilon 细胞产生的两个阶段。
[19]:
scv.pl.velocity_graph(adata, threshold=.1)
此外，该图可用于绘制来自特定细胞的后代/祖先。在这里，内分泌前细胞可以追溯到其潜在的命运。
[20]:
x, y = scv.utils.get_cell_transitions(adata, basis='umap', starting_cell=70)
ax = scv.pl.velocity_graph(adata, c='lightgrey', edge_width=.05, show=False)
ax = scv.pl.scatter(adata, x=x, y=y, s=120, c='ascending', cmap='gnuplot', ax=ax)
最后，根据速率图，可以计算速率伪时间。在从图形中推断出在根细胞上的分布后，它测量从根细胞开始沿着图形行走后到达细胞所需的平均步数。
与diffusion 伪时间分析相反，它含蓄地推断根细胞，并基于定向速率图，而不是基于相似性的diffusion kernel。
[21]:
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
computing terminal states
identified 2 regions of root cells and 1 region of end points
finished (0:00:00) --> added
'root_cells', root cells of Markov diffusion process (adata.obs)
'end_points', end points of Markov diffusion process (adata.obs)
PAGA速率图
PAGA[7]图形概念已作为轨迹推理的性能最高的方法进行基准测试。它提供了数据拓扑图的图形状地图，其加权边缘与两个群集之间的连接相对应。在这里，PAGA 通过速率推断方向性扩展。
[ ]:
# PAGA requires to install igraph, if not done yet.
!pip install python-igraph --upgrade --quiet
[22]:
# this is needed due to a current bug - bugfix is coming soon.
adata.uns['neighbors']['distances'] = adata.obsp['distances']
adata.uns['neighbors']['connectivities'] = adata.obsp['connectivities']
scv.tl.paga(adata, groups='clusters')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
running PAGA
finished (0:00:01) --> added
'paga/transitions_confidence', connectivities adjacency (adata.uns)
'paga/connectivities', connectivities adjacency (adata.uns)
'paga/connectivities_tree', connectivities subtree (adata.uns)
[22]:
reads从左/行到右/列读取，例如分配了从Ductal到Ngn3 low EP的可信过渡。
此表可以通过叠加在 UMAP 嵌入上的定向图进行汇总展示。
[23]:
scv.pl.paga(adata, basis='umap', size=50, alpha=.1,
min_edge_width=2, node_size_scale=1.5)
文中链接
此处: https://scvelo.readthedocs.io/scvelo.datasets.pancreas.html
[2]胰腺数据: https://scvelo.readthedocs.io/scvelo.datasets.pancreas
[3]anndata docs: https://anndata.readthedocs.io/
[4]最佳实践: https://doi.org/10.15252/msb.20188746
[5][scanpy tutorial]: https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html
[6]GIF: https://user-images.githubusercontent.com/31883718/80227452-eb822480-864d-11ea-9399-56886c5e2785.gif
[7]PAGA: https://doi.org/10.1186/s13059-019-1663-x
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注