一种应用于单细胞RNA测序数据的稀有类型细胞检测算法

马 安 琪(复旦大学计算机科学技术学院智能信息处理重点实验室 上海 200433)

0 引 言

近年来，单细胞RNA测序技术得到了迅猛的发展。在这项技术之前，RNA测序使用的是批量测序技术，只能够得到基因在组织样本中所有细胞表达量的平均值，但无法得到更细粒度的信息。与此相比，单细胞RNA测序技术是细胞层级上的测序技术，测量的是基因在每个细胞上的表达量值。因此，单细胞RNA测序技术能够进一步揭示细胞之间的异质性[1]。针对单细胞RNA测序数据的研究也越来越多，并被广泛地应用于免疫学、神经生物学、干细胞和肿瘤研究等多个领域[2]。

稀有类型细胞检测是单细胞RNA测序数据分析中的重要课题之一。其目标是通过分析单细胞RNA测序数据，找到样本组织中所属类别规模占比极少量的稀有类型细胞。尽管这些稀有类型细胞数量很少，但它们并不是可以忽略的。这些稀有类型细胞往往扮演着重要的角色，例如，抗原特异性T细胞和肿瘤干细胞都是稀有类型细胞。抗原特异性T细胞对免疫记忆的形成起关键作用[3]。而干细胞能够替换受损细胞，在治疗帕金森疾病、心脏病和糖尿病等方面有着重要意义[4]。因此，稀有类型细胞检测算法在生物上有着很强的应用需求。

对于生物专家而言，检测稀有类型细胞需要通过实验手段或者基于某些已有的先验知识。例如，已知某个稀有类型细胞的标志性基因，通过找到在这些基因上具有较高表达量的细胞来得到这一稀有类型的细胞。但并不是在所有应用场景中，都具备这种领域性的先验知识。例如，样本中可能包含没有被发现过的新稀有细胞类型，先验知识在这种情况下就不再适用。

在不需要先验知识的算法中，通过K-means[5]、层次聚类[6]等算法先对所有细胞进行聚类，再找到其中规模较小的类作为稀有细胞，是一种较为直观的方法。例如，Giniclust[7]和RaceID[8]都采用了聚类的方法来检测稀有细胞。但是，这些方法都依赖于底层的聚类算法，而大多数聚类算法在类分布空间大小差别悬殊且空间规则性差时都表现不佳。另外，聚类算法需要计算细胞两两之间的距离，因此，随着数据规模的增大，基于聚类的算法具有较高的时间花费，而现在单细胞测序数据的趋势是对越来越多的细胞进行测序。由于目标只是检测出稀有细胞，并不需要对所有细胞进行类别划分，聚类算法的代价显得过于昂贵。因此，寄希望于寻求一种有效的非聚类算法。

已有的一种非聚类思路是将稀有类型细胞看作是数据中的离群点，并使用离群点检测算法，如LOF[9]。但由于离群点和稀有类型细胞并不完全重合，这类方法在大多数应用场景之下，并不是最佳的方案。FiRE[10]是一种专门针对测序数据的稀有类型细胞检测的算法。它基于哈希的思想，快速辨别稀有细胞。但FiRE的精确率并不高。本文受FiRE基本思想的启发，提出了一种改进的稀有类型细胞检测算法。与FiRE不同，该方法用更为简单快速的多轮随机划分代替了FiRE的哈希过程，同时增加了基于最近邻调整预测结果的过程，提高了结果的精确率。

1 数据说明

1.1 数据格式

算法应用于单细胞RNA测序的基因表达量矩阵数据。数据的基本格式如图1所示，为了方便说明，这里只截取了数据的一部分。

整个数据是一个数值矩阵，在此格式的数据中，每一行对应一个基因，每一列对应一个细胞。行名是基因的名字，列名是每个细胞对应的标识码，每个标识码由“A”“T”“C”和“G”四种字符组成。标识码具有唯一性，每个标识码都和某个细胞样本相对应。矩阵中的数值表示了对应基因在对应细胞中测到的转录子的个数。

对单细胞RNA测序而言，表达量矩阵数据通常十分稀疏，大部分的矩阵项为零。另外，数据往往存在着一些噪声。

