ixxmu
commented
9 months ago SIMBA（5）-单细胞多特征统一embedding-基因调控网络推断 by TOP生物信息
【基于python的单细胞分析】回顾：
【基于python的单细胞分析】使用scVI实现批次效应校正
如何进行细胞类型注释
【SIMBA系列教程】回顾：
系列教程：SIMBA（1）-单细胞多特征统一embedding-算法介绍
系列教程：SIMBA（2）-单细胞多特征统一embedding-scRNA-seq分析
系列教程：SIMBA（3）-单细胞多特征统一embedding-scATAC-seq分析
系列教程：SIMBA（4）-单细胞多特征统一embedding-双组学分析
文章简介
《SIMBA: single-cell embedding along with features》于今年5月29日发表在Nature Methods上。我们将用一期推送，介绍SIMBA的文章，然后使用6期推送，系统介绍SIMBA算法的相关应用，以及代码复现，包括scRNA-seq分析、scATAC-seq分析、scRNA+scATAC-seq分析、scRNA+scATAC-seq推断GRN、scRNA-seq批次效应校正和scRNA+scATAC-seq整合分析。
过去我们如何推断基因调控网络（GRN）
最早期，也是最简单的推断GRN的算法是基于加权基因共表达网络算法来推断GRN，这种算法它整个转录组范围内进行成对的相关性分析，以识别共表达基因模块。尽管这种策略有助于无监督地识别基因模块，但缺乏因果调控连接限制了其解释性，并且通常会产生大量的假阳性关联。
再后来，GENIE3以及GRN-Boost方法首先以先验知识的方法从基因中识别出转录因子（TF），然后训练模型，仅基于转录因子的表达来预测目标基因的表达，从而显著减少了要考虑的相互作用数量。同时通过这样做，无向相互作用被转化为有向连接，从而引入了可能的因果关系。
虽然基于bulk RNA-seq的方法可以有效推断出显著的GRN，但在混合样本（如组织）的情况下，它们无法捕捉到GRNs的细胞类型或状态特异性。单细胞测序使得细胞类型特异的TF-基因相互作用，以及这些GRNs在发育和条件下发生的动态变化推断成为可能。
然而，仅从转录组学数据进行推断会引入假阳性，因为许多参与基因调控的其他机制，如染色质可及性被忽略了。单细胞染色质可及性分析（如scATAC-seq）可与单细胞转录组学一起进行测序，双组学的结果可以更加有效地推断GRN。
具体来说，它们从TF基因表达预测基因表达，利用结合Peaks信息将TF分配到开放的顺势调控元件（CREs），并将CREs与受基因组距离限制的目标基因联系起来。
以上就是目前我们推断GRN的方法，不同方法之间大同小异，但是基本都是这套方法。
SIMBA如何推断GRN
确定关键转录因子（TF）
为了先验地定义主调控因子，我们假定其 TF 基序和 TF 基因都应该是细胞类型特异性的，由于活性基因调控涉及 TF 的表达和其结合位点的可及性。因此，TF 基序和 TF 基因应该紧密嵌入在共享的embedding空间中。
为了确定不同细胞类型的特异性的关键转录因子，对于每一个TF的motif，我们首先在embedding空间中计算它与其他所有基因的距离，然后我们可以得到这个TF的motif在所有基因（包括非TF基因）中的距离的排名，通过这个排名，我们可以衡量其显著性。排名越低，则TF和TF motif在embedding空间上的距离更加接近，则更加显著。
此外，我们还可以使用SIMBA的相关指标鉴定细胞类型特异性的关键转录因子。
确定TF调控的靶基因
SIMBA根据以下假设确定TF调控的靶基因：（1）在共享的 SIMBA embedding空间中，靶基因接近 TF motif和 TF 基因，表明靶基因的表达与 TF 的表达和 TF motif的可及性以细胞类型特异性的方式高度相关。（2）靶基因位点附近的可及区域(峰)必须接近 TF 基序和靶 TF 基因，表明靶基因位点附近的顺式调节元件的可及性与 TF 的表达和 TF 基序的可及性以细胞类型特异性的方式高度相关。
对于一个给定的关键TF，SIMBA分别计算TF motif和TF周围k个（k=200）最近邻基因。这些相邻基因的交集是最初候选靶基因的集合，然后根据假定目标基因的 TSS 上游或下游100kb 内的开放区域(峰)必须含有 TF 基序的标准筛选这些基因。
然后，对于每个候选靶基因，在 SIMBA 嵌入空间中计算4种类型的距离: (1)候选靶基因与 TF 基因的嵌入距离; (2)候选靶基因与 TF 基序的嵌入距离; (3)候选靶基因与 TF 基序基因座附近的峰值距离; (4)候选靶基因与候选基因基因座附近的峰值距离。所有的距离(默认的欧几里德距离)被转换为所有基因或所有峰值的等级，以使不同关键TF之间的计算得到距离具有可比性。
目标基因的最终列表使用计算出的等级来决定，使用两个标准: (1)至少有一个与 TF 基因或 TF 基序最接近的峰在预定范围内(默认为前1,000个) ; 和(2)候选目标基因的平均等级在预定范围内(默认为前5,000个)。
代码实战
首先我们下载SIMBA包及其依赖，导入这些包，并进行绘图参数设置
!pip install git+https://github.com/pinellolab/simba_pbg.git
!pip install git+https://github.com/huidongchen/simba
import os
import simba as si
workdir = 'result_simba_atacseq'
si.settings.set_workdir(workdir)
si.settings.set_figure_params(dpi=80,
style='white',
fig_size=[5,5],
rc={'image.cmap': 'viridis'})
# make plots prettier
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina')
我们首先导入数据，adata_G代表生成的Gene的embedding，adata_M代表生成的Motif的embedding，而adata_cmp_CG和adata_cmp_CM则代表经过计算过相应指标的adata_G和adata_M，这些数据都是在第四次推送-多组学分析中计算得到的结果。
adata_G = si.read_h5ad(os.path.join(workdir,'adata_G.h5ad'))
adata_M = si.read_h5ad(os.path.join(workdir,'adata_M.h5ad'))
adata_all = si.read_h5ad(os.path.join(workdir,'adata_all.h5ad'))
adata_cmp_CG = si.read_h5ad(os.path.join(workdir,'adata_cmp_CG.h5ad'))
adata_cmp_CM = si.read_h5ad(os.path.join(workdir,'adata_cmp_CM.h5ad'))
首先我们从这些数据中生成配对的Gene和Motif，这些配对的数据将被用于后续的下游分析。
# find paired TF motifs and TF genes
motifs_genes = pd.DataFrame(columns=['motif', 'gene'])
for x in adata_M.obs_names:
x_split = x.split('_')
for y in adata_G.obs_names:
if y in x_split:
motifs_genes.loc[motifs_genes.shape[0]] = [x,y]
print(motifs_genes.shape)
motifs_genes.head()
我们提取出Gene和Motif的ID，以及这些feature的参数用于后续的分析。
list_tf_motif = motifs_genes['motif'].tolist()
list_tf_gene = motifs_genes['gene'].tolist()
df_metrics_motif = adata_cmp_CM.var.copy()
df_metrics_gene = adata_cmp_CG.var.copy()
有了这些参数，我们就可以从adata_all文件里鉴定关键转录因子了，这里的adata_all就是双组学分析中所有特征的联合embedding。
df_MR是一个dataframe，包含鉴定关键转录因子的函数计算出来的结果，即满足显著特异性的转录因子，此外也包含这些转录因子相应的指标信息，用于后续的分析。
df_MR = si.tl.find_master_regulators(adata_all,
list_tf_motif=list_tf_motif,
list_tf_gene=list_tf_gene,
cutoff_gene_max=1.5,
cutoff_gene_gini=0.35,
cutoff_motif_max=3,
cutoff_motif_gini=0.7,
metrics_gene=df_metrics_gene,
metrics_motif=df_metrics_motif
)
df_MR
接下来我们需要导入之前计算好的Peaks（adata_CP.h5ad）和K-mer（adata_PM.h5ad）的embedding，用于关键转录因子直接调节的基因的鉴定。
adata_CP = si.read_h5ad(os.path.join(workdir,'adata_CP.h5ad'))
adata_PM = si.read_h5ad(os.path.join(workdir,'adata_PM.h5ad'))
#make sure the indices of motifs are the same as in `motifs_genes`
adata_PM.var.index = 'M_' + adata_PM.var.index
为了减小计算量，这里我们手动指定用于后续计算直接调节的基因的关键转录因子，并计算其直接调节的基因。当然你要是对所有转录因子都感兴趣，你也可以都跑一遍。
计算得到的dict_tf_targets是一个字典，键指的是关键转录因子的名字，而值则是这个关键转录因子直接调节的下游基因，且包含这些基因的顺序与显著性得分情况。
master_regulators = ['Lef1', 'Hoxc13', 'Relb', 'Gata6', 'Nfatc1']
list_tf_motif = motifs_genes[np.isin(motifs_genes['gene'], master_regulators)]['motif'].tolist()
list_tf_gene = motifs_genes[np.isin(motifs_genes['gene'], master_regulators)]['gene'].tolist()
# get the overlap between peaks and gene loci
# the field `uns['gene_scores']['overlap']` will be added to adata_CP
_ = si.tl.gene_scores(adata_CP,genome='mm10',use_gene_weigt=True,use_top_pcs=False)
dict_tf_targets = si.tl.find_target_genes(adata_all,
adata_PM,
list_tf_motif=list_tf_motif,
list_tf_gene=list_tf_gene,
adata_CP=adata_CP,
cutoff_gene=5000)
我们也可以手动查看关键转录因子直接调节的下游基因。
dict_tf_targets['M_ENSMUSG00000001655_LINE1151_Hoxc13_D']
总结
本文介绍了SIMBA在解析基因调控网络中的应用，提供了实际操作的代码和相应的结果。作为小编，我很喜欢这篇文章用于推断基因调控网络的方法，我个人认为这篇文章的方法非常的直观，从前到后的逻辑非常的make sense，有很高的参考意义，值得大家学习。
本文为《系列教程：SIMBA-单细胞多特征统一embedding》的第五篇文章，介绍了SIMBA在推断基因调控网络中的应用，在下一篇推送中，我们将介绍SIMBA在消除批次效应中的应用情况。
https://mp.weixin.qq.com/s/fEXomMV9_Q2JiXNwMvidFg