大规模单细胞转录组测序数据的聚类方法比较*

朱晓姝,蒙 霜,龙法宁

(1.广西师范大学计算机科学与工程学院,广西桂林 541004;2.玉林师范学院计算机科学与工程学院,广西玉林 537000)

单细胞转录组测序(single-cell RNA-sequencing,scRNA-seq)技术对单个细胞进行测序,可以准确度量每个细胞的基因表达水平,更清晰地反映它们之间的差异[1,2]。该技术解决了批量细胞(Bulk cell)转录组测序技术对多个细胞测序获得多个细胞的基因表达平均水平时,容易丢失单个细胞独有信息的问题[3,4]。以scRNA-seq数据为基础,分析细胞异质性[5,6],刻画基因表达的动态变化[7],对细胞聚类识别细胞类型[8],可以在细胞发育和细胞分化、疾病早期诊断和预后等精准医疗领域发挥重要的作用[9]。例如,阿尔茨海默症[9,10]、癌症等重大疾病的早期诊断[9,11],对延长病人存活时间、降低家庭和社会负担具有重要的意义[11,12]。

当前,scRNA-seq技术发展迅速,可以对数万个单细胞测序,其样本规模从以前的几十至几百个细胞增加到几千至几万个细胞,导致计算复杂度极大增加[13]。此外,scRNA-seq数据呈现高稀疏性[14]、高噪声[15]、高维度[16]、结构信息和位置信息缺乏等特点,对单细胞准确聚类造成了较大的困难。高稀疏指“0”值占比达65%-95%,是极度稀疏的数据。高噪声指单细胞分离时产生低质量细胞、单细胞扩增时覆盖度不均匀,以及低的测序深度可能导致基因低表达,而且不同测序平台、测序协议和参数得到的测序值范围差异较大,这些都会导致大量的技术噪声。高维度指数据维度超过10 000维,难以准确地度量细胞间相似性,并增加计算开销。结构信息和位置信息缺乏指测序时分离了每个细胞,导致细胞间关联等结构信息[17]、细胞的位置信息丢失,从而降低聚类准确性和鲁棒性[6,18]。当前,单细胞聚类方法包括传统的聚类方法和专门设计的方法[19,20],主要有k-均值聚类(k-means clustering)和层次聚类(Hierarchical Clustering,HC)等经典聚类方法[21,22],以及基于映射[23,24]、基于图划分[25]、基于密度[26,27]、基于集成的单细胞聚类方法[2]。这些方法在“1.3聚类方法”小节中有具体的描述和分析。

对于不同类型、不同规模的scRNA-seq数据,不同的聚类方法在识别细胞类型时,其性能和结果存在较大差异[28,29]。因此,为了便于研究者根据数据特点选择合适的聚类方法,准确识别细胞类型[19,30],本研究分别对采用不同测序协议、数据格式和数据规模的14个scRNA-seq数据集进行分析[31,32],流程见图1。其中,测序协议包括全长测序和双端测序;数据格式包括每百万计数(Counts Per Million,CPM)、每百万转录本(Transcripts Per Million,TPM)、每百万转录本的每千碱基片段(Fragments Per Kilobase of Transcript Per Million,FPKM)、每百万映射读长的每千碱基读长(Reads Per Kilobase Per Million Mapped Reads,RPKM)以及原始读长(Reads);数据规模为124-9 519个细胞;稀疏度为64.11%-94.70%[21,33]。本研究将分析比较6种代表性基因选择方法选择高差异表达基因的情况,以及6种单细胞聚类方法的聚类准确性和鲁棒性。

图1 scRNA-seq数据聚类分析流程

1 方法比较

1.1 质量控制

对scRNA-seq数据进行质量控制,既可以降低数据维度,又可以去除噪声,从而提高单细胞聚类的性能[5,34]。质量控制主要包括两个步骤:过滤和归一化。(1)过滤。通过设置阈值,过滤低质量的细胞和基因。例如,设置某细胞中表达基因数阈值,过滤多细胞、死细胞等低质量细胞;设置某基因表达的细胞数阈值,过滤低表达基因或稀有基因。(2)归一化。通过使用归一化因子或对数转换对基因表达量进行归一化,可以消除scRNA-seq数据的拖尾现象。表1列出了6种典型的基因选择方法的质量控制分析。

表1 6种典型的基因选择方法的质量控制分析

从表1可以看出,6种基因选择方法在质量控制的过滤和归一化中分别设置了不同的阈值参数。例如,过滤包含表达基因数少的低质量细胞时,阈值参数分别设为非0表达基因数、基因表达量总和,以及线粒体基因的占比;过滤在所有细胞中低表达的稀有基因时,阈值参数分别设为表达的细胞数和基因表达量总和。为了过滤低质量细胞,细胞中非0表达基因数的阈值设置为200或2 000,细胞中基因表达计数总和阈值设置为3倍的绝对偏差中位数(Median Absolute Deviations,MADs)。为了过滤稀有基因,非0表达的细胞数阈值设置为10,基因表达均值阈值设置为0.05。使用size.factors因子和对数变换进行归一化,size.factors因子的值分别取决于基因表达总计数均值,或总计数的中位数,或常数10 000;对数转换消除了在scRNA-seq数据中经常出现的拖尾现象。

