ixxmu
commented
8 months ago Nature 复现 | 单细胞Perturbation探寻不同物种的细胞命运(附代码) by 苏州新格元生物科技有限公司
图谱正当时，以人类细胞图谱计划(Human Cell Atlas)为引领的单细胞数据是CNS顶刊的常客。以模式物种为研究材料构建了成熟的单细胞实验与标准分析流程，其他物种的单细胞图谱正也在如火如荼地进行着，一个值得思考的基本问题是：构建完物种的单细胞图谱之后，我们下一步该如何利用这些珍贵的数据资源？
2023年8月17日，JS Fleck等人在Science杂志发表文章《What is a cell type? 》，文章认为一个细胞类型或状态不仅涉及到细胞的发育历史和其当前的分子特征，还与细胞对环境变化或其他perturbations 的响应方式相关。现有的细胞图谱描述了细胞的当前特征，有时也包括其发育历史，但无法提供决定发育和预测未来细胞状态的因素信息。若能从细胞的当前状态预测其响应，则细胞图谱将成为进行计算模拟筛选的丰富资源。因此针对上述的问题，JS Fleck等人建议：A next step for cell atlases should be to chart perturbations in human model systems。
那么，perturbation 是什么？有哪些模型可以来用？今天，就让我们一起在复现CellOracle分析的过程中，体验一把人工智能模型的魅力！目前CellOracle支持小鼠、果蝇、拟南芥、斑马鱼等10多种物种的分析。
扰动(perturbation)是指任何改变细胞表型的输入。它们可以是自然发生的，也可以是实验性的，包括基因变化、化学、物理或电刺激以及病原体感染。扰动引起的变化的复杂性可以有很大差异：某些遗传性功能丧失突变只影响单一信号通路，而其他突变可能会改变多个细胞过程。正如文章新格元推出果蝇预训练大模型：探索动物单细胞研究的新途径指出Perturbation预测已经成为单细胞研究领域的热点之一：
2023年Nature正刊分别发表两篇关于基于单细胞数据进行Perturbation的方法学文章，相关软件为Geneformer和CellOracle。
2023年7月，结合GO知识图谱做单细胞perturbation的方法学文章发表在Nature Biotechnology。
2023年8月30日，机器学习顶会NeurIPS 2023的单细胞挑战赛的题目正式定为：Single-Cell Perturbation Prediction。
新格元数据科学团队推出基于单细胞数据的perturbation系列文章，今天我们先来分享下CellOracle的分析方法。
此次的分析主要是参考2023年2月8日美国圣路易斯华盛顿大学Samantha A. Morris课题组发表在Nature正刊的文章Dissecting cell identity via network inference and in silico gene perturbation，作者主要开发了一种使用从单细胞多组学数据(scRNA-seq和scATAC-seq)推断出的基因调控网络来执行计算机转录因子扰动的方法:CellOracle。CellOracle只需要利用野生型单细胞的多组学数据scRNA-seq和scATAC-seq（非必要），即可准确推断出敲除特定转录因子后而导致人源、小鼠和斑马鱼的相关细胞的表型变化。
@Dissecting cell identity via network inference and in silico gene perturbation
本次复现的主要内容(Fig. 1), 基于小鼠原始的骨髓细胞分化图谱，用CellOracle模拟敲除Spi1和Gata1，得到两者截然不同的“敲除分数”，这与之前发表的关于Gata1和Spi1控制的髓系细胞分化转变实验结果吻合。需要指出的是，尽管Gata1在粒细胞的表达较低，CellOracle也能预测出其敲除后早期粒细胞分化中细胞表型的变化，这是常规的基因差异表达分析无法捕获的。
Fig. 1 小鼠骨髓细胞分化图谱及CellOracle敲除推断
CellOracle支持多种数据格式来进行基础GRN构建，如Fig.2所示，可以是scATAC-seq数据，bulk ATAC-seq数据，物种的promoter database，或者用户提供的自定义转录因子和目标基因配对信息。此处省去数据类型处理的过程，例如scATAC-seq数据处理得到下面的信息，第一列为序列号，第二列为peak id，第三列为基因名字。
Fig. 2 CellOracle构建基于单细胞多组学数据的细胞基因调控网络
Table 1 | 基于scATAC-seq数据的输入格式
# 构建 TFinfo 对象
tfi = ma.TFinfo(peak_data_frame=peaks, ref_genome=ref_genome, genomes_dir=None)
# 利用motifs和peaks寻找转录因子结合位点
tfi.scan(fpr=0.02, motifs=motifs, verbose=True)
# scores值过滤
tfi.reset_filtering()
tfi.filter_motifs_by_score(threshold=10)
tfi.make_TFinfo_dataframe_and_dictionary(verbose=True)
# 按dataframe格式保存结果
df = tfi.to_dataframe()
df.to_parquet("base_GRN_dataframe.parquet")运行上面的代码即可利用peaks信息和motifs信息构建基础的GRN结果，见下表，分别代表了每个peak与转录因子的互作强度。
Table 2 | 基础GRN结果展示
得到基础GRN后，我们加入单细胞数据获得细胞特异的GRN，此处单细胞数据要求输入为raw_counts和细胞需要有机器分群的标签或者细胞注释的标签，如下面数据中基于机器分群的细胞注释结果“louvain_annot”，此外还需有细胞原始（无perturbation）的细胞发育轨迹推断的指标，如下面数据中的“dpt_pseudotime”。
# louvain_annot后续用于作为CellOracle分析的细胞类群
sc.pl.draw_graph(oracle.adata, color="louvain_annot")
# dpt_pseudotime后续用于CellOracle做细胞命运推断
sc.pl.draw_graph(oracle.adata, color="dpt_pseudotime")
Fig. 3 后续用于CellOracle细胞命运推断的细胞类群标签和细胞发育轨迹推断的指标
准备好上述的单细胞数据之后，可以把单细胞表达数据和上述构建的基础GRN放到一个变量“oracle”。值得注意的是CellOracle在进行细胞特异的GRN推断的时候，会先使用与velocyto相同的策略对细胞数据进行imputation，然后再此基础上计算细胞内网络结构。
# 用单细胞表达数据adata初始化oracle
oracle.import_anndata_as_raw_count(adata=adata, cluster_column_name="louvain_annot", embedding_name="X_draw_graph_fa")
# 把之前构建好的GRN加入到oracle
oracle.import_TF_data(TF_info_matrix=base_GRN)
# PCA降维和knn imputation（数据补插）
oracle.perform_PCA()
oracle.knn_imputation(n_pca_dims=n_comps, k=k, balanced=True, b_sight=k*8, b_maxl=k*4, n_jobs=4)
# 根据细胞标签louvain_annot"来进行细胞特异GRN推断
# 此步骤消耗的时间较长
links = oracle.get_links(cluster_name_for_GRN_unit="louvain_annot", alpha=10, verbose_level=10)
执行完上述的代码即可获得细胞类群特异的GRN网络，例如下面以细胞类群“Ery_0”为展示，每一行代表了转录因子（source）和其调控目标基因（target）的连接关系：
Table 3 | CellOracle推断的细胞类群特异的GRN
选择文章的基因Gata1和Spi1进行分析，CellOracle准备好细胞类群特定的GRN后，进行基因敲除perturbation的分析比较简单方便，只需要执行下面的代码，下面以Gata1为例子：
# 设置需要敲除的基因名字
goi = "Gata1"
# 将其表达值设置为0，然后基于GRN来做propagation推断
oracle.simulate_shift( perturb_condition={goi: 0.0}, n_propagation=3 )
# Get transition probability
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
# Calculate embedding
oracle.calculate_embedding_shift(sigma_corr=0.05)
# 可视化结果，scale控制箭头大小
scale = 25
oracle.plot_quiver(scale=scale)
plt.show()
通过执行上面的代码，我们即可获得敲除Gata1后细胞命运分化的推断，如下图的中间umap图，通过加入原始细胞命运分化的轨迹图（需要人为设置root_cell群），如下面的左umap图，按照下图所示的向量方向叠加计算Perturbation Score，就可以获得敲除Gata1对原先细胞命运分化的影响，例如对GMP起到促进的作用，而对MEP起到抑制的作用，使用同样的代码即可获得其他基因敲除对细胞命运分化的影响，例如Spi1。
Fig. 4 CellOracle基于GRN推断敲除Gata1对原先细胞命运分化的影响
值得注意的是，与常规的差异基因表达分析不一样，对细胞命运分化的影响并不是主要决定于该基因的表达情况，而是通过整体的GRN来推断的，譬如Gata1基因在GMP细胞类群的表达属于中低水平（见图Fig. 5a），如果单单从基因表达来是无法正确推断该基因对这个细胞类群的影响（Gata2被认为在早期MEP和GMP群体中发挥着重要作用），而CellOracle能从整个GRN推断出Gata2在GMP细胞类群的基因网络起着重要的调控作用，例如图中Fig. 5c所示Gata2在GMP细胞类群获得高的网络分数值，执行下面的代码即可获得基因在网络的各种分数值(Fig .5c)。
# Visualize Gata2 network score dynamics
links.plot_score_per_cluster(goi="Gata2", save=f"{save_folder}/network_score_per_gene/")
Fig. 5 CellOracle相对常规的表达分析能更加准确推断基因对基因网络的重要性
复现CellOracle关于基因敲除Perturbation的分析之后，基于此分析的结果，我们可以获得不同细胞类型敲除不同基因后对细胞表型的影响，这是常规基因表达差异分析所欠缺的。如果您也想在单细胞水平研究Perturbation，可与新格元当地销售联系或扫描下方二维码。
- THE END -
往期推荐
新格元推出果蝇预训练大模型：探索动物单细胞研究的新途径
啥？玉米单细胞研究还可以这样做
訇然中开 | 海洋生物单细胞测序的新格元方案
供稿：市场支持中心
审核：市场支持中心
想了解更多关于单细胞测序信息，欢迎点击“阅读原文”留下联系方式，我们将安排同事与您对接。
https://mp.weixin.qq.com/s/-E9zcB7nM6DIjz6GIzUIwg