謝林娟，李荔瑄，張少?gòu)?qiáng)

天津師范大學(xué) 計(jì)算機(jī)與信息工程學(xué)院，天津 300387

在過(guò)去的幾年里，為了人類健康以及診斷、監(jiān)測(cè)和治療疾病，單細(xì)胞RNA測(cè)序（scRNA-seq）相關(guān)的研究得到了迅速的發(fā)展。包含數(shù)十萬(wàn)甚至上百萬(wàn)個(gè)單細(xì)胞測(cè)序數(shù)據(jù)的多個(gè)圖譜項(xiàng)目相繼完成[1-2]。這些項(xiàng)目主要?dú)w功于單細(xì)胞測(cè)序技術(shù)的發(fā)展。目前主流的測(cè)序平臺(tái)能夠一次完成數(shù)千到數(shù)萬(wàn)個(gè)細(xì)胞的測(cè)序[3]。但是，由于測(cè)序手段的限制以及基因表達(dá)高度復(fù)雜等原因，scRNA-seq數(shù)據(jù)普遍存在噪聲較大、維度較高、稀疏性較強(qiáng)等特點(diǎn)[4]。細(xì)胞分型是單細(xì)胞數(shù)據(jù)分析的重要流程之一[5]。現(xiàn)在主流的細(xì)胞分型方法先運(yùn)用主成分分析（principal component analysis，PCA）、t分布隨機(jī)近鄰嵌入（t-dis‐tributed stochastic neighbor embedding，t-SNE）或均勻流形近似和投影（uniform manifold approximation and projection，UMAP）[6]等經(jīng)典數(shù)據(jù)降維方法進(jìn)行數(shù)據(jù)降維或者進(jìn)行特征選擇，然后采用DBSCAN、譜聚類、層次聚類、K-means等經(jīng)典聚類方法進(jìn)行分型。例如SC3[7]運(yùn)用PCA和t-SNE降 維及K-means進(jìn)行 一 致性聚類；CI‐DR[8]運(yùn)用PCA降維和層次聚類；SCENA先特征選擇再譜聚類[9]。在幾千細(xì)胞的小規(guī)模測(cè)序數(shù)據(jù)，這些算法均具有良好的性能。

目前針對(duì)大規(guī)模單細(xì)胞RNA測(cè)序數(shù)據(jù)的聚類基本采用構(gòu)建簡(jiǎn)單的細(xì)胞關(guān)系網(wǎng)絡(luò)，并對(duì)該網(wǎng)絡(luò)運(yùn)用社區(qū)檢測(cè)的算法（例如魯汶算法[10]）來(lái)進(jìn)行。而譜聚類和層次聚類算法受限于計(jì)算復(fù)雜性，很難應(yīng)用于大數(shù)據(jù)集。PhenoGraph[11]需要是先構(gòu)造k最近鄰（k-nearest neigh‐bors，KNN）圖然后對(duì)每一對(duì)細(xì)胞的共享最近鄰居（shared nearest neighbors，SNN）重新加權(quán)再用魯汶社區(qū)檢測(cè)算法聚類。SCANPY[12]基于UMAP降維后構(gòu)造KNN圖再調(diào)用PhenoGraph聚類。DropClust[13]使用了對(duì)數(shù)時(shí)間的局部敏感哈希（locality sensitive hashing，LSH）算法來(lái)確定每個(gè)細(xì)胞的近似KNN，然后用魯汶算法聚類。PARC[14]用對(duì)數(shù)時(shí)間復(fù)雜性的“分層導(dǎo)航小世界”（hierarchical navigable small world，HNSW）[15]搜索每個(gè)細(xì)胞的近似KNN，然后用萊頓模塊優(yōu)化（Leiden modularity optimization）[16]的社區(qū)檢測(cè)算法進(jìn)行聚類。這些大規(guī)模聚類工具大都通過(guò)近似構(gòu)造KNN圖來(lái)提高計(jì)算速度。雖然基于社區(qū)檢測(cè)的萊頓算法和魯汶算法能夠保證計(jì)算速度，但是其準(zhǔn)確度低于其他聚類算法[17-18]。scAIDE[19]是一個(gè)先用自動(dòng)編碼器插補(bǔ)和降維，然后用隨機(jī)投影哈希（random projection hashing，RPH）初始化聚類質(zhì)心（centroid）來(lái)加速k-means的方法，該工具可以快速處理大規(guī)模數(shù)據(jù)，但需要事先確定類簇?cái)?shù)目。

