,

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

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

黃莉,賀平安

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

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

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

0 引 言

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

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

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

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

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

1 數(shù) 據(jù)

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

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

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

表1 抗癌藥物反應敏感性數(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ù)預處理

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

第一步，標準化(normalizing)，對第二類和第三類數(shù)據(jù)進行預處理，具體操作如下，

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

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

2 方法和結果

2.1 低秩矩陣填充模型

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

標準的矩陣填充模型如下：

(1)

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

(2)

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

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

(3)

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

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重復步驟2和步驟3，當該算法收斂(迭代誤差小于10-8)或者達到某一給定的迭代次數(shù)(本文設置2000次)停止。

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

(4)

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

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

(5)

(6)

不妨令

則有

=2λ2tr(XTAX)，

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

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Σ.

結合Cai等[12]模型前兩項U和V的梯度，可以得出整個模型中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 算法評估

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

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

2.4 交叉驗證

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

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

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

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

表2 r和λ1取不同值時的RMSE

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

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

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

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

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

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

3 結 論

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

[1] 吳冠青,孫燕.惡性腫瘤的個性化治療[J].癌癥進展,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] 申卯興,鄭武團.矩陣的奇異值分解的應用[J].大學數(shù)學,1996(3):56-61.

[17] 遲彬,葉慶凱.用奇異值分解方法計算具有重特征值矩陣的特征矢量[J].應用數(shù)學和力學,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折交叉驗證法的支持向量選取算法[J].華中師范大學學報(自科版),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)絡出版日期： 2017-05-24

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

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

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

O29

A

1673- 3851 (2017) 06- 0881- 07

(責任編輯：廖乾生)