崔 璐, 劉桂鋒

(1. 吉林大學(xué) 中日聯(lián)誼醫(yī)院醫(yī)療保險管理部, 長春 130033; 2. 吉林大學(xué) 中日聯(lián)誼醫(yī)院放射線科, 長春 130033)

單細胞DNA和RNA測序技術(shù)能幫助人們理解細胞功能和細胞異質(zhì)性, 是生物學(xué)和醫(yī)學(xué)領(lǐng)域的前沿技術(shù)[1]. 隨著單細胞測序數(shù)據(jù)的迅速增長, 數(shù)據(jù)分析與計算方法面臨新的挑戰(zhàn)[2]. 由于生物實驗捕獲率較低, 導(dǎo)致單細胞(RNA-seq)數(shù)據(jù)存在較高比例的缺失值問題. 文獻[3]研究表明, 在小鼠胚胎單細胞測序數(shù)據(jù)中缺失元素比例達30%. 如果直接刪除高比例的缺失值, 則將導(dǎo)致重要信息的缺失.

為了解決上述問題并更好地實現(xiàn)細胞分型和功能分析, 文獻[2]提出了一種兩階段機器學(xué)習(xí)識別估計方法, 該方法能有效識別出缺失值的位置并估計缺失值大小, 但該方法計算復(fù)雜度較高. 文獻[4]將缺失值定義為假陰性誤差, 并提出了一種基于Bayes方法的單細胞差異表達數(shù)據(jù)分析方法, 該方法需先假設(shè)缺失值滿足的分布, 但找到缺失值的準(zhǔn)確分布是一個難題[5]. 因此, 本文提出一種基于非負矩陣分解的單細胞RNA-seq數(shù)據(jù)缺失值補全算法. 該方法通過迭代求解法尋找最佳重構(gòu)基因表達矩陣.

1 單細胞RNA-seq數(shù)據(jù)缺失模型

給定一個單細胞RNA-seq基因表達矩陣E∈n×m, 其中包含n個基因和m個細胞樣本. 設(shè)該表達矩陣中缺失值集合為L, 0≤i

(i,j)=F(li,j，E),ei,j=0,

(1)

E′(i,j)=Q(L(i,j),E),

(2)

其中: (i,j)表示E中的任一缺失值;F(·)表示缺失值識別函數(shù);Q(·) 表示缺失值補全函數(shù);E′為補全后的單細胞RNA-seq基因表達矩陣. 上述模型與傳統(tǒng)數(shù)據(jù)補全模型的不同: 單細胞RNA-seq基因表達矩陣中的缺失值位置和數(shù)量均未知, 因此, 對該矩陣進行補全, 傳統(tǒng)缺失值補全方法已不再適用. 單細胞RNA-seq數(shù)據(jù)缺失模型需同時解決3個問題, 即矩陣中缺失值元素的位置、 數(shù)量和補全值.

2 迭代重構(gòu)基因表達矩陣算法

為解決單細胞RNA-seq數(shù)據(jù)中的缺失值問題, 文獻[2]利用式(1)和式(2)先后調(diào)用不同的機器學(xué)習(xí)模型求解, 較好地解決了缺失值問題, 但上述方法求解過程繁瑣且計算復(fù)雜度較高. 本文將非負矩陣分解算法引入到單細胞RNA-seq數(shù)據(jù)補全問題中求解. 對給定的單細胞RNA-seq基因表達矩陣E, 求解非負矩陣W∈n×t和H∈t×m為min{En×m-Wn×t×Ht×m}. 因此, 該求解過程是一個迭代求解W和H的過程. 其中W和H在非負矩陣分解求解過程中的更新規(guī)則[6]為

(3)

求出W和H后, 二者的矩陣點積為非負矩陣分解算法得到的結(jié)果矩陣. 結(jié)果矩陣中會估計出原始基因表達矩陣中缺失元素的表達值.

經(jīng)過非負矩陣分解算法補全后的重構(gòu)基因表達矩陣計算方法如下：

(4)

其中, 缺失值由非負矩陣分解結(jié)果矩陣補全, 而非缺失值位置的基因表達值E(u,v)保留原始表達值.

