矩陣填充算法在抗癌藥物敏感性研究中的運(yùn)用

,

(浙江理工大學(xué)理學(xué)院，杭州 310018)

矩陣填充算法在抗癌藥物敏感性研究中的運(yùn)用

黃莉,賀平安

(浙江理工大學(xué)理學(xué)院，杭州 310018)

從不完整的數(shù)據(jù)推斷完整有效的數(shù)據(jù)，繼而對原始數(shù)據(jù)給出可靠的分析是一個(gè)重要的數(shù)學(xué)問題。根據(jù)低秩矩陣填充算法，提出一種融合癌細(xì)胞系基因表達(dá)數(shù)據(jù)相似性信息的低秩矩陣填充算法。應(yīng)用該算法對癌細(xì)胞系與抗癌藥物反應(yīng)的敏感性缺失數(shù)據(jù)進(jìn)行恢復(fù)，并對相對反應(yīng)低的數(shù)值進(jìn)行重評估。利用均方根誤差和10倍交叉驗(yàn)證法評估該算法，結(jié)果顯示該算法比已有算法的均方根誤差減少22.7%，說明該算法具有很好的數(shù)據(jù)恢復(fù)效果。

抗癌藥物敏感性；低秩矩陣填充；癌細(xì)胞系；10倍交叉驗(yàn)證；均方根誤差

0 引 言

腫瘤是在體內(nèi)外各種因素的作用下由一系列基因連續(xù)突變導(dǎo)致細(xì)胞生長失去控制所致，因而對于每個(gè)腫瘤患者，即使是同一種腫瘤，其致病因素和體內(nèi)突變的基因可能都不相同。每一個(gè)患者的腫瘤都有獨(dú)特的生物特征，即腫瘤的異質(zhì)性[1]。腫瘤異質(zhì)性要求對不同的腫瘤患者要不同對待，即腫瘤的個(gè)性化醫(yī)療[2]，目前，以患者個(gè)人遺傳信息為基礎(chǔ)的個(gè)性化醫(yī)療已成為醫(yī)學(xué)領(lǐng)域主要研究方向之一。

在臨床試驗(yàn)中，為了研究和獲得針對于每個(gè)腫瘤患者有效個(gè)性化醫(yī)療，通常采用漫長而昂貴的藥物開發(fā)與實(shí)驗(yàn)驗(yàn)證來評估藥物的療效和毒性，但資源的稀缺限制該方案的實(shí)際應(yīng)用。而解決該問題可能性方案是將腫瘤細(xì)胞在體外進(jìn)行二維或三維培養(yǎng)，再直接測定病人腫瘤細(xì)胞的藥物敏感性[3]，或者采用癌細(xì)胞體內(nèi)培養(yǎng)模型，如小鼠異種移植模型和轉(zhuǎn)基因小鼠模型[4]。該方案可以捕獲大部分病人的腫瘤相關(guān)生物學(xué)特性；然而，這種方案代價(jià)大且耗時(shí)，缺乏可擴(kuò)展性，測試藥物數(shù)量最多只能到幾百種。

在過去幾十年中，隨著高通量技術(shù)發(fā)展，研究人員提出一種替代方案，即利用癌細(xì)胞系大面板建立藥物反應(yīng)的基因組預(yù)測因子[5-9]。目前預(yù)測藥物敏感性的方法主要根據(jù)已知癌細(xì)胞株對不同藥物的反應(yīng)特性和功能基因組特征,對藥物反應(yīng)數(shù)據(jù)進(jìn)行回歸或分類[10]。

矩陣填充算法是根據(jù)已有的數(shù)據(jù)對缺失數(shù)據(jù)進(jìn)行預(yù)測和恢復(fù)，Keshavan等[11]對該算法的原理和正確性作了全面闡述和證明，Cai等[12]在此基礎(chǔ)上提出改進(jìn)算法，并將該算法運(yùn)用到1968—2003年H3N2數(shù)據(jù)上，對血凝素抑制試驗(yàn)中缺失數(shù)據(jù)進(jìn)行預(yù)測。

