大規模單細胞轉錄組測序數據的聚類方法比較*

朱曉姝,蒙 霜,龍法寧

(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數據會導致更大的內存需求和計算時間開銷問題,進一步研究基于數據分布的有效去噪方法、適合大規模數據的圖神經網絡降維方法,以提高數據質量并準確度量細胞間相似性,在細胞類型識別時加強生物可解釋性,提升細胞類型識別和下游分析的性能等都是以后的重要工作。