抗乳腺癌候選藥物活性值模型建立和分析

亓凱航，仲梁維

(上海理工大學(xué)，上海 200093)

0 引言

乳腺癌是世界上最常見，致死率較高的癌癥之一。目前，抗激素治療經(jīng)常用于ERα表達(dá)的乳腺癌患者，其通過調(diào)節(jié)雌激素受體活性來控制人體內(nèi)雌激素水平[1]。ERα被公認(rèn)為是治療乳腺癌的重要靶標(biāo)，能夠拮抗ERα活性的化合物是治療乳腺癌的候選藥物[2]。在藥物研發(fā)中，為了節(jié)約時(shí)間及經(jīng)費(fèi)成本，通常采用建立化合物活性預(yù)測模型的方法來篩選潛在活性化合物。本文采用遺傳算法優(yōu)化BP 神經(jīng)網(wǎng)絡(luò)結(jié)合斯皮爾曼等級相關(guān)系數(shù)建立化合物活性預(yù)測模型。該模型能較為準(zhǔn)確地預(yù)測化合物活性值。

1 樣本數(shù)據(jù)處理

生物活性值作為因變量，分子描述符信息作為自變量。分析數(shù)據(jù)可知，將所有操作變量輸入，可能會(huì)導(dǎo)致過擬合，影響結(jié)果。因此采用斯皮爾曼等級相關(guān)系數(shù)對分子描述符和生物活性值進(jìn)行相關(guān)性分析計(jì)算，求得分子描述符和生物活性值之間的斯皮爾曼等級相關(guān)系數(shù)，將系數(shù)由大到小排序并返回索引值，得到前20個(gè)對生物活性最具有顯著影響的分子描述符，詳見表1。

表1 前20個(gè)對生物活性最具顯著影響的分子描述符

2 生物活性值預(yù)測模型

2.1 模型選擇

用前20 個(gè)對生物活性最具顯著影響的分子描述符及pIC50值構(gòu)建BP 神經(jīng)網(wǎng)絡(luò)模型。BP 神經(jīng)網(wǎng)絡(luò)是一種單向傳播多層的前饋神經(jīng)網(wǎng)絡(luò)，過程主要可分為兩個(gè)階段，第一階段是信號的前向傳播，再從輸入層經(jīng)過隱含層，最后到達(dá)輸出層；第二階段是誤差的反向傳播，從輸出層到隱含層，最后到輸入層，依次調(diào)節(jié)隱含層到輸出層的權(quán)重和偏置，輸入層到隱含層的權(quán)重和偏置[4]。

一個(gè)神經(jīng)網(wǎng)絡(luò)包括輸入層、隱含層和輸出層，如圖1所示。輸入層神經(jīng)元個(gè)數(shù)應(yīng)該與輸入數(shù)據(jù)的維數(shù)相同，輸出層神經(jīng)元個(gè)數(shù)與要擬合的數(shù)據(jù)維數(shù)相同，而隱含層神經(jīng)元的個(gè)數(shù)與層數(shù)按照設(shè)計(jì)者的目標(biāo)和需求決定[5]。深度學(xué)習(xí)隱含層的層數(shù)一般為三層，由可修正的權(quán)值互聯(lián)，一般使用的神經(jīng)網(wǎng)絡(luò)都是三層網(wǎng)絡(luò)。BP神經(jīng)網(wǎng)絡(luò)具有較強(qiáng)的非線性映射能力、高度自學(xué)習(xí)和自適應(yīng)的能力、一定的容錯(cuò)能力，但是BP 算法是一種基于梯度的最速下降法，因此存在不可避免的缺陷，學(xué)習(xí)過程收斂速度比較慢，網(wǎng)絡(luò)結(jié)構(gòu)不易確定，甚至網(wǎng)絡(luò)訓(xùn)練很容易陷入局部極小值[6]。故針對BP神經(jīng)網(wǎng)絡(luò)算法存在的缺點(diǎn)，提出利用遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)，避免陷入局部最優(yōu)解的方案，充分利用BP神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)能力和遺傳算法的全局搜索能力，二者相互補(bǔ)償。

圖1 神經(jīng)網(wǎng)絡(luò)

遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)包含種群初始化、適應(yīng)度函數(shù)、選擇操作、交叉操作和變異操作五個(gè)要素[7]。