本文根據(jù)低秩矩陣填充算法，提出一種新型預(yù)測藥物敏感性的方法，將癌細(xì)胞系基因表達(dá)數(shù)據(jù)相似性信息融合進(jìn)已有低秩矩陣填充算法中。利用該算法對癌細(xì)胞系百科全書(cancer cell line encyclopedia，CCLE)中491種癌細(xì)胞系與24種抗癌藥物反應(yīng)的敏感性數(shù)據(jù)進(jìn)行缺失數(shù)據(jù)填充和低反應(yīng)數(shù)據(jù)重估，使得CCLE數(shù)據(jù)庫中抗癌藥物反應(yīng)的敏感性數(shù)據(jù)更加詳盡和完整；獲得模型中參數(shù)r、λ1、λ2的最優(yōu)取值，并通過與已有算法比較證明該算法的有效性。

1 數(shù) 據(jù)

1.1 抗癌藥物反應(yīng)敏感性數(shù)據(jù)

CCLE是由哈佛大學(xué)、麻省理工Broad研究院和諾華生物研究所等研究機(jī)構(gòu)開發(fā)的首個(gè)可免費(fèi)獲取的癌癥基因組數(shù)據(jù)大型公共資源平臺[13]。CCLE整合來自大約1036種人類癌癥細(xì)胞系的基因表達(dá)、染色體拷貝數(shù)等大規(guī)模數(shù)據(jù)，并且還包含504種細(xì)胞系與24種抗癌藥物反應(yīng)的敏感性數(shù)據(jù)，并覆蓋36種常見癌癥類型[14]。CCLE旨在為癌癥研究提供數(shù)據(jù)支持以獲得更多的發(fā)現(xiàn)，通過理解癌癥細(xì)胞系與藥物之間關(guān)系獲取抗癌藥物潛在的敏感性相關(guān)標(biāo)志物,最終為癌癥治療尋找合適藥物[15]。

