01

概況

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

02

介紹

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

細胞聚類是單細胞轉錄組數據分析中識別細胞亞群結構的必要步驟，然而目前仍然存在一些挑戰。首先，由細胞的自身特征（如細胞所處周期階段、細胞大小）和技術（捕獲方法、捕獲效率、PCR擴增、測序深度等）引入的技術噪音和偏差。這些噪音和偏差將導致轉錄組的基因覆蓋極度不均勻，從而造成零覆蓋區域和dropouts的產生。另外，當一個隊列的多個樣本同時進行分析時，樣本間的技術偏差和變異將會主導細胞的聚類，導致細胞群體的形成更偏向于不同樣本來源而非細胞類型，即批次效應。

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

03

數據的預處理

在單細胞轉錄組數據的聚類分析中，數據預處理對于減少技術變異和噪聲（如捕獲效率低、擴增偏差、GC含量、總RNA含量和測序深度的差異等）以及建庫和測序過程中產生的dropouts至關重要。高維的基因表達矩陣通常需要經過標準化及降維映射到低維空間中，一些計算方法還利用到統計學和數學方法來解決dropouts事件。

標準化

原始的單細胞轉錄組數據通常從兩個層面進行標準化：細胞的標準化和基因的標準化。細胞的標準化是為了消除擴增偏差和其他細胞特異性的效應，可以通過常用的reads計數標準化方法實現，如FPKM、RPKM、TPM等。基于UMI建庫的實驗方案，理論上已經避免了與擴增或測序深度相關的誤差，因為被相同UMI標記的reads只會統計一次。然而，由于測序文庫通常是不飽和的，標準化對于該類型的數據也是有效的。細胞標準化的另一個方法是使用“spike-in”，它的基本思想是，由技術原因帶來的誤差對于內外源基因的影響是相同的。另外，使用對數轉換進行原始計數值的處理也非常常見。

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

在SINCERA包中，對基因的標準化方法即是z-score，對細胞的標準化則是使用截尾均值（Trimmed mean）。一些工具會執行更為特殊的標準化。例如，BISCUIT通過學習代表技術誤差的參數，在聚類過程中進行迭代標準化；RaceID將每個細胞內的總表達計數標準化到所有細胞表達計數的中位值。

此外，如果基因或者細胞顯現出極低的表達信號（基因表達值過低或者細胞表達基因過少），通常會將其移除，因為它們往往代表著虛假信號。在不同的研究中，為去除低表達基因和細胞建立了不同的閾值，這主要根據分析中囊括的細胞和基因的數量而有所不同。例如，scVDMC對PBMC樣本的處理中，表達值低于3的基因和總表達計數值小于200的細胞都將被去除。

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

通過下面的例子，我們可以認識到標準化的重要性。



Figure 1. 巨噬細胞群t-SNE圖 來自Zilionis等人數據集的巨噬細胞群t-SNE圖。(A)依據總計數值上色；(B)依據基因S100A9原始計數值上色；(C)依據標準化后的S100A9的表達值上色。

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

Dropout

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

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

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

降維

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

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

MDS：也稱為主坐標分析(PCoA)。MDS將數據點映射到低維空間，通過最小化所有配對數據點的原始空間中的距離與投影空間中的距離之間的差值，從而在低維嵌入保持原始高維空間中的數據點之間的距離。CIDR便是使用MDS來計算細胞的不相似矩陣。MDS的優點是在低維空間中保持原始的成對距離，易于實現非線性特征嵌入。然而，MDS不能擴展到大規模數據，因為必須計算成對距離來最小化目標函數。

t-SNE：是一種將距離轉換為概率的方法。t-SNE構造一個與原始空間及映射后的低維空間中數據點之間的相似性相關的概率分布，然后最小化兩個分布之間的Kullback-Leibler散度。t-SNE被廣泛應用于單細胞數據分析中的數據可視化。

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

Figure 2. Seurat CCA數據整合示意圖

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

相似度及核函數

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

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

審核編輯：劉清