經(jīng)典非負矩陣分解求解前需結(jié)合人工經(jīng)驗輸入rank值, 該值大小會直接影響結(jié)果矩陣的誤差. 通過前期實驗表明, 即使rank值發(fā)生很小的變化, 都會對結(jié)果產(chǎn)生較大影響. 針對該問題, 本文提出一種迭代重構(gòu)基因表達矩陣算法, 算法步驟如下.

1) 算法輸入: 輸入單細胞RNA-seq基因表達矩陣E和初始誤差閾值ε0;

2) 初始化參數(shù): 根據(jù)數(shù)據(jù)設(shè)定rank值搜索范圍(2～θ)(50≤θ≤m)及當(dāng)前誤差ε=ε0;

3) 迭代求解非負矩陣: 根據(jù)式(3)和當(dāng)前rankt值, 迭代求解W和H;

4) 計算誤差εt+1, 若εt+1<εt則轉(zhuǎn)步驟3), 否則執(zhí)行步驟5);

5) 輸出結(jié)果矩陣: 根據(jù)式(4)計算E′和缺失值元素集合L;

6) 用t-SNE(t-distributed stochastic neighbor embedding)方法繪制細胞分型圖譜.

3 實驗和分析

為檢驗算法的有效性, 本文收集并整理了慢性粒細胞白血病單細胞測序數(shù)據(jù). 該數(shù)據(jù)來源于美國國家生物信息中心(NCBI)的基因表達數(shù)據(jù)庫, 數(shù)據(jù)檢索號為GSE76312[7]. 該原始數(shù)據(jù)共包含2 287個干細胞和23 384個基因, 本文選擇其中1 102個不含絡(luò)氨酸激酶抑制劑的細胞, 并利用偽發(fā)現(xiàn)率和差異倍數(shù)選取前234個差異表達基因. 將該數(shù)據(jù)分為5類： 急變期慢性粒細胞白血病細胞(BC-CML)、 慢性期慢性髓性白血病細胞(CP-CML)、 人紅白血病細胞系(k562)、 正常造血干細胞(normal)和前急變期慢性粒細胞白血病細胞(preBC).

圖1為慢性粒細胞白血病單細胞測序原始數(shù)據(jù)分型與迭代重構(gòu)基因表達矩陣數(shù)據(jù)細胞分析比較結(jié)果. 用t-SNE方法進行分型數(shù)據(jù)的可視化, 該方法是一種非線性降維和可視化方法[8], 目前已被廣泛應(yīng)用于腫瘤亞型分析[9]、 雷達目標(biāo)識別[10]及人物動作識別[11]等領(lǐng)域, 并取得了較好的可視化效果. 由圖1(A)可見, 該數(shù)據(jù)經(jīng)過細胞和差異表達基因的過濾篩選, 呈現(xiàn)了數(shù)據(jù)的可分性, 但5種類型的細胞簇被劃分的較分散, 且正常造血干細胞簇被分割開, 圖1(A)中左下側(cè)藍色橢圓簇類中包含了4種不同類別的細胞. 該分類結(jié)果表明, 慢性期慢性髓性白血病細胞、 前急變期慢性粒細胞白血病細胞和正常造血干細胞被混淆在一起. 由圖1(B)可見, 大部分正常造血干細胞聚集在右上角的藍色橢圓簇類中, 且該簇類中的前急變期慢性粒細胞白血病細胞和慢性期慢性髓性白血病細胞都明顯少于圖1(A). 因此, 迭代重構(gòu)基因表達矩陣數(shù)據(jù)細胞分型的結(jié)果更準(zhǔn)確.

圖1 原始數(shù)據(jù)與迭代重構(gòu)基因表達矩陣數(shù)據(jù)細胞分型結(jié)果比較Fig.1 Comparison of cell typing results between original data and iterative reconstruction matrix of gene expression data

綜上所述, 針對單細胞RNA-seq數(shù)據(jù)中存在缺失值的問題, 本文提出了一種單細胞RNA-seq數(shù)據(jù)缺失元素補全算法. 首先定義了單細胞RNA-seq基因表達矩陣數(shù)據(jù)缺失模型, 然后提出了迭代重構(gòu)基因表達矩陣算法, 并將該算法應(yīng)用到慢性粒細胞白血病單細胞測序原始數(shù)據(jù)中, 取得了更準(zhǔn)確的細胞分型結(jié)果.