本文從CCLE數(shù)據(jù)庫(http://www.broadinstitute.org/ccle)下載504種癌細(xì)胞系與24種藥物反應(yīng)的藥物敏感性數(shù)據(jù)及1036種癌細(xì)胞系的基因表達(dá)數(shù)據(jù)。并將504種癌細(xì)胞系與1036種癌細(xì)胞系進(jìn)行匹配，發(fā)現(xiàn)其中只有491種癌細(xì)胞系相同。因此，本文實(shí)際運(yùn)用數(shù)據(jù)為491種癌細(xì)胞系與24種藥物反應(yīng)的藥物敏感性數(shù)據(jù)以及這491種癌細(xì)胞系的基因表達(dá)數(shù)據(jù)。表1為部分抗癌藥物反應(yīng)敏感性數(shù)據(jù)。

表1 抗癌藥物反應(yīng)敏感性數(shù)據(jù)

注：表中“NA”表示敏感值缺失。

表1數(shù)據(jù)表明，部分藥物敏感性數(shù)據(jù)缺失，此外部分藥物敏感性數(shù)據(jù)值為0，因此，藥物敏感性數(shù)據(jù)大致可以被分為3種類型：第一類數(shù)據(jù)，大于0的敏感值；第二類數(shù)據(jù)，等于0的敏感值；第三類數(shù)據(jù)，敏感值缺失。

1.2 數(shù)據(jù)預(yù)處理

為了方便分析，本文將491種癌細(xì)胞系與24種藥物反應(yīng)的藥物敏感性數(shù)據(jù)簡稱CD-491數(shù)據(jù)集，其中的可觀測數(shù)據(jù)有11360個(gè)。假設(shè)mij表示第i種癌細(xì)胞系與第j種藥物反應(yīng)敏感值，需要對原始的藥物敏感性數(shù)據(jù)進(jìn)行預(yù)處理，預(yù)處理包括兩步：

第一步，標(biāo)準(zhǔn)化(normalizing)，對第二類和第三類數(shù)據(jù)進(jìn)行預(yù)處理，具體操作如下，

其中，min(mij)表示CD-491數(shù)據(jù)集中非零最小值。

第二步，修整(trimming)，為了避免奇異向量高度集中在高權(quán)重的列(或行)，需要將矩陣中一些觀察值隨機(jī)設(shè)為0[11]。對于CD-491數(shù)據(jù)集中每一行，如果觀察值個(gè)數(shù)大于|E|/m(|E|表示矩陣E中觀察值個(gè)數(shù)，m表示矩陣E的行數(shù))，就隨機(jī)設(shè)置一些觀察值為0。同樣地，對于每一列，如果觀察值個(gè)數(shù)大于|E|/n(n表示矩陣E的列數(shù))，隨機(jī)將一些觀察值設(shè)置為0。

2 方法和結(jié)果

2.1 低秩矩陣填充模型

矩陣填充就是對矩陣中缺失的數(shù)據(jù)進(jìn)行恢復(fù)和已有數(shù)據(jù)的矯正。假設(shè)M=(mij)m×n的矩陣和一個(gè)集合E，矩陣M的行表示是m種癌細(xì)胞系，列表示n種抗癌藥物，mij表示第i種癌細(xì)胞系與第j種藥物反應(yīng)的敏感值,反應(yīng)值可以分為3種類型，而E為第一類型數(shù)據(jù)與第二類型數(shù)據(jù)位置的集合，即(i,j)∈E?[m]×[n]。假設(shè)矩陣M是秩為r的低秩矩陣，且r?m,n，則由奇異值分解[16-17]可知：存在矩陣Um×r、矩陣Vn×r和矩陣Σr×r，使得M=Um×rΣr×r(Vn×r)T。

標(biāo)準(zhǔn)的矩陣填充模型如下：

(1)

Cai等[12]在構(gòu)建流感病毒血凝素抑制試驗(yàn)數(shù)據(jù)恢復(fù)模型時(shí)認(rèn)為，上述矩陣填充模型只是充分運(yùn)用第一類型數(shù)據(jù)信息，雖然第二類型數(shù)據(jù)比第一類數(shù)據(jù)信息少，但能提高矩陣填充后數(shù)據(jù)準(zhǔn)確性，于是，引入一個(gè)閾值θij，當(dāng)矩陣中數(shù)據(jù)屬于第二類型數(shù)據(jù)時(shí)，θij被設(shè)置成常數(shù)C；若該數(shù)據(jù)不屬于第二類型數(shù)據(jù)，θij=-∞。在模型(1)基礎(chǔ)上，提出一個(gè)修正的融合第二類數(shù)據(jù)信息的低秩矩陣填充模型：

(2)

但是，上述兩個(gè)模型都僅僅是從數(shù)據(jù)的本身出發(fā)，忽略數(shù)據(jù)中對象的本身特征信息學(xué)，例如蛋白質(zhì)序列相似性信息，基因表達(dá)信息等。本文在模型(2)基礎(chǔ)上提出一種融合數(shù)據(jù)對象的特征信息的低秩矩陣填充模型：

s.t.X=Um×rΣr×r(Vn×r)T

(3)

其中：sij表示癌細(xì)胞系i與癌細(xì)胞系j基因表達(dá)數(shù)據(jù)的Person相關(guān)性系數(shù)，Xi表示矩陣X的第i行。此模型是基于以下觀點(diǎn)提出的：如果兩個(gè)癌細(xì)胞系的基因表達(dá)數(shù)據(jù)相關(guān)性系數(shù)越高，那么它們的藥物反應(yīng)敏感性數(shù)據(jù)越相似。即當(dāng)這些癌細(xì)胞系被投射到一個(gè)幾何空間中，相關(guān)性系數(shù)越高的癌細(xì)胞系，它們的空間距離越近。

2.2 模型算法

為了求解模型(3)，本文提出了一種基于梯度下降的算法，算法迭代步驟如下:

步驟3通過梯度下降法更新U(k+1)和V(k+1):即U(k+1)=U(k)+t·▽U(k)，V(k+1)=V(k)+t·▽V(k)，▽U和▽V分別表示U和V的梯度；

步驟4重復(fù)步驟2和步驟3，當(dāng)該算法收斂(迭代誤差小于10-8)或者達(dá)到某一給定的迭代次數(shù)(本文設(shè)置2000次)停止。

上述模型算法中需要分別計(jì)算U和V的梯度。由于本文模型比Cai等[12]提出的模型(2)增加一項(xiàng):

(4)

因此，需要對Cai等[12]的算法加以改進(jìn)。為了求解模型中U和V的梯度，本文首先引入引理1。

引理1[18]設(shè)A=(aij)∈Rm×m和B=(bij)∈Rm×m都是常數(shù)矩陣，X=(xij)∈Rm×m是一個(gè)變量矩陣，則有

(5)

(6)

不妨令

則有

=2λ2tr(XTAX)，

根據(jù)跡的性質(zhì)，有

tr(XTAX)=tr((UΣVT)TAUΣVT)

=tr(VΣTUTAUΣVT)

=tr(UTAUΣVTVΣT)

=tr(VTVΣTUTAUΣ).

由式(5)可知：

=2λ2(AUΣVTVΣT+ATU(ΣVTVΣT)T)

=4λ2AUΣVTVΣT.

由式(6)可知：

=2λ2(VΣTUTAUΣ+V(ΣTUTAUΣ)T)

=4λ2VΣTUTAUΣ.

結(jié)合Cai等[12]模型前兩項(xiàng)U和V的梯度，可以得出整個(gè)模型中U和V的梯度分別是：

▽U=((UΣVT)E-ME)·IVΣT+UQU+λ1f(U,2e(QU1-I1)2(QU1-I1))+4λ2AUΣ(VTV)ΣT,

▽V=((UΣVT)E-ME)T·IUΣ+VQV+λ1f(V,2e(Qv1-I2)2(Qv1-I2))+4λ2VΣTUTAUΣ.

2.3 算法評估

評價(jià)一個(gè)模型的好壞，有許多參數(shù)標(biāo)準(zhǔn)，比如和方差(SSE)、均方根誤差(RMSE)、確定系數(shù)(R-square)等。本文利用均方根誤差來評估低秩矩陣填充模型的優(yōu)劣。一般來說，均方根誤差值越小，表示預(yù)測值越接近真實(shí)值[19]。

假設(shè)有兩個(gè)向量X和Y,且X=(x1,x2,…,xk)，Y=(y1,y2,…,yk)，xi、yi分別代表觀察值和相應(yīng)的估計(jì)值，那么將RMSE定義為：

2.4 交叉驗(yàn)證

在本文矩陣填充算法中，利用10倍交叉驗(yàn)證法，通過訓(xùn)練參數(shù)λ1，λ2，從而獲取最小的RMSE值。所謂的10倍交叉驗(yàn)證法，就是均分樣本數(shù)據(jù)為10組，選取其中的一組數(shù)據(jù)作為測試集，其余9組數(shù)據(jù)作為訓(xùn)練集，重復(fù)10次，將10次結(jié)果均值作為對算法精度估計(jì)的依據(jù)[20]。

本文將矩陣中已有的數(shù)據(jù)均分為10等份，每一次，選取其中的9份數(shù)據(jù)作為觀察值進(jìn)行矩陣的填充，然后計(jì)算填充完后的矩陣與觀察矩陣在剩下1份數(shù)據(jù)集上的RMSE，重復(fù)10次，將10次結(jié)果均值作為對該次算法精度的估計(jì)值，即為需要RMSE值。每次僅對測試集中第一類型數(shù)據(jù)計(jì)算RMSE，RMSE值也被稱作局部RMSE值。

2.5 CD-491數(shù)據(jù)集中缺失值的填充

首先將模型(2)運(yùn)用到CD-491數(shù)據(jù)集中，用以預(yù)測數(shù)據(jù)集中的缺失值和重估值為0的數(shù)據(jù)。表2給出10倍交叉驗(yàn)證時(shí)，r，λ1取不同值情況下的RMSE。從表2可以看出，當(dāng)參數(shù)r=3，λ1=1×10-3時(shí)，模型(2)算法的RMSE最小，其值為0.8654。

表2 r和λ1取不同值時(shí)的RMSE

進(jìn)一步地，將模型(3)運(yùn)用到CD-491數(shù)據(jù)集中，為了得到完整矩陣，需要對3個(gè)未知參數(shù)進(jìn)行訓(xùn)練，除了矩陣X的秩r，模型中第二項(xiàng)系數(shù)λ1，還有第三項(xiàng)系數(shù)λ2，通過10倍交叉驗(yàn)證，分別獲取參數(shù)r、λ1、λ2最優(yōu)取值r=3，λ1=1×10-4，λ2=1×10-9，此時(shí)RMSE值達(dá)到最小0.6688，表3-6展示r、λ1、λ2取不同值情況下的RMSE。

表3 r=3時(shí)，λ1和λ2取不同值時(shí)的RMSE

表4 r=4時(shí)，λ1和λ2取不同值時(shí)的RMSE

表5 r=5時(shí)，λ1和λ2取不同值時(shí)的RMSE

表6 r=6時(shí)，λ1和λ2取不同值時(shí)的RMSE

結(jié)合這3個(gè)參數(shù)的最優(yōu)取值，利用模型(3)，得到一個(gè)比較完整的CD-491數(shù)據(jù)集，原先由于受到試驗(yàn)條件的限制和干擾而無法直接獲取的數(shù)據(jù)通過本文模型可以進(jìn)行預(yù)測。此外，通過比較可以發(fā)現(xiàn)，模型(4)由于增加第三項(xiàng)，均方根誤差由0.8654降到0.6688，即均方根誤差減少22.7%，這說明融合癌細(xì)胞株的基因表達(dá)數(shù)據(jù)有利于矩陣填充。

3 結(jié) 論

在抗癌藥物敏感性數(shù)據(jù)獲取實(shí)驗(yàn)中，實(shí)驗(yàn)條件限制和外界因素的干擾會導(dǎo)致部分?jǐn)?shù)據(jù)不準(zhǔn)或部分?jǐn)?shù)據(jù)缺失，利用數(shù)學(xué)方法從不完整數(shù)據(jù)中推斷出可靠與完整抗癌藥物敏感性數(shù)據(jù)是一個(gè)非常有意義的研究。本文根據(jù)低秩矩陣填充算法，提出一種合理、更具有生物意義低秩矩陣填充模型，該模型不僅僅從實(shí)驗(yàn)數(shù)據(jù)出發(fā)，而且充分利用數(shù)據(jù)對象本身的生物學(xué)特征，通過將癌細(xì)胞系基因表達(dá)數(shù)據(jù)融合進(jìn)低秩矩陣填充模型，提高矩陣填充的準(zhǔn)確性。

