01

概況

單細(xì)胞轉(zhuǎn)錄組測序（scRNA-seq）技術(shù)能夠?qū)?xì)胞群中的每一個細(xì)胞進(jìn)行大規(guī)模的全轉(zhuǎn)錄組分析。它的核心分析是將單細(xì)胞聚類，以揭示細(xì)胞亞型，并根據(jù)細(xì)胞之間的關(guān)系推斷細(xì)胞譜系。本文綜述了在過去幾年間發(fā)展起來的，用于單細(xì)胞轉(zhuǎn)錄組分析中聚類的機(jī)器學(xué)習(xí)和統(tǒng)計方法，重點(diǎn)介紹了如何將一些常見的聚類方法，如層次聚類、基于圖的聚類、混合模型、k-means、集成學(xué)習(xí)、神經(jīng)網(wǎng)絡(luò)和基于密度的聚類等加以調(diào)整及應(yīng)用，從而解決單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析中的獨(dú)特挑戰(zhàn)，例如低表達(dá)基因的缺失，轉(zhuǎn)錄本的不均勻覆蓋，以及由技術(shù)偏差和不相關(guān)的混雜生物變異所帶來的細(xì)胞標(biāo)記的失真。我們評價了標(biāo)準(zhǔn)化、dropouts推測以及降維等預(yù)處理步驟如何提高聚類效果。此外，我們還將介紹一些能夠?qū)r間序列樣本和多個細(xì)胞群進(jìn)行聚類并且檢測罕見細(xì)胞類型的新方法。最后，本文對部分開發(fā)用于單細(xì)胞轉(zhuǎn)錄組聚類分析的軟件進(jìn)行了實(shí)驗(yàn)和比較，以評估其性能和效率，為未來的數(shù)據(jù)分析提供一定的指導(dǎo)和方向。

02

介紹

細(xì)胞的轉(zhuǎn)錄組分析可以捕捉基因的表達(dá)活性，從而揭示細(xì)胞的身份和功能。在傳統(tǒng)的bulk-RNA測序中，轉(zhuǎn)錄組是通過從生物樣本中收集的大量細(xì)胞轉(zhuǎn)錄水平的平均值來測量的，這些平均后的表達(dá)值被用于基因共表達(dá)模塊的識別和樣本聚類。由于忽略了單個細(xì)胞的特性，這些傳統(tǒng)的方法無法在單細(xì)胞分辨率上研究重要的生物學(xué)問題，如細(xì)胞在早期發(fā)育過程中的不同功能角色、復(fù)雜組織中的不同細(xì)胞類型和細(xì)胞譜系關(guān)系。目前，scRNA-seq技術(shù)已廣泛用于量化單個細(xì)胞中的mRNA水平。在單細(xì)胞轉(zhuǎn)錄組的實(shí)驗(yàn)操作中，使用不同的捕獲方法（如FACS，F(xiàn)luidigm C1，microdroplet microfluidics）分離單細(xì)胞，然后對RNA進(jìn)行逆轉(zhuǎn)錄并擴(kuò)增測序。單細(xì)胞轉(zhuǎn)錄組的應(yīng)用已經(jīng)帶來了重要的生物學(xué)見解和發(fā)現(xiàn)，例如，對癌癥中腫瘤異質(zhì)性的理解。

細(xì)胞聚類是單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析中識別細(xì)胞亞群結(jié)構(gòu)的必要步驟，然而目前仍然存在一些挑戰(zhàn)。首先，由細(xì)胞的自身特征（如細(xì)胞所處周期階段、細(xì)胞大?。┖图夹g(shù)（捕獲方法、捕獲效率、PCR擴(kuò)增、測序深度等）引入的技術(shù)噪音和偏差。這些噪音和偏差將導(dǎo)致轉(zhuǎn)錄組的基因覆蓋極度不均勻，從而造成零覆蓋區(qū)域和dropouts的產(chǎn)生。另外，當(dāng)一個隊(duì)列的多個樣本同時進(jìn)行分析時，樣本間的技術(shù)偏差和變異將會主導(dǎo)細(xì)胞的聚類，導(dǎo)致細(xì)胞群體的形成更偏向于不同樣本來源而非細(xì)胞類型，即批次效應(yīng)。