因此如果能夠構(gòu)造更符合細(xì)胞關(guān)系的網(wǎng)絡(luò)圖并能夠在聚類過(guò)程中自動(dòng)確定類簇?cái)?shù)目，將能夠提升單細(xì)胞數(shù)據(jù)的聚類準(zhǔn)確度。此外，在聚類之前耗費(fèi)時(shí)間最長(zhǎng)的是計(jì)算細(xì)胞間的相似度。為此，本文結(jié)合openMP和CUDA來(lái)進(jìn)行CPU+GPU的異構(gòu)并行計(jì)算，并采用KNN和細(xì)胞相似度閾值相結(jié)合的方法構(gòu)建相似網(wǎng)絡(luò)，最后通過(guò)改進(jìn)馬爾可夫聚類算法（Markov clustering，MCL）[20]來(lái)進(jìn)行快速聚類。該算法命名為SCMC（single-cell Markov clustering），其C++代碼在Github開(kāi)源共享，其代碼及算法說(shuō)明文檔詳見(jiàn)https：//github.com/shaoqiang‐zhang/SCMC/。

1 方法與數(shù)據(jù)

1.1 方法概述

給定由細(xì)胞集合C組成的單細(xì)胞表達(dá)矩陣X，其中m為基因數(shù)，n為細(xì)胞數(shù)。SCMC首先對(duì)輸入表達(dá)矩陣X的每個(gè)細(xì)胞的元素進(jìn)行標(biāo)準(zhǔn)化處理；然后計(jì)算所有基因的變異系數(shù)，并根據(jù)變異系數(shù)閾值進(jìn)行特征基因的選擇；再運(yùn)用并行計(jì)算技術(shù)計(jì)算細(xì)胞之間的相似度構(gòu)建細(xì)胞相似網(wǎng)絡(luò)；最后通過(guò)采用修改的馬爾科夫聚類算法對(duì)細(xì)胞相似度網(wǎng)絡(luò)進(jìn)行快速聚類，得到細(xì)胞類型。算法流程見(jiàn)圖1，具體方法步驟見(jiàn)1.2節(jié)至1.5節(jié)描述。

圖1 SCMC算法流程圖Fig.1 Flow chart of SCMC

1.2 數(shù)據(jù)預(yù)處理和特征基因選擇

數(shù)據(jù)預(yù)處理分為以下兩個(gè)步驟：基因過(guò)濾和細(xì)胞數(shù)據(jù)對(duì)數(shù)標(biāo)準(zhǔn)化（log-normalization）。即先從X中移除在所有細(xì)胞表達(dá)值均為0的基因（行），然后對(duì)矩陣的每個(gè)元素除以該元素所在細(xì)胞列的元素總和再乘以一個(gè)比例因子（默認(rèn)為10 000），最后加上偽數(shù)1取對(duì)數(shù)，見(jiàn)公式（1）：

對(duì)對(duì)數(shù)標(biāo)準(zhǔn)化后的矩陣(x'ij)，計(jì)算每個(gè)基因的變異系數(shù)（coefficient of variation）。變異系數(shù)為基因所在行元素的標(biāo)準(zhǔn)差除以均值。根據(jù)所有基因變異系數(shù)的長(zhǎng)尾分布（heavy-tailed）曲線設(shè)置變異系數(shù)變化趨緩的長(zhǎng)尾閾值α，根據(jù)閾值α選取變異系數(shù)最大的前T個(gè)基因作為特征基因。注意不同的數(shù)據(jù)集會(huì)選擇不同的閾值，從而選擇不同數(shù)量的特征基因。

1.3 計(jì)算細(xì)胞相似度構(gòu)建擴(kuò)充的K近鄰網(wǎng)絡(luò)

構(gòu)建K近鄰（K-nearest neighbors，KNN）網(wǎng)絡(luò)是目前多數(shù)大規(guī)模細(xì)胞聚類的通用方法。在構(gòu)建KNN網(wǎng)絡(luò)之前，需要計(jì)算降維后的兩兩細(xì)胞對(duì)應(yīng)特征表達(dá)向量的歐氏距離。設(shè)細(xì)胞i和j對(duì)應(yīng)的特征基因的表達(dá)向量分別為x和y，按照公式（2）計(jì)算兩細(xì)胞之間的余弦相似度。

根據(jù)公式（3）的推導(dǎo)，計(jì)算標(biāo)準(zhǔn)化后兩向量的歐氏距離D等價(jià)于計(jì)算向量的余弦距離[21]。余弦相似度相當(dāng)于對(duì)細(xì)胞的特征表達(dá)向量做了標(biāo)準(zhǔn)化處理再求歐氏距離D，因此無(wú)需對(duì)降維后的細(xì)胞表達(dá)向量再標(biāo)準(zhǔn)化。

因?yàn)橛?jì)算一對(duì)細(xì)胞的余弦相似度不需要與別的細(xì)胞交換信息，因此將n個(gè)細(xì)胞的所有可能的C2n種細(xì)胞對(duì)（cell pairs）組合分配到GPU的所有k個(gè)內(nèi)核上進(jìn)行CUDA并行計(jì)算。這樣理論上比CPU單線程的運(yùn)算提高k倍。

