SIMBA方法解读

fig1

预处理

scRNA-seq

过滤掉在少于三个细胞中表达的基因。原始计数按文库大小标准化，然后进行对数转换。可选地，可以执行HVG选择以删除非信息性基因并加速训练过程。在将特征输入限制为通过HVG选择确定的特征输入时，未观察到所得细胞嵌入的显著差异，但不会生成非可变基因的 SIMBA 嵌入，因为它们未在图中编码。

scATAC-seq

过滤掉少于三个细胞中的peak。或者，实施一个可扩展的基于截断 SVD 的程序来选择峰，作为初步步骤，以额外过滤非信息峰并加速训练过程。首先，选择前 k 个主成分 (PC)，其中 k 的选择基于方差图。然后，对于前 k 个 PC 中的每一个，使用由“kneed”实现的拐点检测算法根据载荷自动选择peak。最后，将为每个 PC 选择的峰组合起来并表示为“变量峰”。与使用 scRNA-seq 数据的观察结果类似，变量峰选择的可选步骤对生成的细胞嵌入的影响可以忽略不计。尽管对生成的嵌入的影响微乎其微，但此特征选择步骤在减少训练过程时间方面具有显着的实际优势。

使用 JASPAR2020 中的“Biostrings”和“motifmatchr”包执行 k-mer 和motif扫描。SIMBA 的实现中包含一个方便的 R 命令行脚本“scan_for_kmers_motifs.R”，它将peak列表（格式为 bed 文件）转换为稀疏的peaks-by-k mers和peaks-by-motifs矩阵，该矩阵存储为 hdf5 格式的文件。

图构建（5种场景）

scRNA-seq分析

在构建细胞和基因graph时，如果基因在给定细胞中表达，则在细胞和基因之间添加一条边。为了区分每条边的强度，提出了一种bins方法，将基因表达值分为不同的级别，同时保留原始分布。不同级别的基因表达由不同类型的关系编码。具体来说，首先使用基于 k-means 的程序近似归一化基因表达矩阵中非零值的分布。首先，将连续的非零值分箱到 n 个区间（默认情况下，n = 5）。使用一维 k 均值聚类定义箱宽，其中每个箱中的值分配给相同的聚类中心。然后将连续矩阵转换为离散矩阵，其中 1、…、n 用于表示 n 个基因表达级别。零值保留在此矩阵中。然后，通过将两种类型的实体（细胞和基因）编码为节点，将具有 n 个不同权重的关系（即 n 个基因表达级别）编码为边来构建图。这 n 个关系权重的范围从 1.0 到 5.0，步长为 5 / n，表示基因表达水平（最低：1.0，最高：5.0），因此与高表达水平相对应的边对嵌入的影响比中等或低表达水平的边更大。正如预期的那样，观察到，随着bins数量的增加，离散化分布接近原始分布。然而，表达分辨率的增加对生成的嵌入影响不大。此离散化是在 SIMBA 包中使用函数“si.tl.discretize()”实现的。

除了关系类型权重外，SIMBA 还支持在构建图时将基因表达值直接编码为边权重。此过程会生成与分箱过程类似的嵌入。这进一步表明离散化bins在捕获生物信息方面是有效的。这种对边权重的支持是在 SIMBA 包中使用函数“si.tl.gen_graph(add_edge_weights=True)”实现的。

scATAC-seq分析

peak-by-cell矩阵被二值化：“1”表示峰内至少有一个read，否则分配“0”。该图是通过将两种类型的实体（细胞和峰）编码为节点，将它们之间的关系（表示给定峰在细胞中的存在）编码为边来构建的。单个关系类型的权重为 1.0。当 DNA 序列特征可用时，它们被使用 k-mer 和motif实体作为节点编码到图中。这是通过首先将peak-by-k mer或peak-by-motif矩阵二值化，然后使用peak、k-mer 和motif作为节点，并使用peak内这些实体的存在作为这些额外节点和峰节点之间的边来构建原始peak-by-cell图的扩展。k-mer 和峰之间的关系被分配了 0.02 的权重，而 TF 基序之间的关系被分配了 0.2 的权重。值得注意的是，根据具体的分析任务，k-mers 和motif可以彼此独立地用作图的节点输入。

