單細胞轉(zhuǎn)錄組數(shù)據(jù)批次效應評測的研究進展

李小剛,高 正,陳佳鋒,馮山入,傅修濤,丁振斌,3*

(1.復旦大學附屬中山醫(yī)院 肝臟外科與移植外科,復旦大學肝癌研究所,上海 200032;2.教育部癌變與侵襲原理重點實驗室,上海 200032;3.上海市徐匯區(qū)中心醫(yī)院,復旦大學附屬中山徐匯醫(yī)院,上海 200032)

單細胞轉(zhuǎn)錄組測序(Single cell RNA sequencing,ScRNA seq)作為新興生物技術(shù),在解析組織復雜性、細胞異質(zhì)性、基因表達調(diào)控等方面發(fā)揮重要作用,很大程度上彌補了普通轉(zhuǎn)錄組測序(Bulk RNA sequencing)的缺陷[1]。批次效應(Batch effects,BEs)是高通量生物學數(shù)據(jù)中源于技術(shù)因素的變異[2],單細胞數(shù)據(jù)通常需要多次實驗來獲得,不同批次實驗的試劑、實驗儀器、實驗員、單細胞捕獲時間,甚至是技術(shù)平臺,均可能存在差異。因此,這些因素也就構(gòu)成了單細胞數(shù)據(jù)的批次,數(shù)據(jù)整合過程中可能會干擾感興趣生物學變異的發(fā)現(xiàn)和解釋,如何有效去除批次效應在單細胞數(shù)據(jù)分析中有充分的必要性。

然而,單細胞數(shù)據(jù)批次效應較為復雜,且和具體的研究有關(guān)[3-4],有效處理能增加數(shù)據(jù)價值,反之則會可能會導致假陽性或假陰性分析結(jié)果[5]。批次效應可能是高度非線性的,常常和真正的生物學變異相互混雜,正確糾正不同批次單細胞數(shù)據(jù)的批次效應,同時正確保留關(guān)鍵生物學變異常常較為困難。為了解決上述問題,此前用于處理芯片數(shù)據(jù)批次效應的工具,包括Combat[6]和Limma等[7],也用于單細胞數(shù)據(jù)的批次效應糾正。由于單細胞數(shù)據(jù)的具有高缺失率(Dropout)以及基因捕獲隨機性的特點[8],和芯片數(shù)據(jù)存在顯著不同,因此適用于單細胞數(shù)據(jù)的批次效應處理算法也在不斷涌現(xiàn)。其中,代表性算法的包括通過識別相互最近鄰來整合不同的數(shù)據(jù)集的MNNs(Mutual nearest neighbours)[9],建立在MNNs基礎上的Scanorama[10]和BBKNN[11],整合于Seurat包中的MultiCCA算法[12],處理多種來源的批次效應具有顯著優(yōu)勢的Harmony[13],將不同數(shù)據(jù)集間變異全部歸于技術(shù)因素的LIGER[14],以及新型的深度學習算法[15],均在不同程度上去除了批次效應,促進了生物學規(guī)律的發(fā)現(xiàn)。

目前批次效應的處理已經(jīng)是單細胞數(shù)據(jù)分析的常規(guī)流程(單細胞數(shù)據(jù)分析流程見圖1),但批次效應的評估是有效去除批次效應的基礎,也常常是被忽視的一步[16]。實際上,選擇合適的指標評測批次效應的來源、對數(shù)據(jù)變異的貢獻度和對數(shù)據(jù)分析的影響有利于判斷處理批次效應的必要性及選擇合適的處理算法。如今跨平臺、跨物種以及多模式單細胞數(shù)據(jù)的整合,包括單細胞表觀組、基因組、轉(zhuǎn)錄組、蛋白組等[17],在研究生物學規(guī)律及疾病發(fā)生發(fā)展機制方面顯示出巨大優(yōu)勢[18]。同時,數(shù)據(jù)量和數(shù)據(jù)來源的增加也使得批次效應更加復雜,準確地評估批次效應的來源、批次效應處理后的效果,對于數(shù)據(jù)整合尤為關(guān)鍵[19]。

圖1 單細胞數(shù)據(jù)分析流程Fig.1 Analysis pipeline of single cell data

