在生信分析中，算法的开发和创新是一个非常重要的环节。随着单细胞技术的飞速发展，现有的一些分析算法存在的不足已经日渐暴露出来了。就在近日，nature communications（14.913）上新发表了一篇关于单细胞数据聚类的新算法，就让我们一起来看看这个新的算法具体又是怎样创新的吧。

DUBStepR——一种用于单细胞数据准确聚类的可扩展的基于相关性的特征选择方法

摘要

特征选择，也就是标记基因的选择，被广泛认为是提高聚类准确度的关键步骤。而目前现有的方法都未考虑到基因相关性的影响，因此本文提出了一种利用基因之间的相关性的新的特征选择算法——DUBStepR（Determining the Underlying Basis using Stepwise Regression）。它不但可以鉴定出类风湿关节炎患者PBMC中与疾病相关的常见和罕见的细胞类型，也可以应用于大型的数据集中，并且还可能可以直接应用于其他单细胞数据类型中，如scATAC-seq。

背景介绍

scRNA-seq数据的无监督聚类过程通常包括质控、标准化、特征选择、主成分分析（PCA）、基于距离聚类以及细胞类型定义。尽管在这个过程中特征选择是一个很关键的步骤，但是目前而言，对于scRNA-seq中的特征选择的方法研究还很少。其中一个比较好的特征选择算法就是选择细胞类型特异性（DE）基因作为特征，但这种算法还需要优化不同细胞类群之间在生物学上的不同特征的区分。目前使用最多的特征选择方法是均值方差建模（mean-variance modeling）。本文中提到了M3Drop和GiniClust两种方法，但它们都未将基因表达之间的相关性考虑进去。

基于以上，研究人员提出了DUBStepR这种基于基因之间相关性的特征选择算法，具体流程见图1中。同时，为了评估这种新的算法的适用性和性能，研究人员纳入了四种不同类型的数据同时与其他7种常用的特征选择算法进行了比较。最后还将此算法应用于类风湿关节炎患者PBMC中的单细胞ATAC-seq数据上。

PBMC，也就是外周血单核细胞(Peripheral blood monoculear cell)，起源于位于骨髓(bone marrow)的造血干细胞( hematopoietic stem cells, HSCs)，主要包括淋巴细胞(T细胞，B细胞和自然杀伤(NK)细胞)和单核细胞。通常情况下，淋巴细胞约占70 – 90%，单核细胞约占10 – 30%，不同个体之间这些细胞群的比例可能存在差异。

图1 DUBStepR流程图

结果

1. 利用基因之间的相关性预测细胞类型特异性DE基因

DUBStepR的第一步是根据细胞类型特异性DE基因（标记基因）的已知特征选择一组初始候选特征。相同细胞类型的特异性DE基因往往是高度正相关，而不同细胞类型的特异性DE基因可能是负相关的（图2a，b）。相比之下，非DE基因之间可能只是弱相关性（图2c）。因此研究人员假设，由基因最强正相关系数和最强负相关系数之间的差异得出的相关范围分数在DE基因中会显著升高。事实上，DE基因的相关范围得分确实显著高于非DE基因（图2d）。

图2 DE基因的表达相关性

2. 逐步回归法（Stepwise regression）识别最小冗余特征子集

在基因-基因相关性（GGC）矩阵中可以观察到大小不一的相关性模块（图3a），每个模块可能代表着细胞中不同的差异表达模式。为了确保候选特征集中不同细胞类型的标记更均匀的分布，通过在GGC矩阵中进行回归算法（图3b-d），识别出最具有代表性的最小冗余子集（“seed”基因）。利用这种方法选择具有不同细胞类型特异性标记的seed基因（图3e-h），然后通过逐步回归树状图的拐点来确定最佳步骤数，也就是seed基因集的大小（图3i，j）。

图3 识别最小冗余特征子集的逐步回归

3. 扩展特征集（Guilt-by-association）

基于上述方法选择的特征集虽然特征表现很显著，但是每个特征相关的基因只有少数几个基因（大多情况下是2-5个基因），这是远远不够的。因此研究人员提出了通过Guilt-by-association的方法来扩展特征集，直到DUBStepR达到特征基因的最佳数量。

4. 基准测试

基准测试是指通过设计合理的测试方法，选用合适的测试工具和被测系统，实现对某个特定目标场景的某项性能指标进行定量的和可对比的测试。本文研究人员为了测量和评估DUBStepR，将其与其余7种算法在scRNA-seq数据中进行了比较。每种算法都用了7个数据集，涵盖有4种不同的scRNA-seq：10x Genomics, Drop-Seq, CEL-Seq2和Smart-Seq on the Fluidigm C1。为了评估所选择的特征的质量，这里采用的评估指标是SI（Silhouette index），该指数指的是属于同一类的细胞之间的接近度相对于其他类群细胞的距离。可以发现，DUBStepR在整个特征集选择范围内明显优于其他所有方法（图4a）。此外，DUBStepR在7个数据集中的使用，有5个都是在这些算法中排第一的（图4b）。

