引言
本系列讲解 R语言中scRNA-seq数据分析 教程，持续更新，欢迎关注，转发！
细胞轨迹分析
伪时间分析和 RNA 速度分析为揭示了单细胞水平上细胞状态变化的丰富细节。当研究发育中的原代组织或类器官等干细胞模型这样的复杂动态系统时，希望利用这些信息来追踪细胞从初始干细胞状态，经历各种中间状态，最终分化成不同终末细胞类型的过程。此前介绍的 PAGA 方法能在一定程度上实现这一目标，但它基于细胞群水平，灵活性有限，无法自由选择使用基于相似性还是基于速度的连接信息。其他轨迹分析工具，如 Monocle、Slingshot 和 URD，能够推测细胞状态轨迹的分支情况，并将每个细胞分配到某个分支中。但这些方法通常只依赖基于相似性的伪时间，难以调整以融入 RNA 速度信息。鉴于这些局限性，研究人员急需一种能在单细胞水平上操作、灵活整合多种信息的分析方法，包括转录组相似性、伪时间、RNA 速度，以及其他元数据（如时间点和代谢标记）和分析结果（如 CytoTRACE）。
library(reticulate)
sc <- import('scanpy')
cr <- import('cellrank')
adata_DS1 <- sc$read_h5ad('DS1_scvelo.h5ad')
接下来，可以根据伪时间、转录组相似性（通过 PCA 表示）以及 scVelo 的结果来初始化不同的核。在 CellRank 中，核的作用是计算细胞之间的转换矩阵，这些矩阵后续会被用于分析。利用现有的信息，可以构建三种核：伪时间核，它可以根据提供的伪时间（比如之前重构的扩散伪时间）来计算细胞之间的转换概率；连接性核，它通过细胞之间的转录组相似性（通常在降维空间，如 PCA 或整合的潜在表示中以欧几里得距离来估计）来计算；以及速度核，它依据之前进行的 RNA 速度分析来计算转换概率。
pk <- cr$kernels$PseudotimeKernel(adata_DS1, time_key="dpt3")$compute_transition_matrix()
ck <- cr$kernels$ConnectivityKernel(adata_DS1)$compute_transition_matrix()
vk <- cr$kernels$VelocityKernel(adata_DS1)$compute_transition_matrix()
还可以根据这些核，并结合不同的自定义权重来创建组合核。不过，从这一步开始，在 R 中运行代码会遇到一些问题。好在 reticulate 提供了一个功能（repl_python），可以直接在 R 中启动一个 Python 交互式环境。在这个环境中，依然可以访问在 R 环境中创建的对象（通过 r.object）。而且，在 Python 环境中创建的任何对象，在返回到 R 之后也可以继续访问（通过 py$object）。那么，先开启 Python 交互式会话吧。
repl_python()
现在，可以在打开的Python会话中进行分析。
import cellrank as cr
import numpy as np
combined_kernel = 0.5 * r.vk + 0.3 * r.pk + 0.2 * r.ck
已经成功创建了组合核。这个组合核是由三个核组合而成的，其中速度核在最终核中的占比为 50%，而伪时间核和连接性核的占比分别为 30% 和 20%。接下来，可以利用这个组合核来尝试推断终末状态。
g = cr.estimators.GPCCA(combined_kernel)
g.fit(n_states=15, cluster_key="celltype")
g.predict_terminal_states(method="top_n", n_states=6)
g.plot_macrostates(which="terminal")
不同类型的神经元被成功识别为终末状态。其中，有一群神经祖细胞（NPCs）被误判为潜在的终末状态。当然，从科学探索的角度讲，应该深入研究这个细胞群体，看看它是否真的代表了系统中尚未被发现的终末细胞状态，这可能会带来令人兴奋的新发现。但如果认为这只是一个误差，那么也可以手动将终末状态设定为不同已注释神经元细胞类型的随机组合子集（在下面的例子中为30个细胞）。
g = cr.estimators.GPCCA(combined_kernel)
g.set_terminal_states({"Midbrain-hindbrain boundary neuron": r.adata_DS1[r.adata_DS1.obs['celltype'] == "Midbrain-hindbrain boundary neuron"].obs_names[:30],
"Ventral telen. neuron": r.adata_DS1[r.adata_DS1.obs['celltype'] == "Ventral telen. neuron"].obs_names[:30],
"Dorsal telen. neuron": r.adata_DS1[r.adata_DS1.obs['celltype'] == "Dorsal telen. neuron"].obs_names[:30],
"Dien. and midbrain excitatory neuron": r.adata_DS1[r.adata_DS1.obs['celltype'] == "Dien. and midbrain excitatory neuron"].obs_names[:30],
"Dien. and midbrain inhibitory neuron": r.adata_DS1[r.adata_DS1.obs['celltype'] == "Dien. and midbrain inhibitory neuron"].obs_names[:30]})
g.plot_macrostates(which="terminal")
现在，可以开始计算每个细胞最终转化为每个终末状态的概率。
g.compute_fate_probabilities()
g.plot_fate_probabilities(legend_loc="right", basis='umap', same_plot=False)
然后，可以提取命运概率估计，将其保存到数据框架中
import pandas as pd
prob = pd.DataFrame(g.fate_probabilities).set_index(g.terminal_states.index).set_axis(g.terminal_states.cat.categories, axis=1)
现在，可以通过输入 exit 命令结束 Python 的交互式会话，返回到 R 环境中，并将细胞命运概率的数据框存储到 Seurat 对象里。
exit
library(Seurat)
seurat_DS1 <- readRDS(file='DS1/seurat_obj_all.rds')
seurat_DS1@meta.data[,colnames(py$prob)] <- py$prob[colnames(seurat_DS1),]
FeaturePlot(seurat_DS1, colnames(py$prob), ncol=5)
通过在细胞层面估算的命运概率，可以进一步开展分析，比如找出可能影响细胞命运变化的驱动基因。可以根据这些概率，将细胞定位到分化路径的不同分支上。例如，在先前关于 脑类器官 和 视网膜类器官 的细胞命运研究中，将单个细胞的命运概率归纳成元簇，直观展示了系统中各类细胞是如何逐步分化成特定类型的。