對(duì)每個(gè)細(xì)胞x,計(jì)算x與其他所有細(xì)胞間的余弦相似度，并保留相似度最高的前條邊集合Nx和相似度超過(guò)閾值β的邊集Bx,得到與頂點(diǎn)x關(guān)聯(lián)的邊集Nx?Bx。由此構(gòu)造擴(kuò)充的K近鄰網(wǎng)絡(luò)（expanded KNN），其中頂點(diǎn)集合為細(xì)胞集合C，邊集為

1.4 馬爾科夫聚類算法的并行實(shí)現(xiàn)

運(yùn)用馬爾科夫聚類（Markov clustering，MCL）算法對(duì)擴(kuò)充的KNN網(wǎng)絡(luò)進(jìn)行聚類。MCL算法是van Dongen[20]開(kāi)發(fā)的基于圖的聚類的算法，主要思想是通過(guò)概率狀態(tài)轉(zhuǎn)移矩陣進(jìn)行隨機(jī)游走（random walk）。該算法不需預(yù)先設(shè)定聚類數(shù)目，通過(guò)反復(fù)修正矩陣以實(shí)現(xiàn)隨機(jī)流模擬。

首先根據(jù)細(xì)胞相似網(wǎng)絡(luò)構(gòu)造鄰接矩陣，然后對(duì)鄰接矩陣進(jìn)行標(biāo)準(zhǔn)化，即每個(gè)矩陣元素除以所在列的所有元素之和得到概率轉(zhuǎn)移矩陣，如公式（4）：

構(gòu)造鄰接矩陣每列元素的概率用不到其他列的信息，因此針對(duì)每列進(jìn)行并行化設(shè)計(jì)（即矩陣每列指派不同的線程處理）。MCL算法是“擴(kuò)展”（Expansion）操作和“膨脹”（Inflation）操作交替重復(fù)執(zhí)行的過(guò)程。首先將概率矩陣進(jìn)行擴(kuò)展操作，對(duì)矩陣N*進(jìn)行p次冪方得到矩陣M如公式（5）：

隨即進(jìn)行膨脹操作，對(duì)矩陣M內(nèi)元素進(jìn)行r次冪方得到矩陣G，如公式（6）：

判斷矩陣G與擴(kuò)展操作前的矩陣N*是否非常接近，即矩陣G與N*所有相同位置元素之差小于規(guī)定的誤差ε（例如ε=10-4），則收斂，輸出聚類結(jié)果；否則令N*=G繼續(xù)迭代執(zhí)行擴(kuò)展操作和膨脹操作直至收斂。

在實(shí)際聚類中，設(shè)置p=2,r=1.5。其中擴(kuò)展操作（公式（5））實(shí)際就是矩陣的乘法運(yùn)算，矩陣被按行和列分塊進(jìn)行并行計(jì)算；而膨脹操作（公式（6））與（公式（4））類似的并行設(shè)計(jì)，即對(duì)G的每列進(jìn)行并行化處理。本文調(diào)用openMP并行庫(kù)（https：//www.openmp.org/）來(lái)并行化MCL的C++代碼。

1.5 馬爾科夫聚類的改進(jìn)