因此,本篇綜述聚焦單細胞轉(zhuǎn)錄組數(shù)據(jù)的批次效應,依次論述了單細胞批次效應與普通轉(zhuǎn)錄組的區(qū)別,目前常用的評測算法的特點,最后總結(jié)了目前面臨的挑戰(zhàn)和未來發(fā)展方向。

1 單細胞轉(zhuǎn)錄組和普通轉(zhuǎn)錄組批次效應的區(qū)別

分析目的及數(shù)據(jù)結(jié)構(gòu)的不同是構(gòu)成兩者批次效應差異的基礎。普通轉(zhuǎn)錄組主要目的是計算樣本間的差異基因及分子分型,在宏觀層面解析其基因表達改變[20]。單細胞轉(zhuǎn)錄組則是以細胞群為研究單位,通過降維、聚類和細胞類型注釋來識別特定的、同質(zhì)性的細胞群,并在細胞群及基因?qū)用孢M行探索性分析[21-22]。數(shù)據(jù)結(jié)構(gòu)方面,普通轉(zhuǎn)錄組數(shù)據(jù)是低維數(shù)據(jù),數(shù)據(jù)結(jié)構(gòu)較簡單,無需做降維處理;而單細胞轉(zhuǎn)錄組數(shù)據(jù)為高維度數(shù)據(jù),下游分析前需要多步驟降維處理。

基因缺失率(Dropout rates)是兩者最大的不同。高質(zhì)量普通轉(zhuǎn)錄組測序產(chǎn)生的缺失率低于20%,而基于微液滴或者微孔技術(shù)的單細胞轉(zhuǎn)錄組的缺失率可達80%。即使是測序深度高并且支持全長測序的Smart-seq2技術(shù),缺失率也達到了50%。然而,基因缺失的產(chǎn)生并不完全隨機,具有一定的基因偏倚性(Gene-based bias)、細胞偏倚性(Cell-based bias)和批次偏倚性(Batch-based bias)[2, 23]。

此外,普通轉(zhuǎn)錄組測序同質(zhì)性較高,批次效應通常涉及不同研究來源及不同平臺的數(shù)據(jù),處理批次效應的流程較為固定。而單細胞轉(zhuǎn)錄組測序則涉及多個樣本、多種細胞類型和成千上萬的細胞,甚至百萬級細胞量[24],批次效應的評測和處置都相對復雜且缺乏固定的流程。不同樣本的細胞類型差異可能很大,甚至部分細胞類型是某些樣本特有的,因此這種批次來源變異和生物來源變異的相互混雜是普通轉(zhuǎn)錄組未曾碰到過的難題。同時,為了消除多模式的單細胞數(shù)據(jù)的整合分析中批次效應的影響,則需要對來源于不同平臺和不同樣本的數(shù)據(jù)進行多步驟、多來源批次效應的評測和處置。

綜上,由于兩者測序技術(shù)、分析目的、數(shù)據(jù)結(jié)構(gòu)以及數(shù)據(jù)量之間的差異,適用于普通轉(zhuǎn)錄組的工具并不能直接用于單細胞轉(zhuǎn)錄組數(shù)據(jù)。因此,在單細胞多組學時代,開發(fā)新型的批次效應評測及處置算法來整合海量單細胞數(shù)據(jù),是正確揭示生物學規(guī)律和疾病發(fā)生發(fā)展機制的基礎。

2 批次效應的定性評測:特征選擇法

特征選擇可以評估批次效應對差異基因的影響。目前常用的差異基因工具包括t-test,limma[25],DESeq2[26]和edgeR[27],通過比較批次效應糾正前后的差異基因、差異基因的功能富集與某個已知表型的一致性來實現(xiàn)對批次效應的評測。特征選擇并不限于基因或者蛋白質(zhì)、代謝物等生物單位,也適用于主成分(Principle components,PCs),通過可視化解釋度最大的幾個主成分是否與細胞群或批次存在更強關(guān)聯(lián)來評測批次效應,尤其是批次效應處理后。

主成分分析(Principle component analysis,PCA)聯(lián)合散點圖是常用的可視化方式[16]。但PCA對批次效應的評測的前提是前兩個或者前三個主成分可以解釋批次效應來源的變異。PCA是將數(shù)據(jù)的變異分解為正交因子,各個主成分之間相互獨立,意味著用于可視化批次效應的PCA散點圖取決于所使用的PC,并且不保留單細胞高維數(shù)據(jù)的局部結(jié)構(gòu)。因此,通過PCA解析數(shù)據(jù)結(jié)構(gòu)需要較多的數(shù)據(jù)分析技巧作為支撐。PCA是將高維數(shù)據(jù)的變異分解為正交因子而不是將高維數(shù)據(jù)嵌入到低維流形中,那么由K個因子構(gòu)成的單細胞數(shù)據(jù)變異就需要K個主成分來解釋,而K值實際上并不知曉。

