scRNA-seqScience相关 杂单细胞测序

10X单细胞(10X空间转录组)基于马尔可夫过程和RNA速率的细

2021-10-20  本文已影响0人  单细胞空间交响乐

这一篇接上一篇,我们来分享Cellrank的基础分析代码,我们直接开始

使用 RNA 速度和转录组学相似性来估计细胞-细胞转换概率。即使没有 RNA 速度信息,也可以应用 CellRank(但还是推荐大家RNA 速率和相似度方法联合使用)。

加载

import sys

if "google.colab" in sys.modules:
    !pip install -q git+https://github.com/theislab/cellrank@dev
    !pip install python-igraph
import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np

scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo")
cr.settings.verbosity = 2
import warnings

warnings.simplefilter("ignore", category=UserWarning)
warnings.simplefilter("ignore", category=FutureWarning)
warnings.simplefilter("ignore", category=DeprecationWarning)

首先,需要获取数据。 以下命令将下载 adata 对象并将其保存在 datasets/endocrinogenesis_day15.5.h5ad 下(示例数据,大家可以下载下来自己研究)。 我们还将显示拼接/未拼接读数的比例,需要用它来估计 RNA 速度。

adata = cr.datasets.pancreas()
scv.pl.proportions(adata)
adata
图片.png

AnnData object with n_obs × n_vars = 2531 × 27998
obs: 'day', 'proliferation', 'G2M_score', 'S_score', 'phase', 'clusters_coarse', 'clusters', 'clusters_fine', 'louvain_Alpha', 'louvain_Beta', 'palantir_pseudotime'
var: 'highly_variable_genes'
uns: 'clusters_colors', 'clusters_fine_colors', 'day_colors', 'louvain_Alpha_colors', 'louvain_Beta_colors', 'neighbors', 'pca'
obsm: 'X_pca', 'X_umap'
layers: 'spliced', 'unspliced'
obsp: 'connectivities', 'distances'

前处理数据

过滤掉没有足够剪接/未剪接计数的基因,对数据进行归一化和对数变换,并限制在高度可变的基因上。 此外,计算速度估计的主成分和矩。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=30)
scv.pp.moments(adata, n_pcs=None, n_neighbors=None)

Run scVelo

我们将使用来自 scVelo 的动力学模型来估计速度。
scv.tl.recover_dynamics(adata, n_jobs=8)
一旦有了参数,就可以使用这些参数来计算速度和速度图。 速度图是一个加权图,它指定了两个cell在给定速度向量和相对位置的情况下转换为另一个cell的可能性。
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(
    adata, basis="umap", legend_fontsize=12, title="", smooth=0.8, min_mass=4
)
图片.png

运行Cellrank

CellRank 提供了多种将方向性注入单细胞数据的方法。 在这里,方向信息来自 RNA 速度,使用这些信息来计算胰腺发育动态过程的初始和终止状态以及fate probabilities。

Identify terminal states

cr.tl.terminal_states(adata, cluster_key="clusters", weight_connectivities=0.2)

The most important parameters in the above function are:

cr.pl.terminal_states(adata)
图片.png

Identify initial states

cr.tl.initial_states(adata, cluster_key="clusters")
cr.pl.initial_states(adata, discrete=True)
图片.png

Compute fate maps

一旦知道终端状态,就可以计算相关的命运图——对于每个细胞,寻求细胞朝着每个确定的终端状态发展的可能性有多大。
cr.tl.lineages(adata)
cr.pl.lineages(adata, same_plot=False)
图片.png
可以将上述内容聚合成一个单一的全局命运图,其中将每个终端状态与颜色相关联,并使用该颜色的强度来显示每个单个细胞的命运:
cr.pl.lineages(adata, same_plot=True)
图片.png

Directed PAGA

我们可以使用具有有向边的 [Wolf et al., 2019] 的改编版本将个体命运图进一步聚合成集群级别的命运图。 我们首先使用 CellRank 识别的 root_key 和 end_key 计算 scVelo 的潜伏时间,它们分别是初始状态或终止状态的概率。
scv.tl.recover_latent_time(
    adata, root_key="initial_states_probs", end_key="terminal_states_probs"
)
Next, we can use the inferred pseudotime along with the initial and terminal states probabilities to compute the directed PAGA.
scv.tl.paga(
    adata,
    groups="clusters",
    root_key="initial_states_probs",
    end_key="terminal_states_probs",
    use_time_prior="velocity_pseudotime",
)
cr.pl.cluster_fates(
    adata,
    mode="paga_pie",
    cluster_key="clusters",
    basis="umap",
    legend_kwargs={"loc": "top right out"},
    legend_loc="top left out",
    node_size_scale=5,
    edge_width_scale=1,
    max_edge_width=4,
    title="directed PAGA",
)
图片.png

Compute lineage drivers

可以计算所有谱系或部分谱系的驱动基因。 还可以通过指定cluster=...将其限制为某些cluster。 在生成的数据框中,还看到了 p 值、校正后的 p 值(q 值)和相关统计量的 95% 置信区间。
图片.png
Afterwards, we can plot the top 5 driver genes (based on the correlation), e.g. for the Alpha lineage
cr.pl.lineage_drivers(adata, lineage="Alpha", n_genes=5)
图片.png

Gene expression trends

上面演示的功能是 CellRank 的主要功能:计算初始和终止状态以及概率命运图。 现在可以使用计算出的概率来例如 沿谱系平滑基因表达趋势。

从细胞的时间顺序开始。 为了得到这个,可以计算 scVelo 的潜伏时间,如前所述,或者,我们可以只使用 CellRank 的初始状态来计算 。

# compue DPT, starting from CellRank defined root cell
root_idx = np.where(adata.obs["initial_states"] == "Ngn3 low EP")[0][0]
adata.uns["iroot"] = root_idx
sc.tl.dpt(adata)

scv.pl.scatter(
    adata,
    color=["clusters", root_idx, "latent_time", "dpt_pseudotime"],
    fontsize=16,
    cmap="viridis",
    perc=[2, 98],
    colorbar=True,
    rescale_color=[0, 1],
    title=["clusters", "root cell", "latent time", "dpt pseudotime"],
)
图片.png

We can plot dynamics of genes in pseudotime along individual trajectories, defined via the fate maps we computed above.

model = cr.ul.models.GAM(adata)
cr.pl.gene_trends(
    adata,
    model=model,
    data_key="X",
    genes=["Pak3", "Neurog3", "Ghrl"],
    ncols=3,
    time_key="latent_time",
    same_plot=True,
    hide_cells=True,
    figsize=(15, 4),
    n_test_points=200,
)
图片.png
还可以在热图中可视化上面计算的谱系驱动程序。 下面,对 Alpha 谱系执行此操作,即在伪时间中平滑假定的 Alpha 驱动程序的基因表达,将 Alpha 命运概率用作细胞级别的权重。 根据它们在伪时间中的峰值对基因进行排序,从而揭示了基因表达事件的级联。
cr.pl.heatmap(
    adata,
    model,
    genes=adata.varm['terminal_lineage_drivers']["Alpha_corr"].sort_values(ascending=False).index[:100],
    show_absorption_probabilities=True,
    lineages="Alpha",
    n_jobs=1,
    backend="loky",
)
图片.png

生活很好,有你更好

上一篇 下一篇

猜你喜欢

热点阅读