⑴種群初始化

種群初始化運(yùn)用實(shí)數(shù)進(jìn)行編碼，將每個(gè)個(gè)體當(dāng)作一個(gè)實(shí)數(shù)串，包括輸入層與隱含層連接權(quán)值、隱含層閾值、隱含層與輸出層連接權(quán)值以及輸出層閾值四部分。種群個(gè)體都包括神經(jīng)網(wǎng)絡(luò)的全部權(quán)值和閾值。在已知網(wǎng)絡(luò)節(jié)點(diǎn)結(jié)構(gòu)的情況下，構(gòu)建結(jié)構(gòu)、權(quán)值、閾值確定的BP神經(jīng)網(wǎng)絡(luò)[7]。

⑵適應(yīng)度函數(shù)

將個(gè)體得到BP神經(jīng)網(wǎng)絡(luò)的初始權(quán)值和閾值，用訓(xùn)練集訓(xùn)練BP神經(jīng)網(wǎng)絡(luò)后預(yù)測系統(tǒng)輸出，把預(yù)測輸出和期望輸出間的誤差絕對值和作為個(gè)體適應(yīng)度值F，計(jì)算公式為：

其中，n為網(wǎng)絡(luò)輸出節(jié)點(diǎn)數(shù)；yi為BP 神經(jīng)網(wǎng)絡(luò)第i個(gè)節(jié)點(diǎn)的期望輸出；oi為第i個(gè)節(jié)點(diǎn)的預(yù)測輸出；k為系數(shù)[7]。

⑶選擇操作

常見的遺傳算法的選擇操作，大多數(shù)使用輪盤賭法、錦標(biāo)賽法等方法，本文考慮到適應(yīng)度比例，因此選用輪盤賭法，得到個(gè)體i的選擇概率pi為：

其中，fi為個(gè)體i的適應(yīng)度值，為了使適應(yīng)度值盡量小，會(huì)采用在個(gè)體選擇前對適應(yīng)度值求倒數(shù)的形式；k為系數(shù)；N為種群個(gè)體數(shù)目[7]。

⑷交叉操作

因?yàn)榉N群個(gè)體編碼采用的是實(shí)數(shù)編碼，所以交叉操作選用實(shí)數(shù)交叉法，第k個(gè)染色體ak和第l個(gè)染色體al在j位的交叉操作方法如下[8]：

其中，b是[ 0,1 ]之間的隨機(jī)數(shù)。

⑸變異操作

選取第i個(gè)個(gè)體的第j個(gè)基因aij進(jìn)行變異，變異操作方法如下：

其中，amax為基因aij的上界；amin為基因aij的下界；r2為一個(gè)隨機(jī)數(shù)；g是當(dāng)前迭代次數(shù)；Gmax為最大進(jìn)化次數(shù)；r為[ 0,1 ]間的隨機(jī)數(shù)[8]。

利用遺傳算法來優(yōu)化BP 神經(jīng)網(wǎng)絡(luò)得到最優(yōu)個(gè)體的權(quán)值和閾值，種群中的每個(gè)個(gè)體都包含了一個(gè)網(wǎng)絡(luò)所有的權(quán)值和閾值，個(gè)體利用適應(yīng)度函數(shù)計(jì)算個(gè)體適應(yīng)度值，遺傳算法使用選擇、交叉和變異操作找到最優(yōu)適應(yīng)度值對應(yīng)的個(gè)體[9]。BP神經(jīng)網(wǎng)絡(luò)預(yù)測就是使用遺傳算法得到最優(yōu)個(gè)體對網(wǎng)絡(luò)初始權(quán)值和閾值賦值，網(wǎng)絡(luò)經(jīng)過訓(xùn)練后對函數(shù)輸出預(yù)測。遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)流程如圖2所示。

圖2 遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)流程

2.2 模型建立

建立BP神經(jīng)網(wǎng)絡(luò)模型，步驟如下。

⑴數(shù)據(jù)輸入及處理

針對前20 個(gè)對生物活性最具顯著影響的分子描述符及生物活性值，在1974 個(gè)化合物中預(yù)留出474 個(gè)化合物的數(shù)據(jù)進(jìn)行之后的BP 神經(jīng)網(wǎng)絡(luò)系統(tǒng)的測試和檢驗(yàn)。

