毛玉杰， 魏東，2， 李玉雙

(1. 燕山大學 理學院， 河北 秦皇島 066004；2. 河北數(shù)港科技有限公司， 河北 秦皇島 066004)

近年來,國際上開展了一系列大型的藥物篩選實驗,其中,規(guī)模最大的兩個研究項目是癌癥細胞系百科全書(CCLE)[1]和癌癥基因組項目(CGP)[2],它們分別識別與單個藥物有關聯(lián)的基因組分子特征.然而，在癌癥的治療中，往往都是多個藥物聯(lián)合使用，因為一個藥物對應的靶標基因不止一個,同時,一個靶標基因所靶向的藥物也不止一個[3].通過探索多個藥物與多個基因之間的聯(lián)系,即模塊化分析,可以更直觀地發(fā)現(xiàn)高維數(shù)據(jù)之間隱含的生物信息.胡尊勝等[4]基于復雜網(wǎng)絡,研究蛋白質界面網(wǎng)絡中的模體和模塊,有利于研究蛋白質界面網(wǎng)絡的形成機制;Zhang等[5]提出稀疏網(wǎng)絡正則約束的非負矩陣分解模型，對microRNA-基因共模塊進行識別;Chen等[6]提出帶有網(wǎng)絡正則約束的稀疏偏最小二乘模型(SNPLS)，對基因-藥物共模塊進行識別，通過對各模塊的分析,有利于發(fā)現(xiàn)不同藥物之間相互作用的分子機理.受SNPLS模型的啟發(fā),本文提出伴有基因和藥物關聯(lián)網(wǎng)絡正則約束的稀疏偏最小二乘(SGDPLS)算法,并將其應用于癌癥藥物敏感性基因組學數(shù)據(jù)庫中.

1 材料和方法

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

從癌癥藥物敏感性基因組學數(shù)據(jù)庫(genomics of drug sensitivity in cancer，GDSC)中下載了細胞系的基因表達數(shù)據(jù)X1∈R641×13 321和細胞系藥物響應數(shù)據(jù)X2∈R641×95.從Pathway Commons數(shù)據(jù)庫(http:∥www.pathwaycommons.org/)中下載基因關聯(lián)網(wǎng)絡A,該網(wǎng)絡對應的鄰接矩陣為A.從NCBI PubChem數(shù)據(jù)庫中下載藥物的化學結構SDF文件，使用OpenBabelGUI工具箱讀取藥物SDF文件，并將藥物特征轉化為分子指紋,通過計算Jaccard相關系數(shù)[7],得到藥物相似性矩陣D.將D中小于0.8的位置賦值為零,其余位置賦值為1，得到藥物關聯(lián)網(wǎng)絡B，該網(wǎng)絡對應的鄰接矩陣為B.為了保證數(shù)據(jù)都處于同一數(shù)量級,且避免離群數(shù)據(jù)的出現(xiàn)，對模型中的初始數(shù)據(jù)X1和X2作了Z-score處理.

1.2 伴有基因和藥物關聯(lián)網(wǎng)絡正則約束的稀疏偏最小二乘算法

在SNPLS算法的基礎上，增加了藥物關聯(lián)網(wǎng)絡信息,構建SGDPLS算法,其流程圖如圖1所示.圖1中：H是藥物Docetaxel對應的化學結構；目標函數(shù)(1)中的L和C是由鄰接矩陣A和B經(jīng)過拉普拉斯矩陣變換得到的[8].SGDPLS算法對權向量g和d進行稀疏約束,可以使模型在迭代優(yōu)化過程中不陷入過擬合狀態(tài),增加識別模塊的可解釋性.

圖1 SGDPLS算法流程圖Fig.1 SGDPLS algorithm flow chart

1.3 SGDPLS模型的求解

圖1中的目標函數(shù)(1)是一個非凸函數(shù)，利用傳統(tǒng)的優(yōu)化算法很難得到模型的全局最優(yōu)解.Chen等[6]采用交替坐標下降法，通過交替更新g和d,求得模型的局部極大值.文中將該方法用于SGDPLS模型的求解,SGDPLS模型的求解過程有以下3個步驟.

步驟2交替更新g和d.

固定變量g,更新變量d,則

v=X2d.

固定變量d，更新變量g,則

u=X1g.

步驟3重復步驟2直到u收斂，算法終止.

對目標函數(shù)(1)做上述迭代運算,得到權向量g和d識別共模塊之前，先對g和d作Z-score標準化處理.通過設定閾值T，篩選共模塊的成員.對于基因(藥物)模塊,選取g(d)中比閾值T1(T2)大的位置所對應的基因(藥物)作為基因(藥物)模塊的成員.對于樣本模塊,首先計算u=X1g和v=X2d；然后,對u+v作標準化處理,取向量中比閾值T3大的位置所對應的樣本作為樣本模塊的成員.在篩選樣本、基因和藥物共模塊時，設置的閾值分別為T1=3,T2=3,T3=2.閾值過大(小)，會使識別出的共模塊過小(大).模塊過大,會導致模塊中包含的信息太多,有用的信息不易被挖掘出來；模塊過小,則包含的信息過少,不具有生物可解釋性[4].

