单细胞分析工具--celldancer细胞速率分析

cellDancer是美国休斯顿卫理公会研究所的助理教授Guangyu Wang团队开发的RNA速率分析新工具，于2023年4月发表于nature biotechnology。该工具基于深度学习框架预测细胞特异的速率参数(α, β and γ)，相较于之前的scvelo等RNA速率分析工具可有效预测 transcriptional boost, multi-lineage forward, multi-lineage backward等复杂情况下的单细胞数据发育预测分析。

Paper：https://www.nature.com/articles/s41587-023-01728-5

Github：https://github.com/GuangyuWangLab2021/cellDancer

Tutorial：https://guangyuwanglab2021.github.io/cellDancer_website

文章概述

（1）模型算法简介

首先根据velocity算法计算出每个细胞的每个基因的unspliced与spliced丰度；
然后将其作为DNN神经网络的输入层，模型输出是3个速率参数；
再根据预测的参数计算该基因未来下一状态的unspliced与spliced丰度；
将计算的unspliced与spliced丰度与该基因所在细胞的邻居unspliced与spliced丰度相似度作为损失函数指标；
在不断优化迭代过程中，得到每个细胞中每个基因的最佳速率参数；
基于上述得到的速率参数，进行后续的下游分析，例如细胞速率推测，伪时间分析等。

（2）复杂情况数据–transcription boost

Transcriptional boost refers to a boost in the expression induced by a change in the transcription rate.即基因在发育过程出现异常的激活
Barile鉴定了89个在红系细胞分化过程中出现transcription boost现象的基因。
这些基因会导致传统工具(scVelo)导致错误的预测，而cellDancer而进行有效的处理

（3）复杂情况数据–multi-lineage异质性

小鼠海马体有5个发育分支谱系，分别是
- dentate gyrus granule neurons
- pyramidal neurons in subiculum and CA1
- pyramidal neurons in CA2/3/4
- oligodendrocyte precursors (OPCs)
- astrocyte
由于在不同分支中，基因可能有不同的表达模式，称为branching gene
branching gene在不同分支中具有multiple dynamics；与之相对的则只有mono-kinetic
如下图所示 cellDancer不仅能够正确预测细胞发育轨迹，也能正确鉴定出branching genes

（4）衍生用法的介绍

小鼠胰腺可发育成4种细胞类型，简称为alpha-, beta-, delta-，epsilon-cells
cellDancer在进一步预测发育轨迹的基础上，可使用每个细胞的每个基因的3个速率参数作为鉴定细胞类型的方式之一
此外cellDancer的分析结果也可无缝衔接到dynamo分析框架，进行多样的后续下游分析。

此外文章进一步结合scEU-seq数据等角度验证cellDancer分析结果的可靠性、鲁棒性等，在此就不记录了，具体看阅读原文。

软件用法

作者提供了一个全面的教程资源以及示例数据可供学习；

https://guangyuwanglab2021.github.io/cellDancer_website/
在此以其中一个示例数据学习其主要的使用方法。

根据教程，搭建conda环境，并下载示例数据（Case 1）

1
2
3
4


conda create -n cellDancer python==3.7.6
conda activate cellDancer
pip install celldancer
# jupyter lab 环境下分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


import os
import sys
import glob
import pandas as pd
import math
import matplotlib.pyplot as plt
import celldancer as cd
import celldancer.cdplt as cdplt
from celldancer.cdplt import colormap
os.mkdir("./tmp")

（1）读取上游velocyto/Seurat分析结果

(1) 每个细胞的每个基因的 unsplice与splice比例
(2) 细胞的分群与UMAP降维结果

1
2
3
4
5
6
7
8
9


cell_type_u_s_path='GastrulationErythroid_cell_type_u_s.csv'
cell_type_u_s=pd.read_csv(cell_type_u_s_path)
cell_type_u_s.head()
# 	gene_name	unsplice	splice	cellID	clusters	embedding1	embedding2
# 0	Sox17	0.0	0.043971	cell_363	Blood progenitors 2	3.460521	15.574629
# 1	Sox17	0.0	0.000000	cell_382	Blood progenitors 2	2.490433	14.971734
# 2	Sox17	0.0	0.018161	cell_385	Blood progenitors 2	2.351203	15.267069
# 3	Sox17	0.0	0.000000	cell_393	Blood progenitors 2	5.899098	14.388825
# 4	Sox17	0.0	0.000000	cell_398	Blood progenitors 2	4.823139	15.374831

（2）每个细胞的每个基因速率分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32


# 以其中的4个基因为例
gene_list=['Smarca2', 'Rbms2', 'Myo1b', 'Hba-x']
loss_df, cellDancer_df=cd.velocity(cell_type_u_s,\
                                   gene_list=gene_list,\
                                   permutation_ratio=0.125,\
                                   n_jobs=-1,
                                   save_path = "./tmp")