在本文中，我們回顧了最近發(fā)展的用于提升單細(xì)胞轉(zhuǎn)錄組聚類效果或其相關(guān)的統(tǒng)計和機(jī)器學(xué)習(xí)方法。這些方法涉及：（1）用于基因表達(dá)值的標(biāo)準(zhǔn)化、dropouts推測、數(shù)據(jù)降維以及細(xì)胞特異Marker鑒定的數(shù)據(jù)預(yù)處理方法；（2）傳統(tǒng)的聚類算法，包括基于劃分的聚類、層次聚類、混合模型、基于圖的聚類、基于密度的聚類、神經(jīng)網(wǎng)絡(luò)、集成聚類和近鄰傳播聚類等；（3）在時間序列樣本和多個批次的細(xì)胞群中進(jìn)行聚類并檢測罕見細(xì)胞類型的新方法。我們還討論了單細(xì)胞轉(zhuǎn)錄組聚類分析中的幾個重要方面，包括細(xì)胞間相似性度量，特征值提取和單細(xì)胞聚類結(jié)果的評估。此外，我們對十多個軟件包進(jìn)行了比較，以評估它們在大規(guī)模單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集上的聚類性能和效率。最后，我們對聚類分析中存在的一些挑戰(zhàn)進(jìn)行了討論。

03

數(shù)據(jù)的預(yù)處理

在單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的聚類分析中，數(shù)據(jù)預(yù)處理對于減少技術(shù)變異和噪聲（如捕獲效率低、擴(kuò)增偏差、GC含量、總RNA含量和測序深度的差異等）以及建庫和測序過程中產(chǎn)生的dropouts至關(guān)重要。高維的基因表達(dá)矩陣通常需要經(jīng)過標(biāo)準(zhǔn)化及降維映射到低維空間中，一些計算方法還利用到統(tǒng)計學(xué)和數(shù)學(xué)方法來解決dropouts事件。

標(biāo)準(zhǔn)化

原始的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)通常從兩個層面進(jìn)行標(biāo)準(zhǔn)化：細(xì)胞的標(biāo)準(zhǔn)化和基因的標(biāo)準(zhǔn)化。細(xì)胞的標(biāo)準(zhǔn)化是為了消除擴(kuò)增偏差和其他細(xì)胞特異性的效應(yīng)，可以通過常用的reads計數(shù)標(biāo)準(zhǔn)化方法實(shí)現(xiàn)，如FPKM、RPKM、TPM等。基于UMI建庫的實(shí)驗(yàn)方案，理論上已經(jīng)避免了與擴(kuò)增或測序深度相關(guān)的誤差，因?yàn)楸幌嗤琔MI標(biāo)記的reads只會統(tǒng)計一次。然而，由于測序文庫通常是不飽和的，標(biāo)準(zhǔn)化對于該類型的數(shù)據(jù)也是有效的。細(xì)胞標(biāo)準(zhǔn)化的另一個方法是使用“spike-in”，它的基本思想是，由技術(shù)原因帶來的誤差對于內(nèi)外源基因的影響是相同的。另外，使用對數(shù)轉(zhuǎn)換進(jìn)行原始計數(shù)值的處理也非常常見。

基因標(biāo)準(zhǔn)化的目的是為了防止一些高表達(dá)基因主導(dǎo)了分析。常用的基因標(biāo)準(zhǔn)化方法如，在PCA中包含的z-score標(biāo)準(zhǔn)化。從過往的經(jīng)驗(yàn)中可以看到，基因的標(biāo)準(zhǔn)化可以提高算法的收斂和聚類效果。值得注意的是，數(shù)據(jù)的標(biāo)準(zhǔn)化處理將會使其失去原本基因表達(dá)的相對尺度，并且由于表達(dá)值的平移，造成表達(dá)矩陣變得不那么稀疏，這可能會影響到大規(guī)模數(shù)據(jù)集的聚類結(jié)果。

