scRNA-seq支持在单细胞分辨率呈现基因表达谱，这提高了对已知细胞类型的检测，以及异质细胞与疾病失调的理解。然而，scRNA-seq技术的广泛应用产生了许多庞大而复杂的数据集，这给集成来自不同批次和平台的数据集带来了计算挑战。
来自：Batch alignment of single-cell transcriptomics data using deep metric learning
目录
背景
scRNA-seq数据分析的基本任务是将细胞聚类成不同簇，然后作为候选细胞类型。对于单源数据集来说，该任务可能比较简单。但对于多源数据集来说，由于批次效应的特性，尤其是对于一些小簇的检测，该任务就非常困难。
尽管已经开发了多种方法来消除scRNA-seq分析中的批次效应，但是大多数方法旨在消除embedding空间中的批次效应，而没有考虑数据集中的聚类结构或局部结构。Seurat和MNN等常用方法都是采用相互最近邻的方法来消除批次效应，MNN一次只能分析两个批次，因此其性能受到批次校正顺序的影响。并且当批次数量增加时，MNN很容易变得不可行。另外两种方法Scanorama 和BBKNN 也在降维空间中搜索MNN，并以相似性加权的方式使用它们来指导批次集成。此外，开发了两种监督MNN方法（SMNN，iSMNN ）用于scRNA-seq的批次校正，但这两种方法要求不同批次之间的细胞类型完全相同。Liger 的目标是使用集成非负矩阵分解来消除技术变异，但其过程需要选择参考数据集。
大多数现有方法首先去除批次效应，然后对细胞进行聚类。然而，这些方法的缺点在于去除批次效应可能导致原始稀有细胞类型信息的丢失。因此，作者从原始数据的先验聚类信息开始，在具有triplet loss的深度度量学习框架下，利用批次内和批次间的最近邻（NN）信息，通过学习数据的低维表示来正确恢复真实细胞类型并消除批效应。最重要的是，scDML不受批次集成顺序的影响。在初始聚类中，首先对细胞进行高分辨率聚类，以保证初始聚类包含所有细微的和潜在的新细胞类型，然后提出一个合并准则来优化最终聚类数。该算法结合了图聚类和层次聚类方法的优点，通过将具有相同标签的点拉到一起，将具有不同标签的点推开，从而消除了批次效应。
方法
scDML的设计目标是对齐多批单细胞转录组数据。scDML的工作流如图1所示。在对scRNA-seq数据进行预处理（包括归一化、log1p变换、寻找高变基因、缩放数据、PCA嵌入）后，首先使用了基于图的高分辨率聚类算法。然后，利用批内和批间的k-最近邻（KNN）和互近邻（MNN）信息评估细胞簇之间的相似性，并建立了具有层次结构的对称相似性矩阵。
在给定的高度（给定的簇数）上切割树，通常会以选定的精度生成分区的簇。为了成功去除批次效应，作者通过考虑Hard triplets，这有助于学习一个低维嵌入，同时去除批次效应。
- 图1：scMDL工作流程。MNN表示相互最近的邻居，KNN表示k-最近的邻居。建立聚类相似度矩阵的过程和所提出的归并规则旨在保留原始的聚类信息，包括一些细微的生物聚类。三元组网络的目标是最小化anchor-positive pair（来自同一簇）之间的距离，同时最大化anchor-negative pair（来自不同簇）之间的距离。
预处理
预处理中需要完成五个重要任务：筛选低质量细胞和基因、细胞归一化、对数归一化、检测高变基因、z-score归一化。以上所有步骤都在scanpy中实现。
令 X X X为 n c × p nc\times p nc×p的scRNA-seq矩阵， n c nc nc为细胞数， p p p为基因数，令 x i j x_{ij} xij表示细胞 i i i的基因 j j j含量。在过滤步骤中，首先去除nGene <10的低质量细胞，然后去除nCells <3的基因。在细胞归一化中，将每个细胞的计数除以所有基因的总计数，并使用sc.pp.normalize_total函数乘以常数10000，得到归一化表达式 y i j y_{ij} yij。然后对 y i j y_{ij} yij进行对数变换，以及高变基因筛选，并缩放数据。
基于PCA空间进行初始聚类
令 Y H V G Y_{HVG} YHVG为预处理后的 n c × p H V G nc\times p_{HVG} nc×pHVG矩阵，仅包含 p H V G p_{HVG} pHVG个高变基因。在 Y H V G Y_{HVG} YHVG上执行PCA得到 X e m b X_{emb} Xemb（主成分 n p c a = 100 n_{pca}=100 npca=100）。然后scDML在PCA嵌入空间上应用基于图的聚类方法Louvain方法得到初始化的聚类结果。这个过程在scanpy中，更高的分辨率意味着找到更多和更小的簇。
然而，真实的细胞类型的数量在真实数据中通常是未知的，因此如何为Louvain算法找到正确的分辨率是一个挑战。scDML在Louvain算法中设置了一个相对较大的分辨率（默认为3.0），这可能有助于在数据集中发现更微妙的细胞类型。
在PCA空间找到NN pairs
为了合并从Louvain算法获得的初始化聚类，我们需要计算聚类之间的相似度。scDML首先通过在批内和批间寻找 NN（最近邻，nearest neighbor）pair 来构建所有簇之间的联合图（joint graph）。
寻找KNN pairs intra batch（批次内）
令 X e m b = ( X 1 , . . , X M ) X_{emb}=(X^{1},..,X^{M}) Xemb=(X1,..,XM)为 M M M个批次PCA数据的连接，形状为 n × n p c a n\times n_{pca} n×npca。记 S k S_{k} Sk为第 k k k批次中所有KNN pair的集合，细胞 i i i和细胞 j j j形成一个KNN pair的情况如下：
- i , j ∈ B k i,j\in B_{k} i,j∈Bk， B k B_{k} Bk为batch k k k的所有细胞；
- “ ( i , j ) ∈ S k (i,j)\in S_{k} (i,j)∈Sk以及 ( j , i ) ∈ S k (j,i)\in S_{k} (j,i)∈Sk” 等价于 “ i ∈ K N N ( x j k ) i\in KNN(x_{j}^{k}) i∈KNN(xjk)或者 j ∈ K N N ( x i k ) j\in KNN(x_{i}^{k}) j∈KNN(xik)”，元组 ( i , j ) (i,j) (i,j)和 ( j , i ) (j,i) (j,i)都是KNN pairs。 K N N ( x j k ) KNN(x_{j}^{k}) KNN(xjk)表示第 k k k batch中细胞 i i i的 k k k个最近邻居的集合。
为了在batch内找到KNN，我们为scDML设置邻居的数量 K i n = 5 K_{in} = 5 Kin=5，并且用cosine距离作为度量。
寻找MNN pairs inter batch（批次间）
许多研究证明基于MNN（mutual nearest neighbor）的方法可以有效去除scRNA-seq数据中的批次效应，如MNN，BBKNN，Scnorama。因此，在这里作者也使用MNN pairs来构建不同批次之间的簇相似度。
为了与KNN pairs的定义对应，令 S a , b S_{a,b} Sa,b为batch a a a 和batch b b b之间的MNN pairs集合，细胞 i i i和细胞 j j j形成一个MNN pair的情况如下：
- i ∈ B a , j ∈ B b i\in B_{a},j\in B_{b} i∈Ba,j∈Bb；
- “ ( i , j ) ∈ S a , b (i,j)\in S_{a,b} (i,j)∈Sa,b以及 ( j , i ) ∈ S a , b (j,i)\in S_{a,b} (j,i)∈Sa,b” 等价于 “ i ∈ M N N ( x j b , X a ) i\in MNN(x_{j}^{b},X^{a}) i∈MNN(xjb,Xa)或者 j ∈ M N N ( x i a , X b ) j\in MNN(x_{i}^{a},X^{b}) j∈MNN(xia,Xb)”， M N N ( x i a , X b ) MNN(x_{i}^{a},X^{b}) MNN(xia,Xb)表示batch b b b中最接近batch a a a中细胞 i i i的细胞集合。
scDML默认设置邻居数 K b w = 10 K_{bw} = 10 Kbw=10，使用cosine距离计算MNN pairs。
计算簇相似度，构造分层聚类树
令 S i n S_{in} Sin表示所有KNN pairs， S b w S_{bw} Sbw表示所有MNN pairs， S S S表示所有NN pairs： S i n = S 1 ∪ S 2 ∪ . . . ∪ S M S b w = ∪ j = 1 M ∪ i = j + 1 M S i , j S = S i n ∪ S b w S_{in}=S_{1}\cup S_{2}\cup ... \cup S_{M}\\ S_{bw}=\cup_{j=1}^{M}\cup_{i=j+1}^{M}S_{i,j}\\S=S_{in}\cup S_{bw} Sin=S1∪S2∪...∪SMSbw=∪j=1M∪i=j+1MSi,jS=Sin∪Sbw令 N = ∣ S ∣ N=|S| N=∣S∣为NN pairs的总数，基于初始聚类的簇标记，和NN pairs集合，可以得到一个对称矩阵： A = ( a 1 , 1 a 1 , 2 ⋯ a 1 , C a 2 , 1 a 2 , 2 ⋯ a 2 , C ⋮ ⋮ ⋱ ⋮ a C , 1 a C , 2 ⋯ a C , C ) A= \left ( \begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,C} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,C} \\ \vdots & \vdots & \ddots & \vdots \\ a_{C,1} & a_{C,2} & \cdots & a_{C,C} \end{matrix} \right ) A= a1,1a2,1⋮aC,1a1,2a2,2⋮aC,2⋯⋯⋱⋯a1,Ca2,C⋮aC,C 其中， a i , j a_{i,j} ai,j表示簇 i i i和簇 j j j的NN pairs数量。值得注意的是，scDML删除了属于同一簇的所有NN pairs（kNN和MNN pairs），也就是说， A A A对角线上的元素为0。
显然，每个簇的规模越小，能找到的NN pairs数量就越少。所以，在用 A A A表示簇之间的相似性时，应该考虑簇大小（每个簇中的细胞数量）。为了实现上述目标，scDML采用了一种简单而直观的方法来计算簇之间的相似矩阵： S = ( s 1 , 1 s 1 , 2 ⋯ s 1 , C s 2 , 1 s 2 , 2 ⋯ s 2 , C ⋮ ⋮ ⋱ ⋮ s C , 1 s C , 2 ⋯ s C , C ) \textbf{S}= \left ( \begin{matrix} s_{1,1} & s_{1,2} & \cdots & s_{1,C} \\ s_{2,1} & s_{2,2} & \cdots & s_{2,C} \\ \vdots & \vdots & \ddots & \vdots \\ s_{C,1} & s_{C,2} & \cdots & s_{C,C} \end{matrix} \right ) S= s1,1s2,1⋮sC,1s1,2s2,2⋮sC,2⋯⋯⋱⋯s1,Cs2,C⋮sC,C 其中： s i , j = a i , j m i n ( m i , m j ) , i , j = 1 , . . . , C s_{i,j}=\frac{a_{i,j}}{min(m_{i},m_{j})},\thinspace i,j=1,...,C si,j=min(mi,mj)ai,j,i,j=1,...,C其中， m i m_{i} mi表示簇 i i i的细胞数，根据 S \textbf{S} S可以构造一个新的graph，节点为 C C C个簇，再次聚类得到新的 K K K个簇（找到graph中所有连通的部分，重新分配聚类的标签）。
度量学习移除批次效应
MNN pairs可以在一定程度上消除批次效应。因此，可以使用上述MNN引导的信息来构建更好的低维嵌入。此外，度量学习是一种直接基于距离度量的方法，旨在从弱监督数据自动构建特定任务。
虽然scDML重新分配了数据集的聚类标签，但原始数据实际上并没有对批处理效应进行校正。在这里，作者使用了具有triplet loss的度量学习（DML，deep metric learning）方法来捕获更准确的低维表示。从广义上讲，目标是学习一种距离度量，将标签相同的点拉近，同时将标签不同的点推远，同时考虑批处理效应的影响。
triplet定义
为了充分利用聚类信息，根据以下准则构造三元组（anchor，positive，negative）。给定一个细胞 a a a（anchor），从与 a a a标签相同的细胞簇中随机选择一个细胞 p p p作为positive，从与 a a a标签不同的细胞簇中随机选择一个细胞 n n n作为negative， ( a , p , n ) (a,p,n) (a,p,n)可视为一个三元组。数据集中的任何细胞都可以用作anchor。
triplet loss
三元组损失定义如下： L ( a , p , n ) = m a x ( d ( a , p ) − d ( a , n ) + m , 0 ) L(a,p,n)=max(d(a,p)-d(a,n)+m,0) L(a,p,n)=max(d(a,p)−d(a,n)+m,0)其中， d d d是距离度量，采用欧几里得距离。 m m m是边界参数，默认为0.2。通过使anchor-positive pairs之间的距离最小化和anchor-negative pairs之间的距离最大化来优化。根据三元组损失的定义，三元组有三种可能的类别：Easy triplets（损失为0的三元组loss），Hard triplets（负点比正点更接近anchor点的三元组loss），Semi-hard triplets（负点并不比正点更接近anchor点，但loss值仍然是正数）。考虑到时间和内存消耗，作者选择嵌入空间中的Hard triplets来训练scDML模型。
scDML设计
scDML采用简单的神经网络，输入层节点数为1000，隐含层节点数为256，嵌入层节点数为32。根据三元组的定义，anchor点是独立的，因此可以用小批量策略优化三元组损失。scDML的实现基于PyTorch框架，该框架充分利用了一个名为pytorch_metric_learning的可伸缩包。
指标
作者通过常用的UMAP方法将scDML和所有比较方法的结果可视化。此外，采用三个指标（ARI 、NMI 、ASW_celltype）来评估聚类性能，并采用另外三个指标（iLISI 、BatchKL 、ASW_batch）来评估去除批次效应的能力。
ARI（Adjusted rand index）用于量化聚类精度，可以通过python模块sklearn.metrics.cluster中的adjuststed_rand_score函数计算。ARI度量两个聚类结果之间的相似性。通过计算ARI，将整合后数据的聚类结果与预定义的细胞类型进行比较（比如，比较同时位于cluster A和ground truth class B的细胞数）。ARI的取值范围为[0,1]，值越大表示相似度越高。
NMI（Normalized mutual information）用于度量聚类精度，通过sklearn.metrics.cluster中的函数normalized_mutual_info_score计算。NMI的输入与ARI相同，NMI在[0,1]范围内且值越高，表明聚类结果与真实细胞类型的相似性越高。
ASW_celltype，细胞类型的平均轮廓宽度（Average silhouette width for cell type）用于评估细胞类型的纯度，可以通过sklearn.metrics中的函数silhouette te_score计算。细胞 i i i的细胞类型标签的轮廓宽度定义为： s i = b i − a i m a x { a i , b i } s_{i}=\frac{b_{i}-a_{i}}{max\left\{a_{i},b_{i}\right\}} si=max{ai,bi}bi−ai其中， a i a_{i} ai是细胞 i i i到所有相同标签细胞的平均距离， b i b_{i} bi是细胞 i i i到所有不同标签细胞的最小距离。ASW_celltype是所有细胞轮廓宽度的平均值（范围[0,1]），其中较高的值表示细胞更接近具有相同标签的细胞，而远离具有不同标签的细胞。作者基于低维（默认维度为32）PCA嵌入空间中定义的细胞来计算ASW_celltype。
ASW_batch，批次的平均轮廓宽度（Average silhouette width for batch）用于评估批次的全局混合程度。一个批次中的细胞 i i i的轮廓宽度定义为： s i = b i − a i m a x { a i , b i } s_{i}=\frac{b_{i}-a_{i}}{max\left\{a_{i},b_{i}\right\}} si=max{ai,bi}bi−ai其中， a i a_{i} ai是从细胞 i i i到同一批次的所有细胞的平均距离， b i b_{i} bi是细胞 i i i到分配给其他批次的同簇细胞的最小平均距离。对于ASW_batch，较高的值表示细胞更接近同一批次的细胞，而远离不同批次的细胞。较高的ASW_batch意味着批次之间的互斥性更强，相反，较低的ASW_batch通常意味着更好的批次混合和批次效应校正。ASW_batch的范围也在[0,1]内。
iLISI，逆辛普森整合指数（Inverse Simpson’s index of integration）用于评估整合后批次的局部混合程度。局部逆辛普森指数（LISI，local inverse Simpson’s index）可用于基于预选困惑度选择的局部邻居来测量批分布（iLISI）。使用选定的细胞邻居，然后在批次标签上为iLISI指数计算LISI，接近预期批次数的分数表示批次混合良好。使用R包lisi中的函数compute_lisi为数据集中的所有细胞计算iLISI，输出平均分数，较高的iLISI表示批次混合性能更好。
BatchKL，KL散度用于评估基于嵌入表示的批次效应去除方法的性能，值越低表示混合性能越好。
