基于loess 回歸加權的單細胞RNA-seq 數據預處理算法

高美加

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

0 引 言

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

另外,在單細胞RNA-seq 數據的降噪方法中常用的有基因篩選和降維。 基因篩選即篩選出在細胞中表達量變化大的基因,這樣可以去除低變化高表達量基因對后續分析的影響［2］。 在RNA-seq 數據分析中常用的降維方式有PCA、KPCA 和t-SNE 等。

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

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

1 預處理方法研究

1.1 回歸加權

圖1 Pollen 數據集Fig. 1 Pollen dataset

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

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

1.2 標準化

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

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

2 實驗結果

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

圖2 Zeisel 數據集Fig. 2 Zeisel dataset

2.1 對可視化效果影響

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

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

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

2.2 對無監督聚類的影響

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

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

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

F1-score 是精確率和召回率的調和平均數,最大為1,最小為0。

聚類分析結果如表1 所示,其中lm＋loess 為經標準化后的回歸加權方法。 經過比較發現兩種預處理方法對這三種聚類方法的實驗結果都有一定的提升作用,其中LM＋LOESS＋SC3、SC3、LOESS＋SIMLR在多數數據集中表現都比較好,說明回歸加權的方法是對后續的無監督聚類分析有一定的提升作用。

表1 聚類結果比較Tab. 1 Comparison of clustering results

3 結束語

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