基于loess 回歸加權(quán)的單細胞RNA-seq 數(shù)據(jù)預處理算法

高美加

(哈爾濱工業(yè)大學 計算機科學與技術學院, 哈爾濱150040)

0 引 言

相比于傳統(tǒng)的細胞測序方法,單細胞RNA-seq的測序數(shù)據(jù)提供了研究細胞異質(zhì)性和基因差異表達的機會,但是單細胞RNA-seq 通常表現(xiàn)出比來自大量細胞群的RNA-seq 數(shù)據(jù)更高水平的噪聲和更多的零值。 scRNA-seq 數(shù)據(jù)的計算分析包括質(zhì)量控制、定位、定量、標準化、聚類幾個步驟,用于鑒定差異表達的基因。 上游的步驟可能對結(jié)果產(chǎn)生實質(zhì)性的影響。 scRNA-seq 的大多數(shù)分析,如基因差異表達分析、細胞類型特異性基因的鑒定、分化軌跡的重建等,都依賴于基因表達測量的準確性。 目前,對單細胞RNA-seq 得到的矩陣的預處理方法主要是對矩陣進行插值,以此減輕過多零值對后續(xù)的影響。此方法利用單細胞基因表達數(shù)據(jù)的結(jié)構(gòu),通過利用相關細胞或者基因表達之間的相似性來校正基因的表達量［1］。 例如：sclmpute 是利用混合模型來定位可能的缺失值,之后對其進行插補；MAGIC 和SAVER 是對矩陣去噪,生成一個新的矩陣,以上都是通過線性來對矩陣進行去噪［1］。 另外,也有一些使用神經(jīng)網(wǎng)絡的方法來進行插補的算法,使用自編碼器可以通過無監(jiān)督的方式,最小化重建數(shù)據(jù)和原始數(shù)據(jù)之間的誤差,來進行非線性的差值,同時也可以進行有效的數(shù)據(jù)壓縮,例如：DCA 算法。

另外,在單細胞RNA-seq 數(shù)據(jù)的降噪方法中常用的有基因篩選和降維。 基因篩選即篩選出在細胞中表達量變化大的基因,這樣可以去除低變化高表達量基因?qū)罄m(xù)分析的影響［2］。 在RNA-seq 數(shù)據(jù)分析中常用的降維方式有PCA、KPCA 和t-SNE 等。

常用的基因篩選算法有Seurat 包里的disp、vst、mvp 等。 但是,雖然一些管家基因的表達信息對于細胞的分類并不能起到什么關鍵的作用［3］,降低這些基因的影響,可能會對后續(xù)分析(細胞聚類等)有一些提升。 在單細胞表達矩陣的預處理過程中,通常會先回歸擬合基因在細胞中表達量的標準差與平均值的變化曲線來對基因進行篩選,但是這樣會損失一部分信息,從而影響后續(xù)的分析質(zhì)量。

基于以上問題,本文提出一種基于Loess 回歸加權(quán)的單細胞轉(zhuǎn)錄組數(shù)據(jù)預處理算法,通過Loess回歸曲線定量計算基因表達偏移水平,并基于偏移水平構(gòu)造基因加權(quán)系數(shù),達到基因軟篩選與數(shù)據(jù)降噪的目的。 本文選擇6 組單細胞RNA-seq 數(shù)據(jù)從可視化和聚類兩方面對算法預處理效果進行測試,實驗證明該方法可以有效降低低質(zhì)量基因?qū)Ψ治鲞^程的影響,提升下游分析的精準水平,顯示出較好應用價值。

1 預處理方法研究

1.1 回歸加權(quán)

圖1 Pollen 數(shù)據(jù)集Fig. 1 Pollen dataset

通過量化基因在每個細胞里表達的高變異度,對表達量矩陣進行加權(quán)來降低變化度的基因?qū)罄m(xù)分析的影響。 在此使用局部加權(quán)回歸(LOESS)擬合基因在細胞中的表達量的標準差與平均值的變化曲線,使用實際的標準差和預測的標準差之間的差值作為每個基因的權(quán)重,然后生成新的表達矩陣,如式(1)～式(4)所示。

其中,meanisdi為基因i 的表達的平均值和標準差,xij為矩陣中的元素,x′ij為新生成的表達值。 以Pollen 數(shù)據(jù)集為例,Pollen 數(shù)據(jù)集經(jīng)變換后,如圖1所示。 在圖1(a)中,可以看出,經(jīng)PCA 降維后,預處理后數(shù)據(jù)集的可視化效果要好一些,有幾類細胞在圖中被有效分離開,圖1(b)是經(jīng)過TSNE 降維后的效果,各個簇也更聚集一些。

1.2 標準化

由于技術原因,單細胞RNA-seq 數(shù)據(jù)中基因表達顯示出明顯的細胞差異,可能是由于生物學和技術上的雙重原因造成的。 在此使用了Hafemeister等人提出的一個標準化方法來降低測序深度對基因表達造成的影響,公式(5)和公式(6)如下。

其中：m 為細胞j 中基因的總的表達量,分別對每個基因在細胞中的表達量和細胞總的表達量做線性回歸,計算出每個基因在細胞中期望的表達值,根據(jù)新的矩陣使用上述的方法做預處理。 圖2 中,以Zeisel 數(shù)據(jù)集為例,圖2(a)為未處理過的數(shù)據(jù),圖2(b)為使用LOESS 加權(quán)處理過的新的矩陣,圖2(c)為標準化后加權(quán)處理形成矩陣的可視化效果。 可以看出,在圖2(a)中形成了3 個比較大的簇,在圖2(b)和(c)中都有其它的簇分裂出來,可視化效果較好。

