馬敬山，毛玉杰，張 杉，3

(1.石家莊郵電職業(yè)技術(shù)學(xué)院，石家莊 050022； 2.燕山大學(xué) 理學(xué)院，河北 秦皇島 066004；3.中車唐山機車車輛有限公司，河北 唐山 063000)

許多研究者陸續(xù)完成了CCLE[1]、GDSC[2]等多項癌癥基因組計劃。研究表明，在高維基因表達數(shù)據(jù)和抗癌藥物反應(yīng)數(shù)據(jù)中，識別具有統(tǒng)計學(xué)和生物學(xué)意義的基因-藥物共模塊，有助于理解抗癌藥物的分子機制，篩選潛在的藥物靶點。Kutalik等使用NCI-60數(shù)據(jù)[3]完成了初步實驗，獲得了Ping-Pong算法[4]。Chen等提出了多矩陣分解算法(NetNMF)[5]，基于基因組數(shù)據(jù)構(gòu)建的相似網(wǎng)絡(luò)矩陣，使用三元非負矩陣分解，來尋找公共模塊和模塊之間的關(guān)聯(lián)。Wang等在非負矩陣分解模型的基礎(chǔ)上，通過向分解因子添加L1-范數(shù)規(guī)范約束，提出了RNMF算法[6]。Zhang等提出了具有網(wǎng)絡(luò)正則性約束的改進的聯(lián)合非負矩陣分解算法(JNMF)[7-8]。

受上述算法的啟發(fā)，提出了一種帶相似約束的稀疏聯(lián)合非負矩陣分解算法(SSJNMF)，并將其應(yīng)用于基因藥物共模塊識別的GDSC數(shù)據(jù)集。

1 數(shù)據(jù)和算法

1.1 數(shù)據(jù)源和預(yù)處理

下載同一細胞系對應(yīng)的最新基因表達數(shù)據(jù)和藥物反應(yīng)數(shù)據(jù)，發(fā)現(xiàn)藥物反應(yīng)數(shù)據(jù)缺失值，數(shù)據(jù)預(yù)處理如下：

刪除缺失值大于30%的列，其余203列對應(yīng)203種藥物，用mice包進行基因填充[9]，以獲得完整的藥物反應(yīng)矩陣。

記基因表達矩陣為G1∈R915×17 737,記藥物響應(yīng)矩陣為G2∈R915×203，利用皮爾遜相關(guān)系數(shù)作為工具，獲得基因相似性矩陣X11、藥物相似性矩陣X22、基因與藥物相似性矩陣X12。取矩陣X11、X22和X12的絕對值，同時將輸入數(shù)據(jù)完成[0,1]均值處理[10]，以保證輸入數(shù)據(jù)(G1,G2)和相似性數(shù)據(jù)處于相同的數(shù)量級。

1.2 算法構(gòu)造

基于文獻[7]中識別多維模塊的思想，將基因相似性數(shù)據(jù)、藥物相似性數(shù)據(jù)和基因藥物相似性數(shù)據(jù)添加到聯(lián)合非負矩陣分解算法中，完成“帶相似約束的稀疏聯(lián)合非負矩陣分解算法”的構(gòu)造，簡稱SSJNMF，根據(jù)降維后的數(shù)據(jù)確定模塊數(shù)，計算后限制分解因子，最后篩選模塊數(shù)據(jù)，目標函數(shù)如下：

1.3 算法求解過程

Lee和Seung提出的乘法器更新迭代算法[11]優(yōu)化SSJNMF模型，以保證變量B和C的凸性，并獲得全局最優(yōu)解：

擴展SSJNMF模型如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

SSJNMF算法的求解過程：

第一步：用B、C1和C2為初始值，代入模型第一部分，得到SSJNMF的初始迭代矩陣。

第二步：按照(6)、(7)和(8)交替更新矩陣B,C1和C2。

第三步：重復(fù)第二步，終止條件為誤差小于10-6或次數(shù)到達500次。

2 結(jié)果分析與比較

2.1 SSJNMF算法的識別過程及識別結(jié)果

識別的第一步是分析統(tǒng)計顯著性，通過比較常用模塊中基因與藥物信息的相關(guān)性，分析算法識別結(jié)果的統(tǒng)計顯著性，具體步驟如下：

一是以共模塊為單位，從G1和G2中篩選出子矩陣sG1和sG2，按列求皮爾遜相關(guān)系數(shù)：

系數(shù)之和表示為：