1.2 基因选择

scRNA-seq数据在大于10 000的维度中,实际上只有2 000-3 000个基因对单细胞聚类有作用[35]。因此,使用特征选择方法过滤对聚类作用不大的大部分基因,可以同时去除噪声和降低计算复杂度。特征选择方法是从高维特征中选择一组具有统计意义的原始特征的方法,降维后的特征仍然是原始特征,没有引入新的噪声,因此该方法越来越受到关注。但是,特征选择方法存在如何设计合适的选择策略,以发现具有实际意义特征子集的问题。

scRNA-seq数据的基因选择策略主要有3种类型。(1)基于高表达的基因选择(High Mean Gene,HMG)。该策略通过设置阈值,删除基因表达量均值低于阈值的基因来筛选出基因表达量均值高的基因。比如,Duò等[31]通过设置阈值为10%,筛选出基因表达量均值前10%的基因。(2)基于高差异表达的基因选择(High Variable Gene,HVG)。该策略通过量化每个基因在所有细胞中基因表达水平的差异度,筛选出高差异表达的基因。例如,Satija等[36]通过计算基因表达量的均值和离散度,量化基因表达水平差异,选择高差异表达的基因。(3)基于基因表达分布的基因选择(Drop-out based method)。该策略通过设计统计模型描述基因表达分布,并根据分布特性选择基因。表2列出了6个经典的基因选择方法,并从模型因子、测序平台、方法局限性和优势等方面进行对比分析。从表2可以看出,大多数基因选择方法中,均值是重要的模型因子,说明高表达基因在聚类中起着至关重要的作用。此外,少量或大量的基因识别数会影响聚类的性能。

表2 6种典型的基因选择方法及其特点对比

1.3 聚类方法

scRNA-seq数据缺乏细胞类型标签和类别数等先验知识,因此,无监督学习的聚类方法是识别细胞类型的重要方法。单细胞聚类方法包含传统的聚类方法和专门的聚类方法。

(1)传统的聚类方法

k均值聚类。Macosko等[21]通过使用k均值方法对基因进行聚类,识别具有相似表达的基因子集,在此基础上对细胞周期进行评分排序,根据评分实现单细胞聚类。Shin等[22]使用皮尔逊相关系数计算细胞间的相似性,采用最小生成树连接k个聚类中心得到细胞的发育轨迹,使用k均值聚类方法实现单细胞聚类。

层次聚类[37]。Llorens-Bobadilla等[32]使用欧氏距离计算细胞之间的距离,通过重采样估算类别数,采用层次聚类方法实现单细胞聚类。Darmanis等[38]使用皮尔逊相关系数计算细胞之间的相似性,使用层次聚类方法对单细胞聚类。

(2)基于映射的单细胞聚类方法

SHARP使用“分而治之”策略,将大规模数据分割成块[23]。使用稀疏随机投影(Random Projection,RP)算法,基于随机矩阵R将原始的D维数据映射为d维数据。随机矩阵R中的元素定义为D1/4乘以1(或0、-1),降维后的维数d定义为d=log2(N)/ε2[ε∈(0,1)]。运行k次RP算法得到k个d维矩阵,计算对应的k个相似性矩阵,在每个相似性矩阵上运行层次聚类,得到k个聚类结果。通过加权wMetaC方法集成这k个聚类结果,得到最终的聚类结果。

scDeepCluster融合零膨胀负二项(Zero Inflation Negative Binomial,ZINB)模型和自编码器,实现非线性函数映射并学习低维嵌入表示[24]。在自编码器中,引入随机高斯噪声以增强低维表示;在解码器中构建3个全连接层分别估计均值、离散度和缺失率对应的ZINB损失。使用Kullback-Leibler(KL)散度度量输入数据和重构数据间的分布差异,并构建新的损失函数。在输出的低维空间使用k-means进行聚类。

(3)基于图划分的单细胞聚类方法

Monocle 3使用统一流形逼近和投影(Uniform Manifold Approximation and Projection,UMAP)将scRNA-seq数据映射到低维空间,在此低维空间中使用Louvain社区检测算法实现图划分,对单细胞聚类,将相邻的细胞类合并为“超级类”[39]。最后,推断单个细胞在发育过程中的路径或轨迹,识别每个超级类的分支和合并位置。

(4)基于密度的单细胞聚类方法