2 結果分析

將SGDPLS模型用于GDSC數(shù)據(jù)庫,通過MATLAB編程[10]調參,識別出20個不同的基因-藥物共模塊.SGDPLS模型識別共模塊分布圖,如圖2所示.由圖2可知：有17個樣本模塊其樣本個數(shù)分布在24～34范圍內(圖2(a))；有13個基因模塊其基因個數(shù)分布在111～151范圍內(圖2(b))；12個藥物模塊的藥物個數(shù)為2(圖2(c)).

(a) 樣本個數(shù) (b) 基因個數(shù) (c) 藥物個數(shù)圖2 SGDPLS模型識別共模塊分布圖Fig.2 Distribution map of SGDPLS model recognition common modules

為了解釋文中識別出的共模塊具有生物意義，利用R語言中的clusterProfiler包對識別出的基因模塊進行GO生物功能項和KEGG通路富集分析.對于20個基因模塊,有14個基因模塊(70%)至少富集一種GO生物過程；有12個基因模塊(60%)至少富集一種KEGG通路.例如,第2個共模塊,樣本模塊包含28個樣本,這些樣本主要富集的癌癥類型是小細胞肺癌和小細胞癌.第2個共模塊生物項的富集結果,如圖3所示.圖3中：p表示富集分數(shù)經(jīng)過Benjamini校正后的值.

基因模塊包含107個基因,這些基因主要富集的GO生物過程是DNA構象變化、核染色體分離、姐妹染色體分離(圖3(a))，這些過程都發(fā)生在細胞增殖過程中.因為姐妹染色體分離到子細胞的過程是不可逆的,檢測點的缺失會導致染色體不穩(wěn)定,所以,在細胞周期中，有絲分裂過程最易發(fā)生癌變[11].現(xiàn)代醫(yī)學利用相關藥物制止細胞中紡錘體的出現(xiàn),從而抑制細胞有絲分裂的進行，使細胞分裂停留在G0階段,利用該技術可以有效地遏制癌細胞的惡性增殖和擴散.另外,基因模塊主要富集的KEGG通路為上皮細胞的細菌侵入、小細胞肺癌、細胞周期、焦點粘連、同源重組等細胞生化過程(圖3(b)),這些生化過程對癌癥的產(chǎn)生和轉移起著很重要的作用.藥物模塊包含2個藥物，即Docetaxel和RDEA119.Riichiroh等[12]指出，藥物Docetaxel已在臨床上被證實可以用于小細胞肺癌的治療.藥物RDEA119[13]是一種MEK抑制劑,MEK是一種磷酸化絲裂原活化蛋白激酶(MAPK)的激酶,MEK酶的靶標ERK具有高選擇性,且可以驅動細胞增殖,從而實現(xiàn)抑制腫瘤細胞增殖并誘導細胞凋亡.研究表明,該藥物對組織的選擇性,可以降低其對中樞神經(jīng)的毒副作用,且口服性的MEK抑制劑RDEA119已在臨床上被開發(fā)應用,有效地推動了藥代動力學的發(fā)展[14].

(a) BP (b) KEGG圖3 第2個共模塊生物項的富集結果Fig.3 Enrichment results of the second common module biological item

3 SGDPLS算法與SNPLS算法的比較

為了展示藥物關聯(lián)網(wǎng)絡在SGDPLS算法中所起的重要作用，將SNPLS算法也用于同樣的數(shù)據(jù)集上，分別計算2種算法識別出來的20個基因-藥物共模塊之間的皮爾遜相關系數(shù)，結果如圖4所示.

(a) 相關性分布箱線圖 (b) 單個模塊相關性柱狀圖 圖4 兩種算法共模塊皮爾遜相關性比較Fig.4 Comparison of co-module pearson correlation between two algorithms

由圖4可知：SGDPLS算法識別的大部分共模塊(65%)比SNPLS算法具有更好的相關性，說明由藥物二維化學結構構造的藥物關聯(lián)網(wǎng)絡，在一定程度上，可以提高共模塊中基因模塊和藥物模塊數(shù)據(jù)模式的相似性，使共模塊更具生物意義.

4 結束語

提出了SGDPLS算法，并深入挖掘多個基因和多個藥物之間的對應關系，針對SGDPLS算法和SNPLS算法，獨立識別出的20個基因-藥物共模塊，從共模塊之間的皮爾遜相關系數(shù)和生物意義角度進行分析.結果表明:SGDPLS算法識別的共模塊中，有65%比SNPLS算法識別的共模塊具有更高的相關性.由此可見,藥物關聯(lián)網(wǎng)絡信息的加入增強了SGDPLS算法識別模塊的相關性,有助于發(fā)現(xiàn)潛在的藥物靶標.同時，SGDPLS算法有望應用于其他領域中,進行高維數(shù)據(jù)的降維,如航天遙感數(shù)據(jù)、金融市場交易數(shù)據(jù),用以發(fā)現(xiàn)數(shù)據(jù)的本質結構[15].