1-cellrank2简介

Cellrank是啥?

cellrank是一款基于多视角单细胞转录组数据(包括RNA速率、发育潜能、拟时序、系谱追踪数据、代谢标签等)用于研究例如发育、再生、癌症、重编程等动态生物过程的工具包。从下面这张图来看,cellrank可以说是对RNA速率以及拟时等的拓展。通常我们对于单细胞转录组数据,分析拟时或者RNA速率构建轨迹,模拟生物体细胞状态动态和命运决策,cellrank2相比cellrank1,功能更加全面。


Cellrank2能干啥?
以下内容翻译自cellrank2原文摘要】:单细胞RNA测序使我们能够利用表达相似性或RNA velocity来模拟细胞状态变化和命运决策,从而重建状态变化(state-change)轨迹;然而,轨迹推断并未纳入有价值的时间点信息或利用其他数据模态,而那些针对不同数据视角的方法则往往难以整合,或在大规模数据下难以扩展。在此,我们提出 CellRank 2,一个灵活且可扩展的框架,可在统一的分析体系下利用多视角单细胞数据,对多达数百万个细胞的命运进行研究。CellRank 2 能够在人体造血和内胚层发育过程中,跨不同数据模态稳定地解析末端细胞状态及其命运概率。该框架还支持在实验时间点内及跨时间点整合细胞转变,我们利用这一特性鉴定出在咽内胚层发育过程中促进髓质胸腺上皮细胞形成的关键基因。此外,我们的方法能够基于代谢标记数据推断细胞特异性的转录与降解速率,并将其应用于肠类器官系统,以刻画分化轨迹并解析关键的调控策略。

Cellrank2有哪些应用

基于多种生物学先验信息来估计细胞分化方向,包括pseudotime, developmental potential, RNA velocity、代谢标记(metabolic labels)等。
计算初始、终末以及中间宏状态(macrostates)。
推断细胞命运概率以及驱动基因(driver genes)。
可视化并聚类基因表达趋势。 ……


官网及原文链接:cellrank2官网有超级详细的教程,我们这里使用的是一篇Nature文章的数据,原文也提供了代码,用来演示整个流程,让其更加清晰化!

■ cellrank2 github:https://github.com/theislab/cellrank
■ cellrank2 paper:https://www./articles/s41592-024-02303-9
■ paper info: Weiler, P., Lange, M., Klein, M. et al. CellRank 2: unified fate mapping in multiview single-cell data. Nat Methods 21, 1196–1205 (2024). https:///10.1038/s41592-024-02303-9


Nature文章及数据链接
Data:https:///records/15857303
Paper:https://www./articles/s41586-025-09503-z


2-cellrank2分析流程