SIMLR假设存在C个细胞类,那么细胞间的相似性矩阵应该具有C个近似的对角性块状结构,通过构造加权的高斯核函数,学习多个高斯核函数的权重[26]。定义细胞间的距离,从不同角度度量细胞间距离,构建对称的相似性矩阵S。同时,对相似性矩阵S、S上低秩约束的辅助低维矩阵L和权重ω进行优化学习。对学习得到的相似性矩阵直接使用亲和传播(Affinity Propagation,AP)算法聚类,或在降维后的低维空间使用k-means进行聚类。

Seurat集成了scRNA-seq数据和原位杂交空间转录组数据,通过识别高差异表达基因子集,学习标记基因的表达模型,去除标记基因表达的随机噪声[27]。通过将scRNA-seq数据估计的双峰表达模型与二值化的空间转录组数据对齐,建立基因表达统计模型,推断单细胞的空间位置。构建细胞共享近邻(Shared Nearest Neighbor,SNN)图,在共享近邻图中使用k-means对细胞聚类。

(5)基于集成的单细胞聚类方法

SC3分别用欧氏距离、Pearson相关系数、Spearman相关系数度量细胞间的距离[2]。利用主成分分析得到前d个主成分,在d个相似性矩阵中使用k-means聚类。根据节点对出现在同一类中的概率,将d个聚类结果集成到一个共识矩阵中,最后使用层次聚类对细胞聚类。

代表性的单细胞聚类方法及其特点对比分析见表3。从表3可以看出,基于映射的聚类方法将高维的scRNA-seq数据映射到低维空间,降低计算复杂度,可扩展性好,适用于大规模scRNA-seq数据,但是存在需要大内存的局限性。基于图划分、基于密度和基于集成的聚类方法,聚类准确性好,但存在计算复杂度高的局限性,适用于小规模scRNA-seq数据。

表3 6种典型的单细胞聚类方法及其特点对比

2 实验与结果分析

为了深入探讨不同方法对scRNA-seq数据分析的性能差异,本研究收集了14个scRNA-seq数据集,分别对6种基因选择方法在识别基因数、基因重叠度等方面进行对比分析。此外,在使用不同基因选择方法的基础上,分别对6种单细胞聚类方法的聚类准确性和稳定性等方面进行对比分析。

2.1 数据集