1.2 数据预处理

在执行算法之前，需要对数据进行一些预处理。预处理的目的主要是去除测序过程中因技术原因造成的一些异常数据，并对数据进行标准化处理。具体操作如下：

(1) 只保留至少在3个细胞上转录子数量达到2的基因。

(2) 计算每个细胞对应的文库大小，即列和，并将每一列除以文库大小，然后乘上文库大小的中位数。

(3) 对所有表达量值取自然对数。

(4) 根据标准化后的Fano系数(方差除以均值)挑选出在不同细胞中表达量最多变的前1 000个基因。

2 算法流程

本文提出的稀有类型细胞检测算法整体流程如下。首先，该算法在预处理之后的单细胞RNA测序数据上，进行多轮随机划分；然后综合随机划分的结果，对每个细胞的稀有程度进行打分，并初步预测出稀有类型细胞；最后，算法根据最近邻对结果进行调整，得到最终预测的稀有类型细胞结果。

2.1 随机划分

随机划分过程如下：

Step1随机抽取一个未选择过的基因，并在基因表达量的最小值和最大值之间以均匀概率随机生成一个阈值。

Step2将未分类细胞中在该基因上表达量大于等于这个阈值的细胞分为一类。

Step3重复Step 1-Step 2，直到所有细胞分类完成，或者选取的基因达到了Gmax个，停止这个过程，并将剩余细胞分为一类。

随机划分过程的伪代码如算法1所示。

算法1快速划分

输入：所有基因集合G，所有细胞集合C,表达量矩阵E。

1 已选基因数=0；

2 while已选基因数

3 gi←从G中随机抽取一个基因；

4 已选基因数+=1；

5 G←G-{gi}；

6 threshold_i<-random(Emin,Emax)；

7 Ci←{c:E[gi][c]>=threshold_i}；

8 将Ci中的细胞标记为一类；

9 C←C-Ci；

10 if C为空:

11 结束；

12 end if

13 end while

14 if还有未被分类的细胞:

15 将剩余细胞标记为一类；

16 end if

输出：细胞随机划分结果。

2.2 打 分

重复上一步随机划分过程L轮，根据这L轮的分类结果，对细胞的稀有程度进行打分。本文认为在随机划分情况下，与稀有类型细胞划分在同一类的可能性是比较小的，用plx表示在第l轮中与细胞x分在一起的概率，表示为：

(1)

式中：n代表总细胞数，labelsl代表第l轮快速划分后的结果类别标签。综合L轮的随机划分结果，细胞x的稀有程度分数为：

(2)

该分数越高，则认为细胞的稀有程度越高。

2.3 预 测

在得到连续的分值之后，算法还希望预测出二值化的标签，即预测出一个细胞是否是稀有类型细胞。算法采用和FiRE相同的判定策略，如果满足式(3)，则把x初步标记为稀有类型细胞，否则标记为非稀有类型细胞。

scorex≥(75%分位数)+1.5×(25%分位数)

(3)

2.4 基于最近邻算法调整

经过以上几步，算法已经得到了初步的预测结果。但在初步的结果中，往往包含一些假阳性的噪声点。为了去除掉结果中大部分的假阳性，接下来，算法基于最近邻算法对结果进行调整，以此提高结果的精确率。由于稀有类型细胞并不是一些孤立的离群点，而是按类聚集的。本文认为，如果一个细胞是稀有类型细胞，它最近的邻居也应该是稀有类型细胞。

在求最近邻前，算法先通过PCA[12]降维算法，把数据降至20维。接下来，对上一步中每一个稀有细胞，基于欧氏距离计算降维后最近的k个邻居。算法将检查它的k个最近邻居是否都在上一步标记成稀有类型细胞，如果有至少一个不是，则重新标记为非稀有类型细胞。

3 实验与结果分析

3.1 数据集

实验使用了两个来自10X Genomics测序平台的单细胞RNA测序数据集。这两个数据集都已经根据已有的生物知识打好了细胞类型的标签，作为实验中计算各个指标的标准答案。