原始的馬爾科夫聚類結(jié)果常常會(huì)產(chǎn)生規(guī)模很小的類簇，甚至只包含一個(gè)頂點(diǎn)的孤立點(diǎn)簇（統(tǒng)稱“小類簇”）。對(duì)于較大規(guī)模的單細(xì)胞數(shù)據(jù)，這些小類簇應(yīng)該不是新的細(xì)胞類型。為此，為了提高聚類的準(zhǔn)確度，SCMC將小類簇重新分配到其他類簇（統(tǒng)稱“大類簇”）中。對(duì)每個(gè)小類簇S，在擴(kuò)充的KNN網(wǎng)絡(luò)中標(biāo)記與該小類簇每個(gè)頂點(diǎn)s∈S關(guān)聯(lián)的相似度最高的頂點(diǎn)所屬的大類簇編號(hào)l(s)；從而得到小類簇S關(guān)聯(lián)度最高的大類簇編號(hào)集合L(S)={l (s),s∈S}。取L(S)中出現(xiàn)頻次最高的大類簇作為S的最終目標(biāo)類簇，將S分配給該類簇。如果n中存在多個(gè)頻次最高的大類簇，則選擇與S中所有頂點(diǎn)相似度最高的頂點(diǎn)所屬類簇。

1.6 數(shù)據(jù)集和評(píng)估指標(biāo)

為了測(cè)試算法的性能，本文選取7個(gè)較大規(guī)模具有代表性的單細(xì)胞測(cè)序數(shù)據(jù)集（依次編號(hào)Data1~Data7）。這些數(shù)據(jù)集分別從基因表達(dá)綜合數(shù)據(jù)庫(kù)（gene expres‐sion omnibus，GEO）和10X Genomics官網(wǎng)（https：//www.10xgenomics.com/resources/datasets）免費(fèi)獲取。Data1為3 005個(gè)小鼠大腦皮層（cerebral cortex）細(xì)胞scRNA-seq數(shù)據(jù)集（GEO編號(hào)：GSE60361）[22]；Data2為10X Genom‐ics的4 340個(gè)健康人體外周血單核細(xì)胞（peripheral blood mononuclear cells，PBMCs）數(shù)據(jù)集；Data3、Da‐ta4、Data5分別為三個(gè)混合人體HEK293T細(xì)胞和小鼠NIH3T3細(xì)胞的數(shù)據(jù)集（樣本標(biāo)識(shí)分別為10k_hgmm_3p、10k_hgmm_3p_nextgem_Chromium_X 和20k_hgmm_5pv2_HT_nextgem_Chromium_X）；Data6為24 822個(gè) 小鼠前額葉皮層（prefrontal cortex）細(xì)胞數(shù)據(jù)集（GEO編號(hào)：GSE124952）[23]；Data7為68 579個(gè)健康人體外周血單核細(xì)胞（peripheral blood mononuclear cells，PBMCs）數(shù)據(jù)集。實(shí)驗(yàn)所用數(shù)據(jù)集的具體信息見(jiàn)表1。

表1 實(shí)驗(yàn)數(shù)據(jù)集Table 1 Experimental data set

調(diào)整蘭德系數(shù)（adjusted Rand index，ARI）[24]被廣泛用于評(píng)估單細(xì)胞數(shù)據(jù)聚類效果[25]。給定n個(gè)目標(biāo)（細(xì)胞）的聚類結(jié)果X=( )X1,X2,…,Xr和真實(shí)劃分結(jié)果調(diào)整蘭德系數(shù)的定義如公式（7）：

其中，nij代表Xi和Yj中均有的目標(biāo)數(shù)目，ai和bj分別表示在Xi和Yj中的細(xì)胞數(shù)目。

1.7 方法選擇和參數(shù)設(shè)置

在大規(guī)模單細(xì)胞聚類中應(yīng)用最廣泛的是SEURAT4[26]，該工具均直接調(diào)用PhenoGraph中的魯汶算法進(jìn)行聚類。SCENA是最新的單細(xì)胞聚類算法，該方法基于多特征選擇的譜聚類算法并聲明在大多數(shù)實(shí)驗(yàn)數(shù)據(jù)上優(yōu)于SEURAT4和SC3等主流工具[9]。本文重點(diǎn)將SCMC與SCENA和SEURAT4，以及另一基于魯汶算法的工具SCANPY進(jìn)行比較。所有比較工具的參數(shù)均采用默認(rèn)參數(shù)，其中KNN的參數(shù)K均設(shè)置為相同值。單細(xì)胞數(shù)據(jù)的預(yù)處理（包括對(duì)數(shù)標(biāo)準(zhǔn)化和刪除低表達(dá)基因等）均按照各個(gè)工具默認(rèn)參數(shù)進(jìn)行。所有程序均在Linux工作站（CPU為Intel Xeon E5-2620/2.10 GHz/8×2 cores，GPU為NVidia GTX1080Ti）上運(yùn)行。