为了优化细胞类型聚类，特征选择算法最好只选择DE基因，即细胞类型或亚型特异性基因作为特征标记。因此，研究人员量化了特征选择算法区分DE和非DE基因的能力，引入了AUROC（ROC曲线下的面积）。图4c中可以看到，DUBStepR在所有7个数据集上的AUROC值都超过了0.97，说明它能够很好的区分DE和非DE基因。也就是说，DUBStepR大大提高了我们选择细胞类型或亚型特异性基因（DE基因）来对scRNA-seq数据进行聚类的能力。

此外，研究人员还对DUBStepR进行了进一步的优化，以便能以较快的速度可以处理大型数据集。

图4 特征选择方法的基准测试

5. 通过DI（Density index，密度指数）预测最佳特征集

从图4a中可以看出，选择太少或者太多的特征基因都会导致聚类效果不佳。由于无法确定特征集的SI，所以本文研究人员定义了一个新的指标，这个指标可以大致模拟SI，且无需知道细胞类型标签；这个指标也就是密度指数（DI），即所有细胞对之间的均方根距离除以细胞与其k个临近细胞之间的平均距离。

直观地来说，当细胞聚在一起时，它们之间的距离应该最小，DI值应该是最大的（图5a，b）。可以发现DI和SI确实是呈现正相关性，并且在大致相同的特征集大小下趋向于达到它们的最大值（图5c）。同时，对于7个基准数据集中的5个，DI最高的特征集也最大化了SI（图5d）。除此之外，DI还有一个优点就是计算相对来说比较简单。在默认条件下，DUBStepR会选择最大化DI的特征集大小。

图5 密度指数

6. DUBStepR在类风湿关节炎样本中可靠地检测出罕见细胞类型和隐匿的细胞状态

前面的定量基准分析主要基于对健康供体细胞细胞系或FACS纯化细胞群中常见的细胞类型（＞所有细胞的10%）的检测。所以研究人员为了证明DUBStepR能够在复杂的原始样本中对细胞进行聚类，这里加入了对四名类风湿关节炎（RA）患者的8312个PBMC的scRNA-seq数据的分析。首先利用SingleR分离出T细胞和NK细胞群（5329细胞）（图8），然后通过DUBStepR对这些细胞进行聚类，进一步识别出了10个具有显著不同的基因表达特征的亚型（图6a）。其中还包括了四个罕见的细胞群（图6b）。总而言之，通过比较其他方法，DUBStepR是唯一一个可以可靠地检测到淋巴细胞群中常见和罕见的细胞类型和亚型的方法。

图6 类风湿关节炎患者外周血淋巴细胞群分析

7. DUBStepR应用于scATAC-seq数据

常用的特征选择通常不适用于scATAC-seq数据上，但是特征相关性却是适用的。研究人员进一步将新的算法DUBStepR推广应用于人骨髓细胞的scATAC-seq数据上，结果发现DUBStepR的结果更清楚地揭示了造血干细胞的三大谱系的出现：淋巴系、髓系和巨核/红细胞系（图7）。具体来说，使用Monocle 3分析得出的拓扑结构仅在DUBStepR中与已报道的造血分化层次匹配（图7d-f）

图7 已发表文献中的scATAC-seq数据特征选择的比较

图8 SingleR scores热图

讨论

DUBStepR是一种基于细胞类型特异性标记基因往往彼此之间存在比较好的相关性的算法。随着单细胞技术的飞速发展，使得可扩展性成为新型单细胞算法的一个基本特征，而DUBStepR的扩展范围是比较大的，其中的主要因素是一旦构建了基因-基因相关矩阵，下游步骤的时间和内存复杂性相对于细胞数量是恒定的。本研究中定义了一个特征空间中细胞的不均匀性或“聚集性”的度量，也就是密度指数（DI）。DI的算法相对来说也更简单。有趣的是，尽管DUBStepR不是专门设计用于检测稀有细胞类型的，但是它的检测效果要明显优于其他算法。最后，研究人员还推测DUBStep可能也可以用于scATAC-seq、scChIP-seq和单细胞甲基组测序等数据中。

本文所开发的算法已经写成了R包，源代码也提供在了GitHub上可供我们获取。此外，所采用的分析数据也是可以直接下载或申请，那我们何不赶紧动动手去用一下这新的算法呢？