在SINCERA包中，對基因的標(biāo)準(zhǔn)化方法即是z-score，對細(xì)胞的標(biāo)準(zhǔn)化則是使用截尾均值（Trimmed mean）。一些工具會執(zhí)行更為特殊的標(biāo)準(zhǔn)化。例如，BISCUIT通過學(xué)習(xí)代表技術(shù)誤差的參數(shù)，在聚類過程中進(jìn)行迭代標(biāo)準(zhǔn)化；RaceID將每個細(xì)胞內(nèi)的總表達(dá)計數(shù)標(biāo)準(zhǔn)化到所有細(xì)胞表達(dá)計數(shù)的中位值。

此外，如果基因或者細(xì)胞顯現(xiàn)出極低的表達(dá)信號（基因表達(dá)值過低或者細(xì)胞表達(dá)基因過少），通常會將其移除，因?yàn)樗鼈兺碇摷傩盘?。在不同的研究中，為去除低表達(dá)基因和細(xì)胞建立了不同的閾值，這主要根據(jù)分析中囊括的細(xì)胞和基因的數(shù)量而有所不同。例如，scVDMC對PBMC樣本的處理中，表達(dá)值低于3的基因和總表達(dá)計數(shù)值小于200的細(xì)胞都將被去除。

雖然基因和細(xì)胞的標(biāo)準(zhǔn)化在目前大多數(shù)的單細(xì)胞數(shù)據(jù)分析流程中是常見的，但關(guān)于其對聚類結(jié)果的影響仍存在一些爭論。一項(xiàng)研究的分析表明，基于bulk的標(biāo)準(zhǔn)化方法在單細(xì)胞上的應(yīng)用可能會對其分析產(chǎn)生嚴(yán)重的不良后果，例如在聚類前進(jìn)行的高變基因的檢測。相同的，也有研究表明，通過中位數(shù)或者“spike-in”進(jìn)行標(biāo)準(zhǔn)化無法解決dropouts存在的問題，反而可能消除每種細(xì)胞類型特有的生物隨機(jī)性，這兩者都會導(dǎo)致潛在的細(xì)胞類型的不恰當(dāng)聚類或表征。

通過下面的例子，我們可以認(rèn)識到標(biāo)準(zhǔn)化的重要性。



Figure 1. 巨噬細(xì)胞群t-SNE圖 來自Zilionis等人數(shù)據(jù)集的巨噬細(xì)胞群t-SNE圖。(A)依據(jù)總計數(shù)值上色；(B)依據(jù)基因S100A9原始計數(shù)值上色；(C)依據(jù)標(biāo)準(zhǔn)化后的S100A9的表達(dá)值上色。

從圖1A,B很容易看出，S100A9的原始表達(dá)值與總計數(shù)高度相關(guān)，兩個圖的中心區(qū)域計數(shù)和表達(dá)量都較低，而外圍區(qū)域計數(shù)和表達(dá)量較高。我們能得出的唯一結(jié)論是，當(dāng)細(xì)胞中捕獲的轉(zhuǎn)錄本總量增加時，S100A9轉(zhuǎn)錄本的數(shù)量也會增加。這顯然沒什么意義。而在圖1C中，經(jīng)過標(biāo)準(zhǔn)化后的S100A9表達(dá)值與總計數(shù)之間似乎沒有相關(guān)性。我們可以說，S100A9表達(dá)的差異不依賴于測序深度等技術(shù)噪音，而應(yīng)該來自(主要)生物因素。

Dropout