[1] 吳冠青,孫燕.惡性腫瘤的個(gè)性化治療[J].癌癥進(jìn)展,2008,6(6):562-578.

[2] CARNEY K. Personalized medicine[J]. Journal of the California Dental Association,2003,4(6):548-558.

[3] GRIFFITH L G, SWARTZ M A. Capturing complex 3D tissue physiology in vitro[J]. Nature Reviews Molecular Cell Biology,2006,7(3):211-224.

[4] RICHMOND A, SU Y. Mouse xenograft models vs GEM models for human cancer therapeutics[J]. Disease Models & Mechanisms,2008,1(2/3):78-82.

[5] GARNETT M J, EDELMAN E J, HEIDORN S J, et al. Systematic identification of genomic markers of drug sensitivity in cancer cells[J]. Nature,2012,483(7391):570-575.

[6] SHOEMAKER R H. The NCI60 human tumour cell line anticancer drug screen[J]. Nature Reviews Cancer,2006,6(10):813-823.

[7] HEISER L M, WANG N J, TALCOTT C L, et al. Integrated analysis of breast cancer cell lines reveals unique signaling pathways[J]. Genome Biology,2009,10(3):31.

[8] KUTALIK Z, BECKMANN J S, BERGMANN S. A modular approach for integrative analysis of large-scale gene-expression and drug-response data[J]. Nature Biotechnology,2008,26(5):531-539.