2 實驗結(jié)果

分別從可視化效果和無監(jiān)督聚類效果兩方面來對預處理效果進行評價。 在這里選取了6 個數(shù)據(jù)集：pollens,Biase,Yan,Goolams,Deng,Zeisel。 其中Zeisel 數(shù)據(jù)集中的標簽為SC3 算法得出的標簽。 數(shù)據(jù)來源：https：／ ／hemberg-lab.github.io／scRNA.seq.datasets／。 本實驗使用R 語言中的scater 包進行大部分的單細胞RNA-seq 數(shù)據(jù)分析。

圖2 Zeisel 數(shù)據(jù)集Fig. 2 Zeisel dataset

2.1 對可視化效果影響

選取了在單細胞RNA-seq 分析中比較常用的3種可視化方法來進行對比實驗,這3 種方法分別是：PCA、TSNE、UMAP。 其中TSNE 和UMAP 由于實現(xiàn)過程中有隨機性,所以重復了500 次實驗進行對比。使用輪廓系數(shù)(Silhouette Coefficient)來對圖中同種細胞的聚集程度和不同種細胞的離散程度進行量化。輪廓系數(shù)就是針對樣本空間中的一個特定樣本,計算它與所在聚類其它樣本的平均距離a,以及該樣本與距離最近的另一個聚類中所有樣本的平均距離b,該樣本的輪廓系數(shù)為(b - a)／max(a, b),將樣本空間中所有樣本的輪廓系數(shù)取算數(shù)平均值,作為聚類劃分的性能指標s。 在前兩個維度里,類間距離越大,類內(nèi)距離越小,就認為這個可視化效果是好的。

可視化實驗結(jié)果如圖3 所示(每張圖中前3 個為TSNE 可視化的輪廓系數(shù)對比,后3 個為UMAP 可視化后輪廓系數(shù)對比；這三列分別為原始矩陣、標準化后又加權(quán)的矩陣和只回歸加權(quán)的可視化輪廓系數(shù)對比。 圖3(b)從左到右,從上到下分別為：biase、deng、gool、pollen、yan、zeisel 數(shù)據(jù)集)。 通過比較發(fā)現(xiàn),對于PCA 來說,在biase、deng 和yan 數(shù)據(jù)集中標準化之后再進行預處理的效果較好,在pollen 和zeisel 數(shù)據(jù)集中直接進行預處理效果較好。 對于tSNE 和UMAP算法來說,biase、gool、zeisel 數(shù)據(jù)標準化之后再進行預處理效果較好,在pollen、yan、deng 這3 個數(shù)據(jù)集中直接進行預處理效果更好。 從實驗結(jié)果來看,基因表達矩陣經(jīng)算法處理過之后,經(jīng)過降維,前兩維中數(shù)據(jù)同類樣本更集中,不同類樣本之間也更加分散,它對后續(xù)的可視化效果是有一定提升的。

圖3 可視化效果對比Fig. 3 Visual effect comparison

2.2 對無監(jiān)督聚類的影響

選取3 個常用的單細胞測序數(shù)據(jù)的聚類算法：SC3、SIMLR、和Seurat。 使用F1 - score 來對聚類結(jié)果進行評價。 F1 - score 具體定義為公式(9)：

通過公式(7) 和(8) 計算每個類別下的precision 和recall：

其中, TP (True Positive)預測答案正確； FP(False Positive)錯將其他類預測為本類； FN (False Negative)本類標簽預測為其他類標。 最后通過公式(10)計算各個類別下的F1 - score 的平均值：

F1-score 是精確率和召回率的調(diào)和平均數(shù),最大為1,最小為0。

聚類分析結(jié)果如表1 所示,其中l(wèi)m＋loess 為經(jīng)標準化后的回歸加權(quán)方法。 經(jīng)過比較發(fā)現(xiàn)兩種預處理方法對這三種聚類方法的實驗結(jié)果都有一定的提升作用,其中LM＋LOESS＋SC3、SC3、LOESS＋SIMLR在多數(shù)數(shù)據(jù)集中表現(xiàn)都比較好,說明回歸加權(quán)的方法是對后續(xù)的無監(jiān)督聚類分析有一定的提升作用。

表1 聚類結(jié)果比較Tab. 1 Comparison of clustering results

3 結(jié)束語

本文提出的基于Loess 回歸加權(quán)單細胞RNA-seq數(shù)據(jù)的預處理算法。 可以看出,在一些數(shù)據(jù)集中,預處理之后的可視化和無監(jiān)督聚類過程都有一定的提升作用,數(shù)據(jù)經(jīng)過PCA 或者t-SNE 降維后,經(jīng)處理后的數(shù)據(jù)同類細胞間往往表現(xiàn)的更加聚集,不同類之間更加分散,這同樣會加強后續(xù)的聚類效果,使聚類算法表現(xiàn)更好。 gool|、deng、yan 等人的數(shù)據(jù)集經(jīng)過預處理后,聚類結(jié)果準確度明顯有了很大的提升。 但是此算法也有一定的局限性,預處理之后的數(shù)據(jù)產(chǎn)生的值并不符合矩陣中元素為基因在細胞中的表達量這一定義,不利于差異表達基因等下游分析的進行,還有待進行一些分析與改進。