2 方法與數(shù)據(jù)

2.1 擴(kuò)充KNN網(wǎng)絡(luò)提高聚類效果

以前的方法大都根據(jù)細(xì)胞相似性構(gòu)建KNN網(wǎng)絡(luò)或者近似KNN網(wǎng)絡(luò)，然后對(duì)網(wǎng)絡(luò)圖進(jìn)行聚類。直接用KNN網(wǎng)絡(luò)進(jìn)行聚類效果并不理想，單純KNN網(wǎng)絡(luò)沒(méi)有充分考慮細(xì)胞之間的相似度。一對(duì)高度相似表達(dá)的細(xì)胞因?yàn)閿?shù)目K的限制，可能在KNN圖中沒(méi)有直接連接，因此通過(guò)保留高相似度（相似度>0.95）的連接來(lái)擴(kuò)充KNN網(wǎng)絡(luò)。如表2所示，分別在7個(gè)單細(xì)胞數(shù)據(jù)集上構(gòu)造單純的KNN網(wǎng)絡(luò)與擴(kuò)充的KNN網(wǎng)絡(luò)并進(jìn)行后續(xù)聚類步驟，發(fā)現(xiàn)對(duì)大部分?jǐn)?shù)據(jù)集，在擴(kuò)充KNN網(wǎng)絡(luò)的聚類比在單純KNN網(wǎng)絡(luò)上聚類效果顯著增強(qiáng)。只有Data2數(shù)據(jù)集效果沒(méi)有顯著變化，這可能該數(shù)據(jù)集本身容易聚類(ARI=0.974)。而在Data6數(shù)據(jù)集上擴(kuò)充的KNN網(wǎng)絡(luò)比單純的KNN提升效果非常顯著。

表2 KNN網(wǎng)絡(luò)與擴(kuò)充KNN網(wǎng)絡(luò)實(shí)驗(yàn)比較Table 2 Experimental comparison between KNN networks and extended KNN networks

SCMC算法的主要參數(shù)指標(biāo)是KNN中的最近鄰數(shù)K。選取三個(gè)不同規(guī)模的數(shù)據(jù)（Data2、Data4和Data6）來(lái)驗(yàn)證不同K值對(duì)最終結(jié)果的影響。分別選取K=10、20、30、40，然后進(jìn)行聚類并計(jì)算ARI。如圖2所示，在這三個(gè)數(shù)據(jù)中SCMC算法的主要參數(shù)指標(biāo)是KNN中的最近鄰數(shù)K。如圖2所示，在這三個(gè)數(shù)據(jù)集中，K≥20表現(xiàn)比較平穩(wěn)，其中最大的數(shù)據(jù)集Data6，隨著K增大ARI也相應(yīng)也有所增加。但K越大，KNN網(wǎng)絡(luò)的密度就會(huì)增大，算法復(fù)雜性就會(huì)增高。為此在所有較大規(guī)模數(shù)據(jù)集，為了減少算法復(fù)雜度，建議K設(shè)置為20。

圖2 不同K值對(duì)ARI的影響Fig.2 Effect of different K values on ARI

聚類算法前的特征基因選擇也會(huì)影響最終的聚類結(jié)果，SCMC算法通過(guò)計(jì)算基因的變異系數(shù)（coefficient of variation）來(lái)選擇高可變的基因作為特征基因。圖3是7個(gè)數(shù)據(jù)集基因變異系數(shù)的分布曲線，這些曲線的形態(tài)基本一致；變異系數(shù)越高的基因越可能是特征基因，根據(jù)帕累托法則（Pareto principle），重尾（heavy-tailed）部分大概率是特征顯著的基因。由于帕累托分布（Pareto distribution）屬于重尾連續(xù)概率分布，與圖3的曲線形態(tài)基本一致，為此，SCMC用帕累托分布來(lái)擬合基因變異系數(shù)的分布曲線（可以選用fitdistrplus的R包實(shí)現(xiàn)），選取面積占比5%的重尾部分的高可變基因作為特征基因，如圖3所示。在算法實(shí)際運(yùn)行中，也可以觀測(cè)曲線選擇變化開(kāi)始趨緩的值作為特征基因選擇的閾值（例如：曲線a點(diǎn)的值v(a)小于下一點(diǎn)a+1值v( a+1)的2倍，v(a)<2v( a+1)，則取值a為閾值來(lái)選擇特征基因）。