單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)中一個重要的技術(shù)誤差被稱為“dropouts”。Dropout事件是指在反轉(zhuǎn)錄過程中由于缺失或轉(zhuǎn)錄本表達(dá)過低而導(dǎo)致基因未表達(dá)的錯誤定量。先前的研究也表明，簡單的數(shù)據(jù)標(biāo)準(zhǔn)化并不能解決該問題。因此，一些聚類算法中包含了特定的機(jī)制以矯正dropouts。例如，Seurat通過跨細(xì)胞的基因共表達(dá)模式，在聚類前進(jìn)行標(biāo)記基因的挑選。

另外，也可以通過計算配對相似性來估算dropouts。CIDR便是在聚類前進(jìn)行缺失值的填補(bǔ)。首先分析單細(xì)胞中可能出現(xiàn)的dropouts，識別每個細(xì)胞中的候選dropout基因，計算每個基因的dropout率；然后使用候選基因的dropout率來估算表達(dá)水平，即當(dāng)dropout事件以高概率被識別時，檢測算法會從其它細(xì)胞的表達(dá)譜中對該基因的表達(dá)值進(jìn)行填充；最后，利用矯正后的值計算細(xì)胞間的不相似度，進(jìn)行層次聚類。Seurat和SNN-Cliq是基于共享最近鄰SNN來度量細(xì)胞相似性。已經(jīng)證明，在稀疏的高維數(shù)據(jù)中，SNN考慮到周圍的近鄰數(shù)據(jù)點(diǎn)，更適合應(yīng)用于存在dropouts的聚類分析。

在一個更復(fù)雜的概率圖模型中，BISCUIT明確估計了每個細(xì)胞中的基因表達(dá)，以及通過數(shù)據(jù)分布和先驗(yàn)分布估算的代表技術(shù)和生物學(xué)變異的參數(shù)。其中，代表著未觀測到的基因真實(shí)表達(dá)水平的隨機(jī)變量被引入圖模型中并通過吉布斯抽樣來估算表達(dá)值。

降維

降維通常用于將高維基因表達(dá)數(shù)據(jù)投射到低維空間，使分析聚焦于低維空間中的相關(guān)信號，從而更好地實(shí)現(xiàn)數(shù)據(jù)的可視化、聚類分析等，幫助進(jìn)行生物學(xué)解釋。當(dāng)維數(shù)大于樣本數(shù)時，降維還有助于解決樣本不足的統(tǒng)計學(xué)問題。許多降維方法已經(jīng)應(yīng)用于單細(xì)胞轉(zhuǎn)錄組聚類算法，包括PCA、多維尺度變換(MDS)、t分布、隨機(jī)近鄰嵌入(t-SNE)、典型相關(guān)分析(CCA)、潛在狄利克雷分布(LDA)以及嵌入其他模型的降維等等。

PCA：將原本數(shù)據(jù)點(diǎn)映射到與協(xié)方差矩陣的最大特征值相關(guān)聯(lián)的特征向量（即主成分），以保留原始數(shù)據(jù)中的大部分方差。例如，pcaReduce在聚類前將表達(dá)矩陣映射到一個含有K-1個主成分的空間中；SC3使用PCA和拉普拉斯變換應(yīng)用于距離矩陣以獲得一致性矩陣并進(jìn)行層次聚類。此外，在聚類之后，PCA也被廣泛應(yīng)用于二維或三維的數(shù)據(jù)可視化。PCA是一種基于假設(shè)數(shù)據(jù)為高斯分布的線性投影方法，為了捕捉數(shù)據(jù)中的非線性結(jié)構(gòu)，可以使用核主成分分析與非線性核映射相結(jié)合。

MDS：也稱為主坐標(biāo)分析(PCoA)。MDS將數(shù)據(jù)點(diǎn)映射到低維空間，通過最小化所有配對數(shù)據(jù)點(diǎn)的原始空間中的距離與投影空間中的距離之間的差值，從而在低維嵌入保持原始高維空間中的數(shù)據(jù)點(diǎn)之間的距離。CIDR便是使用MDS來計算細(xì)胞的不相似矩陣。MDS的優(yōu)點(diǎn)是在低維空間中保持原始的成對距離，易于實(shí)現(xiàn)非線性特征嵌入。然而，MDS不能擴(kuò)展到大規(guī)模數(shù)據(jù)，因?yàn)楸仨氂嬎愠蓪嚯x來最小化目標(biāo)函數(shù)。

