崔 璐, 劉桂鋒

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

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

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

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

給定一個單細(xì)胞RNA-seq基因表達(dá)矩陣E∈n×m, 其中包含n個基因和m個細(xì)胞樣本. 設(shè)該表達(dá)矩陣中缺失值集合為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(·) 表示缺失值補(bǔ)全函數(shù);E′為補(bǔ)全后的單細(xì)胞RNA-seq基因表達(dá)矩陣. 上述模型與傳統(tǒng)數(shù)據(jù)補(bǔ)全模型的不同: 單細(xì)胞RNA-seq基因表達(dá)矩陣中的缺失值位置和數(shù)量均未知, 因此, 對該矩陣進(jìn)行補(bǔ)全, 傳統(tǒng)缺失值補(bǔ)全方法已不再適用. 單細(xì)胞RNA-seq數(shù)據(jù)缺失模型需同時解決3個問題, 即矩陣中缺失值元素的位置、 數(shù)量和補(bǔ)全值.

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

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

(3)

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

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

(4)

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

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

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

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

3) 迭代求解非負(fù)矩陣: 根據(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)方法繪制細(xì)胞分型圖譜.

3 實(shí)驗(yàn)和分析

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

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

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

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