终端创建一个单独的环境,安装cellrankconda create -n cellranks python==3.10conda activate cellrankspip install cellrank -i https://pypi.tuna.tsinghua.edu.cn/simple
import cellrank as crimport scanpy as scimport numpy as npimport pandas as pdfrom matplotlib import rcParamsimport matplotlib as mplimport matplotlib.pyplot as pltimport matplotlib.cm as cm # color mapsfrom matplotlib.pyplot import rc_contextfrom matplotlib.colors import LinearSegmentedColormapimport seaborn as snscr.settings.verbosity = 3
#这个数据使用的是DPT分析的拟时,您自己的数据无论是哪种轨迹构建方法,都可以adata = sc.read_h5ad("../cellrank2/datasets/FAprojection_DPT_final.h5ad")
pheno_col = [    "#66A61E",  # already hex    "#40E0D0",  # turquoise    "red",      # brown2    "#68228B",  # darkorchid4    "#1E90FF",  # dodgerblue    "#00868B",  # turquoise4    "#00868B",  # turquoise4    "#00868B",  # turquoise4    "#FFA500",  # orange]
cmap = plt.cm.get_cmap("Spectral")sc.pl.draw_graph(adata, color=['dpt_pseudotime','cell_type'], legend_loc='right margin', color_map=cmap, palette=pheno_col)

选择kernel,计算转移矩阵:因为演示数据是pseudotime分析,所以选择PseudotimeKernel分析,具体的自己的数据按照自己实际的情况选择,只是构建转移矩阵这里不同,后续的流程都是一样。
from cellrank.kernels import PseudotimeKernelpk = PseudotimeKernel(adata, time_key="dpt_pseudotime"#time_key,adata.obs储存拟时的列
pk.compute_transition_matrix()
#非必要,小提琴图看看celltype pseudotime分布fig, ax = plt.subplots(figsize=(105))sc.pl.violin(        adata,        keys=["dpt_pseudotime"],        groupby="cell_type",        rotation=45,        ax=ax,    )
#通过流线可视化转移矩阵,注意,这不是RNA速率,是指转移矩阵在二维空间的嵌入fig, ax = plt.subplots(figsize=(77))pk.plot_projection(            color="cell_type",            recompute=True,            basis="X_draw_graph_fa",            title="",            legend_loc="on data",            alpha=0.25,            linewidth=2,            ax=ax,        )

Estimator-Computing Initial and Terminal States:这里我们使用GPCCA
g = cr.estimators.GPCCA(pk)
#转移矩阵Schur分解,提取最相关的特征值与特征向量【Schur分解是矩阵分解方法,具体原理感兴趣的可查阅线性代数】g.compute_schur(n_components=20)g.plot_spectrum(real_only=True)

#原文先在后面识别宏状态的时候,cluster对应的是leiden_scVI_1.2细分群。sc.pl.draw_graph(adata, color=['leiden_scVI_1.2'], legend_loc='right margin')
#划分细胞状态,g.compute_macrostates参数中的n_states就是上一步骤中g.plot_spectrum确定选择的值。#但是注意,实际设置中为了得到更多的细胞状态,可以自行设置值的大小,例如原文献中就是自行设定【依据自己实际需求核对结果的解释】。#这个值不宜设置过大,可能会引入噪声。这里设置9,能够覆盖到这个数据中所有的9种celltype。g.compute_macrostates(n_states=9, cluster_key="leiden_scVI_1.2")#把计算好的宏观状态可视化在低维嵌入降维图上。每个宏观状态会用一种颜色标注,让你直观看到哪些细胞属于哪个终点/中间态。g.plot_macrostates(which="all",basis='X_draw_graph_fa',legend_loc="right", s=100)

g.predict_terminal_states(allow_overlap=True)g.plot_macrostates(which="terminal",basis='X_draw_graph_fa')
#识别起始状态g.predict_initial_states(allow_overlap=True)g.plot_macrostates(which="initial",basis='X_draw_graph_fa')

Estimating Fate Probabilities:计算每个细胞“最终”落入每个terminal state的概率。

Cellrank2演示教程:在多模态单细胞数据中的统一细胞命运推断
#命运概率计算及其可视化g.compute_fate_probabilities()g.plot_fate_probabilities(legend_loc="right",basis="X_draw_graph_fa")
#单个可视化命运概率g.plot_fate_probabilities(same_plot=False, basis="X_draw_graph_fa",figsize=(55),ncols=3)

#环形多边形图展示多个向terminal 发展的命运概率palette=[  "#E41A1C"# strong red  "#377EB8"# medium blue  "#4DAF4A"# green  "#984EA3"# purple  "#FF7F00"# orange  "#FFFF33"# yellow  "#A65628"# brown  "#e7298a"# pink  "#666666"# grey  "#66C2A5"# teal  "#FC8D62"# salmon  "#8DA0CB"# soft blue  "#E78AC3"# soft pink (different from 8)  "#A6D854"# light green (but yellowish tint, not green)  "#FFD92F"# lemon yellow  "#E5C494"# light brown  "#B3B3B3"# light grey  "#1B9E77"# deep teal  "#D95F02"# dark orange  "#7570B3"# strong purple  "#66A61E"]  # olive green (NOT same green as before)cr.pl.circular_projection(adata, keys="leiden_scVI_1.2", legend_loc="right", palette=palette)

driver genes:最后,使用“compute_lineage_drivers()”将命运概率与基因表达进行关联,从而发现潜在的驱动基因,分析得到的结果是一个基因 × lineage的相关性表。

#这里计算所有的lineage drivers genesall_drivers = g.compute_lineage_drivers()

3-可视化:


采用广义加性模型(GAM),根据每个细胞在每个轨迹中的贡献程度(该贡献度由其命运概率向量决定)对其进行加权。首先需要初始化一个模型,注意需要安装rpy2version>=3.3.0`。

model = cr.models.GAMR(adata, n_knots=6, smoothing_penalty=10.0)

可视化拟时基因表达趋势 这里我们随便挑选几个基因作为演示。根据命运概率和伪时间,可以绘制出特定轨迹下的基因表达趋势图。

cr.pl.gene_trends(    adata,    model=model,    genes=['Ascl1','Neurod1','Atoh1','Yap1'],    same_plot=True,    ncols=2,    time_key="dpt_pseudotime",    hide_cells=True,    weight_threshold=(1e-31e-3),)

热图可视化某个lineage轨迹基因
# Compute lineage driversNE_drivers = g.compute_lineage_drivers(lineages="0",model=model)genes = NE_drivers.head(100).index
cr.pl.heatmap(    adata,    model=model,    lineages=["0"],    return_models=False,    genes=genes,     cbar=False,    time_key="dpt_pseudotime",    show_fate_probabilities=False,    show_all_genes=True, cluster_key="leiden_scVI_1.2", cluster_genes=False, figsize=(816))

其他感兴趣的可视化可以详细阅读网关可视化教程!
觉得分享有用的点个赞再走呗!