[9] KANDELA I, ZERVANTONAKIS I. Registered report: Discovery and preclinical validation of drug indications using compendia of public gene expression data[J]. Science Translational Medicine,2011,3(96):9677-9687.

[10] ZHANG N, WANG H, FANG Y, et al. Predicting anticancer drug responses using a dual-layer integrated cell line-drug network model[J]. Plos Computational Biology,2015,11(9):1-4.

[11] KESHAVAN R H, OH S, MONTANARI A. Matrix completion from a few entries[J]. IEEE Transactions on Information Theory,2009,56(6):2980- 2998.

[12] CAI Z, ZHANG T, WAN X F. A computational framework for influenza antigenic cartography[J]. Plos Computational Biology,2010,6(10):1922-1928.

[13] BARRETINA J, CAPONIGRO G, STRANSKY N, et al. The Cancer Cell Line Encyclopedia enables predictive modeling of anticancer drug sensitivity[J]. Nature,2012,483(7391):603-607.

[14] VENKATESAN K, STRANSKY N, MARGOLIN A, et al. Prediction of drug response using genomic signatures from the Cancer Cell Line Encyclopedia[J]. Clinical Cancer Research,2010,16:2-5.

[15] DONG Z, ZHANG N, LI C, et al. Anticancer drug sensitivity prediction in cell lines from baseline gene expression through recursive feature selection[J]. BMC Cancer,2015,15(1):1-12.