圖3 各數(shù)據(jù)集基因變異系數(shù)分布曲線Fig.3 Distributi on curveDsa tao 5f coefficien of gene variation for 7 datasets

2.2 改進(jìn)MCL以提高聚類效果

MCL算法是經(jīng)典的圖聚類算法，該算法已經(jīng)被成功應(yīng)用到模體聚類等研究領(lǐng)域。當(dāng)數(shù)據(jù)規(guī)模很大的時(shí)候，該算法復(fù)雜性較高( O(n2))以及容易產(chǎn)生大量的離群孤立或者小規(guī)模類簇，因此在單細(xì)胞聚類中，首先對(duì)MCL算法進(jìn)行openMP并行設(shè)計(jì)來(lái)加快算法，另外對(duì)MCL算法產(chǎn)生的孤立點(diǎn)和小類簇進(jìn)行再分配，以提升聚類效果。

實(shí)驗(yàn)中將類簇規(guī)模小于3的進(jìn)行再分配。如表3所示，分別比較了改進(jìn)前和改進(jìn)后MCL聚類的ARI結(jié)果。除了Data2、Data3和Data4的ARI沒(méi)有明顯變化，算法在其他數(shù)據(jù)上ARI指標(biāo)均有所增加。說(shuō)明對(duì)MCL的改進(jìn)能夠提升聚類的效果。另外，改進(jìn)后的MCL聚類得到的類簇?cái)?shù)與數(shù)據(jù)集提供的細(xì)胞類別數(shù)比較接近。對(duì)于Data7，改進(jìn)后MCL的聚類結(jié)果的類簇?cái)?shù)大幅度減少，由此可以看出算法對(duì)于大規(guī)模數(shù)據(jù)類簇?cái)?shù)減少明顯，說(shuō)明MCL本身產(chǎn)生了較多的小類簇。

2.3 并行異構(gòu)計(jì)算加速聚類算法

由于SCMC沒(méi)有對(duì)細(xì)胞的特征基因進(jìn)行降維，而直接計(jì)算細(xì)胞之間的相似度，因此計(jì)算細(xì)胞之間相似度并構(gòu)建擴(kuò)充的KNN圖會(huì)占用大量的計(jì)算時(shí)間。為此，SCMC通過(guò)CUDA并行化來(lái)計(jì)算細(xì)胞相似度，通過(guò)openMP并行化MCL聚類來(lái)加快算法。從Data7中隨機(jī)取樣不同固定數(shù)量的細(xì)胞（從4 000到24 000個(gè)細(xì)胞）樣本來(lái)分別評(píng)估SCMC在CPU+GPU異構(gòu)并行的運(yùn)行時(shí)間。如果單純用30個(gè)線程的CPU，4 000個(gè)細(xì)胞就耗時(shí)十幾個(gè)小時(shí)。但是如圖4所示，如果采用CPU（30線程）+GPU異構(gòu)計(jì)算，對(duì)于24 000個(gè)細(xì)胞的大規(guī)模數(shù)據(jù)集，SCMC的運(yùn)行時(shí)間僅用時(shí)0.87 h。對(duì)于4 000個(gè)細(xì)胞的數(shù)據(jù)集，SCMC的運(yùn)行時(shí)間則僅需要8 min左右。

0.05表3 MCL改進(jìn)前后聚類結(jié)果比較Table 3 Comparison of clustering results before and after MCL improvement

圖4 SCMC算法不同細(xì)胞數(shù)量的運(yùn)行時(shí)間Fig.4 Running time of SCMC on different numbers of cells

2.4 SCMC與其他聚類方法的性能比較