t-SNE(t-distributed Stochastic neighbor embedding)[28]和UMAP(Uniform manifold approximation and projection)[29]是PCA的拓展和延伸,都是將單細胞高維數(shù)據(jù)嵌入低維度、非線性流形的算法[30]。與PCA不同,t-SNE和UAMP將高維數(shù)據(jù)嵌入低維空間并保留高維數(shù)據(jù)的局部結(jié)構(gòu),有利于直觀展示樣本分布和潛在的細胞群,在解析單細胞數(shù)據(jù)中的異質(zhì)性和批次效應來源的變異方面具有優(yōu)勢。然而,t-SNE對單細胞高維數(shù)據(jù)的全局結(jié)構(gòu)賦予的權(quán)重很小,使得降維后低維流形中細胞群間的相對空間距離和空間位置沒有意義,因此利用t-SNE中的空間位置信息來推斷細胞群之間的相似性或者其他任何關(guān)系均不可行。UMAP與t-SNE類似,但運行的速度更快,在保留單細胞高維數(shù)據(jù)的全局結(jié)構(gòu)方面表現(xiàn)更好。

UMAP和t-SNE將單細胞高維數(shù)據(jù)嵌入到低維流形中時,依賴對細胞間距離的不同定義,而這些距離的定義具有一定的探索性,對于一些單細胞數(shù)據(jù)并不適用。因此應該慎重對待降維結(jié)果,Charia認為,即不加區(qū)別地應用這些算法可能會造成數(shù)據(jù)結(jié)構(gòu)在多個維度上的改變,雖然保留了大部分的數(shù)據(jù)結(jié)構(gòu),但細胞、最近鄰和細胞類型之間的量化關(guān)系被嚴重扭曲[31]。同時,得到的結(jié)果并不穩(wěn)健,會隨著參數(shù)的調(diào)整而發(fā)生顯著的變化[32]。因此,對使用UMAP和t-SNE可視化的細胞群做生物學分析時,應結(jié)合其他局部指標,包括標記基因(Marker genes)的特異性、效應值大小或概率推斷的結(jié)果[31]。

由于細胞被嵌入的空間缺乏生物學解釋度,t-SNE和UMAP可能無法揭示關(guān)于批次效應的可分析信息[32]。t-SNE和UMAP都是建立在鄰接圖(Neighbor graph)基礎上的算法,并不保證單細胞高維數(shù)據(jù)中細胞群間的距離關(guān)系被正確保留。同時,將單細胞的高維數(shù)據(jù)變異嵌入低維空間后,一些細微、不顯著的批次效應干擾可能無法有效識別。

綜上,上述方法均有一定局限性,無偏倚的實驗設計對于減少批次效應的干擾仍是關(guān)鍵。若批次來源和細胞群來源的變異完全混雜,使用上述方法顯然不合適,但當兩者是正交關(guān)系時,PCA仍是一種系統(tǒng)評測批次效應的方法,通過對批次高度相關(guān)的主成分進行系統(tǒng)性檢驗來針對性地糾正或避免批次效應相關(guān)的數(shù)據(jù)特征,而UMAP和t-SNE無法實現(xiàn)這點[33]。

3 批次效應的定量評測:新型算法

單細胞測序技術(shù)的發(fā)展催生出評測批次效應的新方法。盡管PCA散點圖是一種評測批次效應的經(jīng)典可視化算法,但新型算法工具已經(jīng)展現(xiàn)出極大優(yōu)勢。