多模态分析

将上述使用 scRNA-seq 和 scATAC-seq 数据构建图的策略结合起来，构建了多组学图。

批次整合

按照“scRNA-seq分析”中所述构建每个批次的图。通过基于截断随机 SVD 的程序推断不同批次细胞之间的边缘，以链接不同批次的不相交图。更具体地说，在 scRNA 序列数据的情况下，考虑两个基因表达矩阵 $X1_{n_{1}\times m}$ 和 $X2_{n_{2}\times m}$ ，其中 $n_{1}$ 和 $n_{2}$ 分别是两个批次的细胞数量， $m$ 是gene数量。

然后计算： $X=X1\times X2^{T}$ 随后对 $X$ 执行截断随机 SVD： $X=U\times \Sigma\times V^{T}$ 其中， $U$ 是 $n_{1}\times d$ 的矩阵， $\Sigma$ 是 $d\times d$ 的矩阵， $V$ 是 $n_{2}\times d$ 的矩阵，默认 $d = 20$ 。

$U$ 和 $V$ 都进一步进行了 L2 归一化。对于 $U$ 中的每个细胞，我们在 $V$ 中搜索 k 个最近邻居，反之亦然（默认情况下，k = 20）。最终，只有 $U$ 和 $V$ 之间的相互最近邻居被保留为细胞之间的边（注意是推断的边）。推断不同批次细胞之间的边的过程在 SIMBA 包中的函数“si.tl.infer_edges()”中实现。

对于多个批次，SIMBA 可以灵活地推断任意一对batch-pair之间的边。然而，在实践中，边是在最大的数据集或包含最完整预期细胞类型集的数据集与其他数据集之间推断的。

多模态整合

scRNA-seq 和 scATAC-seq 图分别按照“scRNA-seq 分析”和“scATAC-seq 分析”中的步骤构建。为了推断 scRNA-seq 和 scATAC-seq 细胞之间的边，首先计算 scATAC-seq 数据的基因活性分数（gene activity score）。更具体地说，对于每个基因，考虑 TSS （转录起始位点）上游和下游 100 kb 内的peak。与基因体区域重叠或在基因体上游 5 kb 内的peak的权重为 1.0。否则，使用指数衰减函数根据peak value与 TSS 的距离对其进行加权： $exp(\frac{-distance}{5000})$ 。随后，将每个基因的gene score计算为所考虑峰值的加权和。然后将这些基因得分缩放到相应的基因大小。这些步骤由 SIMBA 中的函数“si.tl.gene_scores()”实现。为了方便用户，SIMBA 包整理了几个常用参考基因组的基因注释，包括 hg19、hg38、mm9 和 mm10。一旦获得基因得分，就执行“批次整合”中描述的相同程序，使用 SIMBA 中的函数“si.tl.infer_edges()”推断 scRNA-seq 和 scATAC-seq 分析的细胞之间的边。

生成图的过程在 SIMBA 包中的函数“si.tl.gen_graph()”中实现。

图学习

在构建生物实体之间的多关系图之后，作者采用了知识图谱和推荐系统中的图嵌入技术来为这些实体构建无监督表示。

提供一个input无向图 $G = (V, E)$ ，其中 $V$ 是一组实体（节点） $E$ 是一组边，在源实体 $u$ 和目标实体 $v$ 之间存在通用边 $e = (u, v)$ 。进一步假设每个实体都有不同的已知类型（例如，细胞或peak）。

图嵌入方法通过随机梯度下降优化edge预测目标，为每个 $v \in V$ 学习一个 $D$ 维嵌入向量，其中实验中使用 $D = 50$ 。实体 $v$ 的embedding记为 $\theta_{v}$ 。