[16] 申卯興,鄭武團(tuán).矩陣的奇異值分解的應(yīng)用[J].大學(xué)數(shù)學(xué),1996(3):56-61.

[17] 遲彬,葉慶凱.用奇異值分解方法計(jì)算具有重特征值矩陣的特征矢量[J].應(yīng)用數(shù)學(xué)和力學(xué),2004,25(3):233-238.

[18] PETERSEN K, PEDERSEN M. The Matrix Cookbook[M]. Copenhagen: Technical University of Denmark,2012:12-14.

[19] BARNSTON A G. Correspondence among the correlation, RMSE, and Heidke forecast verification measures; refinement of the Heidke score[J]. Weather & Forecasting,2006,7(4):699-709.

[20] 牛曉太.基于KNN算法和10折交叉驗(yàn)證法的支持向量選取算法[J].華中師范大學(xué)學(xué)報(bào)(自科版),2014,48(3):335-338.

MatrixCompletionforPredictionoftheCancerousDrugSensitivity

HUANGLi,HEPingan

(School of Sciences, Zhejiang Sci-Tech University, Hangzhou 310018, China)

It is an important issue that how to use the incomplete data to obtain the complete and effective data, and then make the reliable analysis of the original data. In the work, based on a low-rank matrix completion, a novel low-rank matrix completion algorithm, which integrated the similarity information of gene expression data in cancer cell lines, was proposed to obtain reliable and complete anticancer drug sensitivity datasets. The model was applied to the observed datasets of cancer cell lines’ responses to chemical compounds, recovering the missing data and re-evaluating the low react value. Compared to previous methods, the root-mean-square error (RMSE) in the method is reduced by 22.7% in a 10-fold cross validation analysis, which demonstrated that the novel algorithm can improve matrix completion quality.

anticancer drug sensitivity; low-rank matrix completion; cancer cell lines; 10-fold cross validation analysis; root-mean-square error (RMSE)

10.3969/j.issn.1673-3851.2017.11.020

2017-03-03 網(wǎng)絡(luò)出版日期： 2017-05-24

國家自然科學(xué)基金項(xiàng)目(61170110,61272312)；浙江省自然科學(xué)基金項(xiàng)目(LY14F020049)

黃 莉(1990-)，女，湖北黃岡人，碩士研究生，主要從事生物信息學(xué)方面的研究。

賀平安，E-mail：pinganhe@zstu.edu.cn

O29

A

1673- 3851 (2017) 06- 0881- 07

(責(zé)任編輯：廖乾生)