二是將基因表達數(shù)據(jù)和藥物應(yīng)答數(shù)據(jù)按列重新排列，選岀1 000個維度相同的隨機共模塊，Srand用于表示每個隨機共模塊各列相關(guān)系數(shù)之和。

三是在1 000個Srand所形成的分布下，估計岀Sreal所對應(yīng)的概率分布統(tǒng)計值P1 (大于Sreal的比例)。若P1<0.05，可以認為它不服從1 000個分布的置信區(qū)間之外的隨機分布，即結(jié)果是可行的和可解釋的，相反共模塊服從隨機分布。此外，非參數(shù)秩和檢驗的實驗結(jié)果也被用來推斷服從的分布，將其置于分布之下并計算出相應(yīng)的概率分布統(tǒng)計值P2，如果P2<0.05/k，實驗結(jié)果具有統(tǒng)計學(xué)意義。

辨識的第二步是確定參數(shù)，并通過控制變量法調(diào)整參數(shù)來選擇最優(yōu)參數(shù)。表1為不同參數(shù)組合下的統(tǒng)計結(jié)果，其中，Num1為非空且不服從隨機分布的共模塊個數(shù)，Num2為藥物模塊為空的共模塊個數(shù)，Num3為非空且通過非秩和檢驗的共模塊個數(shù)。結(jié)果表明，當(dāng)λ1=0.1、λ2=150、λ3=0.5時，識別出的具有統(tǒng)計意義的共模塊最多。

表1 不同參數(shù)組合下的統(tǒng)計結(jié)果Tab.1 Statistical results under the combination of different parameters

表2為當(dāng)λ1=0.1、λ2=150、λ3=0.5時，調(diào)整γ1、γ2、k的部分統(tǒng)計結(jié)果，其中Num4為富集分析后有生物意義的基因模塊個數(shù)，通過綜合分析P1[P1=Num1/(70-Num2)]，P2[P2 =Num3/(70-Num2)]以及具有生物意義的共模塊所占的比例P3[P3=Num4/(70-Num2)]及綜合表1，可以得岀以下結(jié)論：當(dāng)參數(shù)λ1=0.1、λ2=150、λ3=0.5、γ1=0.5、γ2=0.5、k=70時，SSJNMF算法識別岀的70個基因-藥物共模塊為最優(yōu)結(jié)果(模塊個數(shù)由閾值T決定，當(dāng)T=3.7時，識別的70個共模塊指標最優(yōu)[12])。

表2 固定λ1、λ2、λ3調(diào)整γ1、γ2、k的統(tǒng)計結(jié)果Tab.2 Statistical results of adjustment of γ1、γ2、k under fixed λ1、λ2、λ3

2.2 共模塊的生物功能分析

借鑒Mao等人[13]提岀的生物功能分析方法，分析SSJNMF算法識別岀的68個有意義的共模塊，結(jié)果表明，有60個模塊的藥物個數(shù)大于1。表3列舉了第18和第46個共模塊的生物功能分析結(jié)果。表中，ID:共模塊序號；G:對應(yīng)模塊中的基因；D:藥物個數(shù)。

表3 部分基因藥物共模塊的生物功能分析結(jié)果Tab.3 Results of biological function analysis of some gene-drug common models

生物功能相關(guān)分析結(jié)果表明，每個基因模塊的GO生物功能項所富集的生化過程與相應(yīng)藥物模塊中藥物靶向的信息之間存在著很強的相關(guān)性。

2.3 三類算法的識別結(jié)果比較

為了更好地說明SSJNMF的識別性能，在同一數(shù)據(jù)集上執(zhí)行了NetNMF[6]和 JNMF[8]，識別岀各自的70個基因-藥物共模塊。結(jié)果如表4?？梢钥磳?，SSJNMF得到的P1和P2均比NetNMF和JNMF高，說明SSJNMF識別的共模塊有更強的統(tǒng)計意義。

表4 三種算法識別的共模塊的統(tǒng)計意義對比分析Tab.4 Comparative analysis of statistical significance of common modules of three algorithm identification

3 結(jié)論

SSJNMF算法是有效的基因-藥物共模塊的識別工具，識別結(jié)果具有很強的統(tǒng)計意義和生物意義。基因模塊的GO生物功能項所富集的生化過程與對應(yīng)的藥物模塊中藥物所靶向的信息之間具有強關(guān)聯(lián)性。SSJNMF篩選岀的部分結(jié)果可能存在信息相似情況，后期可以進行合并處理，以進一步提高共模塊識別的有效性。