⑵網(wǎng)絡(luò)訓(xùn)練

用歸一化的1500組數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)訓(xùn)練，構(gòu)建該網(wǎng)絡(luò)結(jié)構(gòu)模型需要20個(gè)輸入?yún)?shù)、1個(gè)輸出參數(shù)，所以設(shè)置的BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)為20-＞41-＞1，即輸入層有20個(gè)節(jié)點(diǎn)，隱含層有41 個(gè)節(jié)點(diǎn)，輸出層有1 個(gè)節(jié)點(diǎn)，共有20 × 41+1=821 個(gè)權(quán)值，41+1=42 個(gè)閾值，所以遺傳算法個(gè)體編碼長度為821+42=863，均方誤差為0.065；學(xué)習(xí)率為0.1；訓(xùn)練結(jié)果與實(shí)際結(jié)果對比如圖3。

圖3 GA優(yōu)化BP神經(jīng)網(wǎng)絡(luò)輸出

⑶BP神經(jīng)網(wǎng)絡(luò)的優(yōu)化

本文采用遺傳算法對BP 神經(jīng)網(wǎng)絡(luò)的初始權(quán)值和閾值進(jìn)行優(yōu)化尋求最優(yōu)個(gè)體，在一定范圍內(nèi)對遺傳算法的參數(shù)設(shè)定，進(jìn)化代數(shù)為50，種群規(guī)模是50，交叉概率選擇為0.4，變異概率為0.3；適應(yīng)度曲線和預(yù)測誤差如圖4所示。

圖4 適應(yīng)度曲線和預(yù)測誤差

⑷網(wǎng)絡(luò)性能評估

均方誤差作為評價(jià)點(diǎn)估計(jì)最經(jīng)典的標(biāo)準(zhǔn)，均方誤差越小，越能說明點(diǎn)估計(jì)精確性。所謂均方誤差，就是預(yù)測數(shù)據(jù)和原始數(shù)據(jù)對應(yīng)點(diǎn)誤差的平方和的均值，計(jì)算公式為：

其中，n為樣本個(gè)數(shù)，SSE為和方差[8]。

優(yōu)化前后BP 神經(jīng)網(wǎng)絡(luò)誤差對比如圖5。首先，將未優(yōu)化前的BP 神經(jīng)網(wǎng)絡(luò)與優(yōu)化后的神經(jīng)網(wǎng)絡(luò)性能評估對比，對二者的最小均方誤差值作比較，優(yōu)化前的MSE值為0.22097，優(yōu)化后的MSE值為0.06923，經(jīng)過遺傳算法優(yōu)化后的BP 神經(jīng)網(wǎng)絡(luò)的MSE 值更小。其次，觀察預(yù)測輸出與期望輸出對比圖，遺傳算法優(yōu)化的輸出與期望輸出比較接近，擬合效果比較好。通過預(yù)測誤差得知，當(dāng)樣本數(shù)在100左右時(shí)，預(yù)測誤差最小。因此，經(jīng)過遺傳算法優(yōu)化后的BP神經(jīng)網(wǎng)絡(luò)對生物活性值的預(yù)測結(jié)果更加準(zhǔn)確。

圖5 優(yōu)化前后BP神經(jīng)網(wǎng)絡(luò)誤差對比

⑸神經(jīng)網(wǎng)絡(luò)預(yù)測

通過上述訓(xùn)練好的模型，將數(shù)據(jù)導(dǎo)入，得到50 個(gè)化合物的生物活性值，用公式pIC50=-lg(IC50× 10-9)將生物活性值(pIC50)轉(zhuǎn)化為IC50，最后得到生物活性值，詳見表2。

表2 生物活性預(yù)測值

3 結(jié)論

本文針對BP 神經(jīng)網(wǎng)絡(luò)預(yù)測化合物生物活性值存在學(xué)習(xí)過程收斂速度比較慢，容易陷入局部最優(yōu)解的缺陷，提出了基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型，可以有效降低均方誤差，遺傳算法優(yōu)化的BP 神經(jīng)網(wǎng)絡(luò)能夠較為準(zhǔn)確地預(yù)測化合物的生物活性值，預(yù)測結(jié)果符合預(yù)期值。