較為流行的kBET(k-nearest Neighbor batch effect test)由無監(jiān)督機器學習算法(ML,Machine learning)改造而來[34],計算單個樣本的每個細胞在不同批次間的最近鄰,并評測每個批次間的混合方式及程度。該算法原理較為簡單:若檢測到每個細胞的最近鄰在不同批次間是均勻分布的,那認為是不存在批次效應的。然而,如果批次間的細胞類型組成是不均勻的(比如某細胞類型的大多數(shù)細胞均來源于某一批次樣本),kBET是不適用的。同樣,數(shù)據(jù)存在明顯離群值時(某一個細胞群內(nèi)部,部分細胞所處的空間位置明顯偏離該細胞群的主體空間位置)或數(shù)據(jù)有明顯異質(zhì)性時(多批次的樣本內(nèi)部有大量的細胞亞群,并且不同批次間的細胞亞群組成的異質(zhì)性較大),kBET也無法有效地評測批次效應。LISI[13](Local inverse simpson’s index)算法與kBET的原理類似,不僅評估批次間的細胞混合程度,還計算細胞譜系間的混雜程度。

另外一種方法是ASW(Average silhouette width),同樣由無監(jiān)督機器學習算法延伸而來,是用于評估高質(zhì)量細胞群的聚類驗證指數(shù)[35]。ASW通常與無監(jiān)督聚類的算法聯(lián)合使用(比如k-means)。在ASW中,由一群點聚集而成的輪廓是對高維數(shù)據(jù)的低維展示,代表一個細胞群。該算法通過細胞群內(nèi)部及細胞群間的緊湊度進行量化,將量化值作為該單細胞數(shù)據(jù)的整體分群質(zhì)量,然后判斷整體分群情況與和批次之間是否存在關(guān)聯(lián)。然而,ASW畢竟是對整體分群情況的一種量化,難以評估特定細胞群是否存在批次上的偏倚,并對離群值較敏感。除此之外,即使ASW提示細胞群的分群質(zhì)量較高并與批次無相關(guān)性,也不意味著批次效應可以忽略。即使一些單細胞數(shù)據(jù)的細胞群之間區(qū)分度很高(分群時,群與群之間的細胞數(shù)目差別很大或者群與群之間的距離很近),也可能會產(chǎn)生低ASW值。這種較小的 ASW值往往是錯誤地提示不同批次之間的細胞混合較好,導致細胞群與批次之間存在的關(guān)聯(lián)無法被正確識別。

與上述三種算法的原理差別較大,蘭德指數(shù)(RI,Rand index)評估對同一個數(shù)據(jù)應用不同聚類方法后所產(chǎn)生的聚類結(jié)果之間的關(guān)聯(lián)性[36]。RI通過評估每個樣本中的細胞如何聚到每個細胞群來判斷批次效應糾正后對數(shù)據(jù)的影響,即如果人為產(chǎn)生了與批次效應相關(guān)聯(lián)的細胞群,那就說明批次效應的糾正是對聚類的結(jié)果有影響的。而蘭德矯正指數(shù)(ARI, Adjusted rand index)則是對蘭德指數(shù)(RI)的一種隨機性矯正,通過評估不同方法對之間的比較來獲得他們之間的相似性,從而建立不同方法之間比較的基線標準。

然而,上述算法的前提條件是批次效應對聚類的影響要達到一定的程度才能在低維空間中識別出。如果批次間的樣本混合較好,上述算法也會計算出一個較好的評測結(jié)果,但如果批次效應對聚類的結(jié)果影響很小,通過這些工具難以識別批次效應。相比之下,PCA是評測批次效應對數(shù)據(jù)影響程度的一種魯棒性很高的方法[33],而在局部水平評估混合程度的算法(比如kBET)易受到糾正批次效應和保留有意義生物學信號之間的平衡的影響。

4 單細胞批次效應評測方法總結(jié)

目前評測單細胞批次效應最常用的是PCA、tSNE、UMAP這些可視化的方法,尤其是數(shù)據(jù)量較小、細胞類型較少時,并且他們在直觀展示分群結(jié)果以及發(fā)現(xiàn)新的細胞類型方面具有獨特的優(yōu)勢。然而,當單細胞數(shù)據(jù)量增大、細胞類型更復雜時,這些可視化方法的主觀性和局限性就較為突出了。差異基因分析(DEGs analysis)則更偏于對分析結(jié)果的直觀評價和對生物學意義保留情況的判斷。

不同的是,kBET、LISI、ARI和ASW都是相對客觀的指標,但其計算結(jié)果可能與可視化的結(jié)果存在差異較大甚至相互矛盾的情況。實際應用中,kBET和LISI的結(jié)果較為一致,可能是因為兩者都是對局部水平批次間混合程度的計算;反之,ARI和ASW都是對整體水平批次效應的評測,這可能是導致不同指標的結(jié)果不一致的原因。值得注意的是,這些評測指標對細胞類型純度和各批次混合程度是分開計算的,實際應用中還需要采用合適的策略將細胞類型純度和各批次混合程度進行整合。