从基因表达综合数据库(GEO,https://www.ncbi.nlm.nih.gov/geo/)和欧洲生物信息学研究所网站(EMBL-EBI,https://www.ebi.ac.uk/)下载14个带有真实标签的金标准scRNA-seq数据集,包括5个全长测序数据集和9个双端测序数据集(表4)。这些数据集分别具有不同的测序协议、数据规模和稀疏度,其中5个数据集包含的细胞数大于3 000。这些数据集使用的Smart-seq2、10×genomics、Drop-seq等测序协议具有不同的特点:(1)与10×genomics相比,Smart-seq2具有更高的敏感度,可以检测到更多的基因,但其测序数据呈单峰分布,检测到的低表达基因少;(2)10×genomics数据呈双峰分布,可以检测到大量的0表达,这可能导致有更多的缺失(Dropout)事件,但它可以测序更多的细胞,更有效地检测罕见的细胞类型;(3)Drop-seq捕获效率较低,成本低,速度更快,不适合小样本测序。

表4 14个scRNA-seq数据集信息

2.2 不同基因选择方法检测基因的对比分析

为了观察不同基因选择方法选择基因的情况,使用6种基因选择方法分别在Wang_Lung、Adam数据集上检测基因。图2和图3分别是所检测基因的基因数、基因重叠度的upset图和韦恩图。在upset图中,横坐标是基因选择方法,包括检测到独有基因的基因选择方法,以及检测到共有基因的基因选择方法组合(这些基因选择方法由竖线相连),纵坐标是检测的基因数,其下方左侧是每个基因选择方法经过质量控制后留下的基因数。

图2 6种基因选择方法检测基因的upset图

从图2和图3可以看出,在Wang_Lung数据集中,除了Seurat方法以外,其他5种方法都检测到独有基因,Brennecke检测到120个独有基因,NBDropFS检测到116个独有基因,M3Drop检测到654个独有基因,Scran检测到79个独有基因,Monocle检测到813个独有基因,6种方法检测到124个共有基因。在Adam数据集中,Seurat过滤了绝大部分基因,仅保留了560个基因;Scran过滤了少部分基因,保留了11 515个基因。Scran检测到7 695个独有基因,Monocle检测到880个独有基因,Brennecke检测到14个独有基因,每种方法都与其他5种基因选择方法检测的基因有重叠度,6种方法检测到182个共有基因。Brennecke和Scran在两个数据集中检测的共有基因数分别是1 488和2 209,可以看出,这两种方法检测的基因重叠度比较大。Monocle与其他5种方法检测的基因重叠度相对比较小。

2.3 不同基因选择方法对单细胞聚类性能的影响

为了观察不同基因选择方法检测的基因对不同单细胞聚类方法的性能影响,分别对6种基因选择方法检测到的基因使用6种聚类方法进行单细胞聚类,采用调整的兰德指数(Adjusted Rand Index,ARI)[53]评价聚类性能。ARI度量了在预测类和真实类中都处在相同类的节点对的数量,其值的范围是-1到1。当ARI达到最大值1时,表示预测的类与真实类一致。

绘制ARI均值热图、ARI箱形图,以便更深入地观察和分析聚类性能的差异。融合6种基因选择方法和6种单细胞聚类方法,运行100次,其中5种单细胞聚类方法聚类结果的ARI均值见图4,聚类结果的ARI值箱形图见图5;此外,第6种单细胞聚类方法scDeepCluster取10个不同的随机种子,聚类结果的ARI箱形图见图6。实验中,Monocle分别采用了tSNE和UMAP降维,Seurat分别采用了PCA和ICA降维,其他方法没有进行降维。Monocle 3分别采用了densityPeak和Louvain对单细胞聚类,其他方法则分别采用自带的聚类方法。

图4 结合6种基因选择方法和5种单细胞聚类方法聚类结果的ARI均值

图6 结合6种基因选择方法,scDeepCluster聚类结果的ARI箱形图

从图4和图5可以看出,Seurat、SC3和Monocle 3结合不同的基因选择方法时,聚类性能ARI的稳定性更好;SHARP和SIMLR结合不同的基因选择方法时,聚类性能ARI的差异相对比较大。从图5可以看出,结合不同的基因选择方法,Seurat在所有数据集上的聚类稳定性和准确性最好,SC3次之,Monocle 3也比较好。此外,在Plasschaert、Wang_Kidney、Wang_Lung、Yong和Zeisel2015等5个数据集中,由于相似性计算的开销大,SIMIR方法没有实验结果。从图6可以看出,scDeepCluster结合不同的基因选择方法时,在所有数据集上也表现出比较好的稳定性,在大部分数据集上表现出很好的聚类性能。

3 讨论

为给研究者在选择合适的方法分析scRNA-seq数据时提供借鉴,本研究对比分析了scRNA-seq数据当前典型的质量控制、基因选择和聚类等方法。在对比分析质量控制时,发现通过设置不同的阈值,可以过滤低质量细胞和稀有基因,并且采用对数转换归一化可以消除数据拖尾现象。在对比分析基因选择时,通过比较6种典型的基因选择方法,发现均值是检测基因的重要模型因子,除Seurat以外的5种基因选择方法都使用了均值建模。此外,从实验结果可以看出,不同方法检测到一些相同的共有基因和少量的独有基因。6种基因选择方法在Adam和Wang_Lung数据集分别可以检测到182个和124个共有基因,Scran、Monocle、Brennecke、NBDropFS和M3Drop都检测到独有基因,Seurat则未检测到。检测到的共有基因包含了识别细胞类型的重要信息,检测到的独有基因反映了该方法建模条件下识别细胞类型的重要信息。在检测到的共有基因和独有基因的基础上,可以进一步分析它们在细胞发育过程轨迹推断中的作用。不同方法检测到的基因数有比较大的差异,Seurat检测到的基因数最少(<1 000),而Scran检测到的基因数最多(10 000左右)。在对比分析聚类时,结合6种不同基因选择方法,对6种单细胞聚类方法进行聚类性能比较,发现Seurat、SC3、Monocle 3和scDeepCluster的聚类稳定性较好,而SHARP和SIMLR的聚类稳定性则相对较差;Seurat在所有数据集上的聚类稳定性和准确性最好,scDeepCluster在大部分数据集上有很好的聚类准确性。因此,选择合适的scRNA-seq数据分析方法,需要综合考虑测序平台、数据规模,以及基因表达分布等因素。

随着第三代测序技术的迅速发展,产生了空间转录组(Spatial Transcriptome,ST)测序数据、单细胞基因组测序(single cell DNA sequencing,scDNA-seq)数据、单细胞甲基化测序(single cell methylation sequencing,sc-methyl-seq)数据等多种组学的测序数据,研究不同组学测序数据的对齐方法,有效融合多组学测序数据的重要信息,实现信息对齐和互补,有助于更准确地识别细胞类型。另外,当前的scRNA-seq数据具有长读长、大规模的新特点,长读长scRNA-seq数据存在更多的噪声,大规模scRNA-seq数据会导致更大的内存需求和计算时间开销问题,进一步研究基于数据分布的有效去噪方法、适合大规模数据的图神经网络降维方法,以提高数据质量并准确度量细胞间相似性,在细胞类型识别时加强生物可解释性,提升细胞类型识别和下游分析的性能等都是以后的重要工作。