cellDancer_df.to_csv("cellDancer_estimation.csv",index=False)

# 可视化
ncols=2
height=math.ceil(len(gene_list)/2)*4
fig = plt.figure(figsize=(10,height))

for i in range(len(gene_list)):
    ax = fig.add_subplot(math.ceil(len(gene_list)/ncols), ncols, i+1)
    cdplt.scatter_gene(
        ax=ax,
        x='splice',
        y='unsplice',
        cellDancer_df=cellDancer_df,
        custom_xlim=None,
        custom_ylim=None,
        colors=colormap.colormap_erythroid,
        alpha=0.5,
        s = 5,
        velocity=True,
        gene=gene_list[i])

    ax.set_title(gene_list[i])
    ax.axis('off')
plt.show()

（3）根据上述所有基因的结果进行细胞速率分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22


cellDancer_df=cd.compute_cell_velocity(cellDancer_df=cellDancer_df, 
	projection_neighbor_choice='gene', expression_scale='power10', 
	projection_neighbor_size=10, speed_up=(100,100))
# cellDancer_df=cd.compute_cell_velocity(cellDancer_df=cellDancer_df,
# 									   speed_up=(100,100))
# ## 使用默认参数影响不大

# 可视化
fig, ax = plt.subplots(figsize=(10,10))
cdplt.scatter_cell(ax,
                   cellDancer_df,
                   colors=colormap.colormap_erythroid,
                   alpha=0.5,
                   s=10,
                   velocity=True,
                   legend='on',
                   min_mass=15,
                   arrow_grid=(20,20),
                   custom_xlim=[-6,13],
                   custom_ylim=[2,16], )
ax.axis('off')
plt.show()

（4）伪时间分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


%%capture
# set parameters
import random
# set parameters
dt = 0.05
t_total = {dt:int(10/dt)}
n_repeats = 10
# estimate pseudotime
cellDancer_df = cd.pseudo_time(cellDancer_df=cellDancer_df,
                               grid=(30,30),
                               dt=dt,
                               t_total=t_total[dt],
                               n_repeats=n_repeats,
                               speed_up=(100,100),
                               n_paths = 3,
                               plot_long_trajs=True,
                               psrng_seeds_diffusion=[i for i in range(n_repeats)],
                               n_jobs=8)

# 可视化
fig, ax = plt.subplots(figsize=(6,6))
im=cdplt.scatter_cell(ax,cellDancer_df, colors='pseudotime', alpha=0.5, velocity=False, custom_xlim=(-5,11), custom_ylim=(4,18))
ax.axis('off')
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


## 随着伪时间的变化， spliced RNA的变化水平
gene_list=['Smarca2', 'Rbms2', 'Myo1b', 'Hba-x']
ncols=2
height=math.ceil(len(gene_list)/ncols)*4
fig = plt.figure(figsize=(10,height))

for i in range(len(gene_list)):
    ax = fig.add_subplot(math.ceil(len(gene_list)/ncols), ncols, i+1)
    cdplt.scatter_gene(
        ax=ax,
        x='pseudotime',
        y='splice',
        cellDancer_df=cellDancer_df,
        custom_xlim=None,
        custom_ylim=None,
        colors=colormap.colormap_erythroid,
        alpha=0.5,
        s = 5,
        velocity=False,
        gene=gene_list[i])

    ax.set_title(gene_list[i])
    ax.axis('off')

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27


## 展示每个细胞的的每个基因的速率参数和unspliced/spliced变化
gene_samples=['Smarca2', 'Rbms2', 'Myo1b', 'Hba-x']

for gene in gene_samples:
    fig, ax = plt.subplots(ncols=5, figsize=(12,3))
    cdplt.scatter_cell(ax[0],cellDancer_df, colors='alpha',
                 gene=gene, velocity=False)
    cdplt.scatter_cell(ax[1],cellDancer_df, colors='beta',
                 gene=gene, velocity=False)
    cdplt.scatter_cell(ax[2],cellDancer_df, colors='gamma',
                 gene=gene, velocity=False)
    cdplt.scatter_cell(ax[3],cellDancer_df, colors='splice',
                 gene=gene, velocity=False)
    cdplt.scatter_cell(ax[4],cellDancer_df, colors='unsplice',
                 gene=gene, velocity=False)
    ax[0].axis('off')
    ax[1].axis('off')
    ax[2].axis('off')
    ax[3].axis('off')
    ax[4].axis('off')
    ax[0].set_title('alpha-'+gene)
    ax[1].set_title('beta-'+gene)
    ax[2].set_title('gamma-'+gene)
    ax[3].set_title('spliced-'+gene)
    ax[4].set_title('unspliced-'+gene)
    plt.tight_layout()
    plt.show()

文章概述#

软件用法#

文章概述

软件用法