t-SNE：是一種將距離轉(zhuǎn)換為概率的方法。t-SNE構(gòu)造一個與原始空間及映射后的低維空間中數(shù)據(jù)點(diǎn)之間的相似性相關(guān)的概率分布，然后最小化兩個分布之間的Kullback-Leibler散度。t-SNE被廣泛應(yīng)用于單細(xì)胞數(shù)據(jù)分析中的數(shù)據(jù)可視化。

CCA：是一種基于互協(xié)方差矩陣的降維方法。給定兩個或多個數(shù)據(jù)集，該方法查找每個數(shù)據(jù)集的映射，以最大化數(shù)據(jù)集之間的相關(guān)性。在單細(xì)胞轉(zhuǎn)錄組的數(shù)據(jù)分析中，CCA通常用于不同來源樣本的整合，如Seurat（圖2）。

Figure 2. Seurat CCA數(shù)據(jù)整合示意圖

LDA：該方法最初是在自然語言處理中提出的。LDA假設(shè)一個文檔（document）是通過如下方法生成的：首先從具有狄利克雷先驗(yàn)的話題（topic）的多項(xiàng)分布中對話題進(jìn)行抽樣，然后對文檔中的單詞（word）進(jìn)行抽樣，這些單詞的多項(xiàng)式分布是基于每個話題的狄利克雷先驗(yàn)條件。然后，每個文檔都可以在包含k個話題的低維空間中表示。cellTree使用LDA學(xué)習(xí)“topics”作為潛在特征來表示細(xì)胞，其中“words”是受所選的潛在特征制約的基因表達(dá)水平。LDA的生成過程產(chǎn)生了一組可解釋的潛在特征。

相似度及核函數(shù)

在許多聚類方法的計算過程中，不是使用降維的方法，而是通過核函數(shù)或相似度函數(shù)來計算單個細(xì)胞之間的配對相似性進(jìn)行聚類。核函數(shù)策略將從N × M表達(dá)矩陣中計算獲得N × N相似矩陣，以期望通過核映射或相似函數(shù)在隱式特征映射空間中減少原始特征空間中的差異（如果使用有效的核函數(shù)）。SNN-cliq和Seurat使用SNN作為相似圖。cellTree在用LDA找到的話題直方圖上通過卡方找到細(xì)胞間的距離。DTWscore利用時間序列樣本為每個基因找到細(xì)胞對之間的動態(tài)時間規(guī)整（DTW）距離，以選擇高度可變的基因，其中DTW距離是基于兩個時間序列在最佳規(guī)整路徑上的比對計算的?；赥CC的聚類使用細(xì)胞間的Jensen-Shannon距離作為譜聚類或近鄰傳播聚類的輸入。SIMLR結(jié)合多個核來學(xué)習(xí)得到細(xì)胞相似矩陣，并使用秩約束和圖擴(kuò)散來解決dropouts問題。

大多數(shù)其他方法使用更標(biāo)準(zhǔn)的相似性函數(shù)或距離函數(shù)。BackSPIN，DendroSplit，ICGS和SINCERA在層次聚類策略中使用Pearson相關(guān)來尋找最佳分割點(diǎn)。GiniClust和RaceID也分別使用相關(guān)性矩陣進(jìn)行DBSCAN和k-means聚類。參考成分分析（RCA）計算單個細(xì)胞和參考細(xì)胞之間的表達(dá)譜之間的相關(guān)性，作為聚類的新特征，以最小化技術(shù)差異和批次效應(yīng)。SC3使用斯皮爾曼、皮爾森和歐氏距離來計算細(xì)胞間的配對相似性或距離以獲得一致性矩陣。

審核編輯：劉清