对于edge $e = (u, v)$ ，记 $s_{e}=\theta_{u}\cdot\theta_{v}$ 为 $e$ 的得分，损失为： $L_{e}=-log\frac{exp(s_{e})}{\sum_{e'\in N}exp(s_{e'})}w_{e}$ 其中， $N$ 是通过破坏 $e$ 生成的一组“负样本”候选边， $w_e$ 是边权重，默认情况下是关系权重，但在每种关系类型中可能因边而异。例如，细胞和基因之间的边可以编码为具有不同边权重的单一关系，这些边权重编码标准化的基因表达水平（见“scRNA-seq分析”）。

通过将目标边 $e = (u, v)$ 中的源实体或目标实体替换为随机采样的实体来构建负样本。因此，例如，对于cell-peak的边，仅对cell和peak实体之间的负候选样本进行采样。这种设置至关重要，因为大多数随机选择的边是无效的（例如，峰-峰）。

使用了 PyTorch-BigGraph 框架，该框架可以高效计算多种实体类型的多关系图嵌入，并且可以扩展到包含数百万或数十亿个实体的图。对于 130 万个细胞，PyTorch-BigGraph 训练本身仅需大约 1.5 小时，使用 12 个 CPU 核心，无需 GPU。

SIMBA空间中查询实体

信息丰富的 SIMBA 嵌入空间可用作实体（包括细胞和特征）的数据库。为了在“SIMBA 数据库”中查询给定细胞或特征的邻近实体，我们首先根据其 SIMBA 嵌入构建所有实体的 k-d 树。然后，使用欧几里得距离在树中搜索最近的邻居。为此，SIMBA 查询可以在指定半径内执行 k 最近邻居 (KNN) 或最近邻居搜索。SIMBA 还提供了将搜索限制为某些类型实体的选项，当某种类型的实体数量远远超过其他实体时，这很有用。例如，给定细胞的 k 个最近特征可能都是峰值，而基因是感兴趣的特征。在这种情况下，SIMBA 允许用户添加“过滤器”以确保在指定类型的实体内执行最近邻居搜索。此过程在函数“st.tl.query()”中实现，其可视化在 SIMBA 包中的函数“st.pl.query()”中实现。

识别TF-target genes

为了推断给定主调节因子的靶基因，我们假定，在共享的 SIMBA 嵌入空间中，（1）靶基因靠近 TF 基序和 TF 基因，表明靶基因的表达与 TF 的表达和 TF 基序的可及性高度相关，并且以细胞类型特异性的方式呈现；（2）靶基因位点附近的可及区域（峰）必须靠近 TF 基序和靶 TF 基因，表明靶基因位点附近的顺式调控元件的可及性与 TF 的表达和 TF 基序的可及性高度相关，并且以细胞类型特异性的方式呈现。

给定一个主调节因子，通过比较 SIMBA 共嵌入空间中 TF 基因、TF 基序和候选靶基因基因组位点附近的峰的位置来识别其靶基因。

更具体地说，我们首先分别搜索该主调节因子的基序（TF 基序）和基因（TF 基因）周围的 k 个最近邻基因（默认 k = 200）。这些邻居基因的并集就是初始的候选靶基因集。然后根据以下标准对这些基因进行筛选：假定靶基因 TSS 上游和下游 100 kb 内的开放区域（峰）必须包含 TF 基序。

接下来，对于每个候选靶基因，我们计算了 SIMBA 嵌入空间中的四种距离：(1) 候选靶基因与 TF 基因的嵌入之间的距离；(2) 候选靶基因与 TF 基序的嵌入之间的距离；(3) 候选靶基因与 TF 基序的基因组位点附近的峰之间的距离；以及 (4) 候选靶基因与其基因组位点附近的峰之间的距离。所有距离（默认为欧几里得距离）都转换为所有基因或所有峰之间的等级，以使距离在不同的主调节器之间具有可比性。

最终的靶基因列表由计算出的排名决定，使用两个标准：（1）TF 基因或 TF 基序最近的峰值中至少有一个在预定范围内；（2）候选靶基因的平均排名在预定范围内。此过程在 SIMBA 中的函数“st.tl. find_target_genes ()”中实现。