第一个数据集PBMC[11]是中测试检测稀有类型细胞的一个生成数据，采样自外周血单核细胞测序数据，包含CD14+Monocyte、CD56+NK和CD19+B三种类型的细胞，其中CD14+Monocyte是该采样数据中占比少数的稀有类型细胞。第二个数据集是FiRE工具包中的样例数据，包含293T、Jurkat两种类型的细胞，Jurkat是该数据集中的稀有类型细胞。

各数据集中稀有类型细胞具体情况如表1所示。

3.2 评估指标

实验中使用了三种评估指标，分别是精确率P(Precision)、召回率R(Recall)和综合前两项指标的F1值。具体计算公式如下：

(4)

(5)

(6)

式中：TP、FP和FN分别表示预测结果中真阳性、假阳性和假阴性的个数。

3.3 实验结果

本文实验统一取参数Gmax=50，L=100，k=2。对于对比算法，均使用官方推荐的默认参数。

图2是本文方法在Jurkat-293T数据集上各细胞稀有程度分数的灰度图，该图中细胞的分布根据原数据经过t-SNE[13]降维后得到，稀有程度分数越高，图中对应亮度越低。图2中下方真实稀有类型细胞所在簇的颜色都较暗，而非稀有细胞在灰度图上较亮。另外，尽管有一些非稀有类型的细胞也显示出了较暗的颜色，它们中的大部分会在基于最近邻调整的阶段被从预测结果中剔除。

图2 稀有程度分数分布的灰度图

图3是在Jurkat-293T数据集上基于最近邻调整前后预测结果的对比，其中：深黑色的点代表本文算法预测出的稀有类型细胞，浅灰色的点代表非稀有类型细胞。算法初步预测出的结果中，虽然大部分都是真实的稀有类型细胞，但也包含了一些非稀有类型细胞。而在基于最近邻调整后，结果中的大部分假阳性都被剔除了。尽管右上角仍然存在一些被误判为稀有类型细胞的点，但总体上调整后的结果和数据集真实答案基本一致。

表2是该方法与FiRE、LOF在两个单细胞RNA测序数据集上的对比结果。从整体表现来看，LOF最差，FiRE稍好一些，而本文方法是最好的。

表2 在各数据集上的对比结果

其中，FiRE和本文方法在两个数据集中的Recall都是1，也就是说，这两种方法都能找出两个数据集中所有的稀有类型细胞。而LOF会遗漏掉部分稀有类型细胞。另外，LOF和FiRE的精确率都不够高，因此导致最后的F1值也比本文方法低。在这两个单细胞RNA测序数据集的结果表明，本文方法在检测稀有细胞类型问题上比其他方法更加精确有效。

为了进一步研究稀有类型细胞所占比例对算法的影响，本文调整了Jurkat-293T数据中稀有细胞的比例，记录三种算法在不同稀有类型细胞比例数据集F1值的结果。图4表明，三种算法的效果都随着稀有细胞所占比例的下降而变差，但本文方法始终是表现最好的，且下降幅度比其他两种算法要低。即使在最低比例(0.5%)的数据集上，本文方法依然有着较高的F1值，这表明本文方法在不同稀有类型比例的数据集上，具有一定的稳定性。

图4 三种算法在不同稀有类型细胞比例数据集上的F1值

4 结 语

本文提出了一种应用于单细胞RNA测序表达量数据的稀有类型细胞检测算法。它基于多轮随机划分的结果，对每个细胞的稀有程度进行打分，然后基于最近邻对结果进行调整，得到最终每个细胞是否为稀有类型的预测标签。

这一方法不依赖于生物的领域知识，同时也避免了聚类算法中较高的时间复杂度。从单细胞RNA测序数据集上的实验结果来看，它具有较高的精确率和召回率，表现优于其他非聚类方法。尽管当稀有类型细胞比例下降时，本文方法表现也有所下降，但仍然要比其他方法在同比例的情况下要好。未来工作将考虑把本文方法和聚类相结合，先用本文提出的稀有类型检测方法预测出稀有类型细胞，然后对稀有类型细胞和非稀有类型细胞采用不同的聚类策略，以在单细胞RNA测序数据上得到更加准确的聚类结果。