將SCMC與2種最流行的算法SEURAT4和SCANPY，以及最新聚類算法SCENA進(jìn)行性能比較。分別用這4個(gè)工具對(duì)表1的7個(gè)數(shù)據(jù)集進(jìn)行聚類，并計(jì)算聚類的ARI指標(biāo)。

如圖5所示，SCMC的 ARI 指標(biāo)在這些工具中始終表現(xiàn)最好，特別是SCMC在所有測(cè)試數(shù)據(jù)上的 ARI 均大于0.7，算法具有較好的魯棒性。由于SCMC在大部分實(shí)驗(yàn)數(shù)據(jù)上要比其他三個(gè)算法更接近真實(shí)類簇（細(xì)胞類型）數(shù)目，因此反映在ARI指標(biāo)上要優(yōu)于其他算法。因?yàn)镾CANPY和SEURAT4均使用社區(qū)檢測(cè)算法，在默認(rèn)聚類粒度參數(shù)下，大都屬于粗粒度聚類，因此得到類簇?cái)?shù)與實(shí)際細(xì)胞類型數(shù)相差較大。SCENA采用譜聚類，并用近鄰傳播聚類估計(jì)類簇?cái)?shù)。在Data6數(shù)據(jù)集上SCENA預(yù)測(cè)的聚類數(shù)是14，更接近真實(shí)聚類數(shù)目，由此SCMC與SCENA的ARI指標(biāo)比較接近。但在其他測(cè)試數(shù)據(jù)上性能較明顯低于SCMC，說(shuō)明SCENA的魯棒性較差SCANPY聚類。

作為無(wú)監(jiān)督聚類方法，SCMC在聚類中無(wú)須在聚類前估計(jì)細(xì)胞類簇的數(shù)目。雖然SCMC輸出的類簇?cái)?shù)與數(shù)據(jù)本身提供的類簇?cái)?shù)有所出入，但是較高的ARI指標(biāo)說(shuō)明SCMC可能是對(duì)大類細(xì)胞進(jìn)行了亞群的細(xì)分。例如，SCMC 預(yù)測(cè) Data7 數(shù)據(jù)（pbmc68K）有 13 個(gè)細(xì)胞類型，接近其真實(shí)的 11 個(gè)細(xì)胞亞型。再如圖 6 所示，對(duì)Data2數(shù)據(jù)的聚類結(jié)果通過(guò)UMAP降維后取前兩維作圖發(fā)現(xiàn)，SCMC聚類標(biāo)簽與原始標(biāo)簽吻合度最高，并且發(fā)現(xiàn)右下角存在稀有細(xì)胞類型，但在原始標(biāo)簽中該部分被分到了距離較遠(yuǎn)的細(xì)胞類型中。從圖6（c）、（d）、（e）可以看出另外三個(gè)聚類算法可視化的聚類標(biāo)簽比較雜亂。此外，SCMC算法另一優(yōu)點(diǎn)是占用內(nèi)存低。在構(gòu)建擴(kuò)充KNN圖的時(shí)候，通過(guò)CUDA將每一對(duì)細(xì)胞分別輸入到GPU計(jì)算單元，計(jì)算每對(duì)細(xì)胞的相似度后自動(dòng)釋放內(nèi)存，因此SCMC算法在計(jì)算細(xì)胞相似性時(shí)候只占用小量GPU存儲(chǔ)，而其他基于矩陣運(yùn)算的算法（例如SCENA等）基本要一次讀入所有細(xì)胞數(shù)據(jù)，占用大量的內(nèi)存。

圖6 UMAP降維后Data2數(shù)據(jù)集的可視化結(jié)果Fig.6 Clustering visualization effects of Data2 after UMAP dimensionality reduction

3 結(jié)論與展望

SCMC算法將擴(kuò)充的KNN網(wǎng)絡(luò)與改進(jìn)的MCL聚類算法進(jìn)行結(jié)合，并且基于CPU+GPU進(jìn)行并行化程序設(shè)計(jì)，因此適用于規(guī)模較大的單細(xì)胞RNA測(cè)序數(shù)據(jù)的細(xì)胞分類。通過(guò)在7個(gè)10X單細(xì)胞測(cè)序數(shù)據(jù)上實(shí)驗(yàn)比較，算法比其他幾個(gè)最流行的算法聚類效果要好，因此該算法比較適用于10X單細(xì)胞數(shù)據(jù)分析。