評測單細胞批次效應的較好流程是,將可視化的方法和kBET等這些定量指標相互結(jié)合,綜合判定批次效應的來源、對數(shù)據(jù)分析的影響程度以及批次效應的糾正效果(單細胞數(shù)據(jù)批次效應評測方法總結(jié)和對比見表1)。

表1 單細胞數(shù)據(jù)批次效應評測方法總結(jié)和對比Table 1 Summary and comparison of batch effect evaluation methods for single cell data

5 批次效應評測面臨的挑戰(zhàn)

目前,單細胞中批次效應的評測還面臨以下挑戰(zhàn),并且這些挑戰(zhàn)很可能會延伸到新興起的技術(shù),比如說空間組學技術(shù)(Spatial omics technology)。

1)糾正批次效應的算法通常是對批次間的細胞進行相似度計算,然后對批次間的細胞做混合聚類來糾正批次效應。過度糾正批次來源的變異會埋沒部分生物學變異,導致某些樣本特異的細胞類型或者罕見細胞類型不能通過聚類鑒定出來。然而,PCA及其他算法工具尚無法充分評估批次效應糾正和生物意義保留之間的平衡關(guān)系。同時,現(xiàn)有工具還無法精確評測不均衡實驗設計(不同批次間樣本的細胞類型差異較大)帶來的批次效應,從而為不均衡實驗批次效應的去除提供依據(jù)。

2)糾正批次效應來源的變異與保留生物學變異之間的平衡?，F(xiàn)有的評測批次效應的算法工具不能量化批次效應和生物學變異各自對數(shù)據(jù)總變異的解釋度,以及計算出合適閾值在糾正批次效應的同時最大化保留生物學變異。PCA無法保留高維數(shù)據(jù)的局部結(jié)構(gòu)以及判斷生物學變異和批次的變異對數(shù)據(jù)的影響;t-SNE和UMAP雖然可以保留單細胞高維數(shù)據(jù)的局部結(jié)構(gòu),但是嵌入的低維空間并不能保證對生物學變異的解釋度。

3)對單細胞數(shù)據(jù)高缺失率造成的批次效應研究不足。由于測序深度、mRNA捕獲效率等的限制,大量基因的count是零值,但基因缺失的產(chǎn)生并不完全隨機,同樣有一定的偏倚性。某些基因更易產(chǎn)生零值,導致某些高表達這些基因的細胞類型及富集這些細胞類型的某批次樣本更易產(chǎn)生零值?；蛄阒凳桥涡闹匾獊碓?具體產(chǎn)生原因及數(shù)據(jù)分布特征尚不清楚,因此如何評測基因零值來源的批次效應并將其可視化,是評測批次效應的的重大挑戰(zhàn)之一。

4)不同單細胞批次效應評測指標可能出現(xiàn)結(jié)果不一致情況,并且與PCA、tSNE、UMAP這些可視化結(jié)果也可能差異較大。因此,如何將不同評測方法的評測方面進行統(tǒng)一化是面臨的挑戰(zhàn)之一。同時,如何很好地綜合細胞類型純度和各批次混合情況,也是評測工具未來需要解決的問題。

6 未來展望

多模式的單細胞數(shù)據(jù)不斷出現(xiàn)以及數(shù)據(jù)的開源性進一步增加,大數(shù)據(jù)整合分析已經(jīng)是單細胞技術(shù)領(lǐng)域的重要環(huán)節(jié),批次效應的評測仍然是生物信息學分析的熱點領(lǐng)域,必將直接影響到單細胞數(shù)據(jù)分析質(zhì)量和結(jié)果解讀。目前常用的評測批次效應的定性或者定量算法工具,比如說PCA,t-SNE和UMAP, kBET,ASW,RI等,有效地評測了批次效應,但在應用方面均存在一定的局限性。因此,未來隨著單細胞轉(zhuǎn)錄組測序技術(shù)的進步和分析流程的不斷優(yōu)化,基于機器學習或者深度學習的批次效應評測算法的開發(fā)和優(yōu)化,將有希望開發(fā)出整合批次效應可視化和定量指標、區(qū)分不同來源批次效應的影響、綜合細胞純度和各批次混合程度的評測工具。