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

亓凱航，仲梁維

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

0 引言

乳腺癌是世界上最常見，致死率較高的癌癥之一。目前，抗激素治療經常用于ERα表達的乳腺癌患者，其通過調節雌激素受體活性來控制人體內雌激素水平[1]。ERα被公認為是治療乳腺癌的重要靶標，能夠拮抗ERα活性的化合物是治療乳腺癌的候選藥物[2]。在藥物研發中，為了節約時間及經費成本，通常采用建立化合物活性預測模型的方法來篩選潛在活性化合物。本文采用遺傳算法優化BP 神經網絡結合斯皮爾曼等級相關系數建立化合物活性預測模型。該模型能較為準確地預測化合物活性值。

1 樣本數據處理

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

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

2 生物活性值預測模型

2.1 模型選擇

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

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

圖1 神經網絡

遺傳算法優化的BP神經網絡包含種群初始化、適應度函數、選擇操作、交叉操作和變異操作五個要素[7]。

⑴種群初始化

種群初始化運用實數進行編碼，將每個個體當作一個實數串，包括輸入層與隱含層連接權值、隱含層閾值、隱含層與輸出層連接權值以及輸出層閾值四部分。種群個體都包括神經網絡的全部權值和閾值。在已知網絡節點結構的情況下，構建結構、權值、閾值確定的BP神經網絡[7]。

⑵適應度函數

將個體得到BP神經網絡的初始權值和閾值，用訓練集訓練BP神經網絡后預測系統輸出，把預測輸出和期望輸出間的誤差絕對值和作為個體適應度值F，計算公式為：

其中，n為網絡輸出節點數；yi為BP 神經網絡第i個節點的期望輸出；oi為第i個節點的預測輸出；k為系數[7]。

⑶選擇操作

常見的遺傳算法的選擇操作，大多數使用輪盤賭法、錦標賽法等方法，本文考慮到適應度比例，因此選用輪盤賭法，得到個體i的選擇概率pi為：

其中，fi為個體i的適應度值，為了使適應度值盡量小，會采用在個體選擇前對適應度值求倒數的形式；k為系數；N為種群個體數目[7]。

⑷交叉操作

因為種群個體編碼采用的是實數編碼，所以交叉操作選用實數交叉法，第k個染色體ak和第l個染色體al在j位的交叉操作方法如下[8]：

其中，b是[ 0,1 ]之間的隨機數。

⑸變異操作

選取第i個個體的第j個基因aij進行變異，變異操作方法如下：

其中，amax為基因aij的上界；amin為基因aij的下界；r2為一個隨機數；g是當前迭代次數；Gmax為最大進化次數；r為[ 0,1 ]間的隨機數[8]。

利用遺傳算法來優化BP 神經網絡得到最優個體的權值和閾值，種群中的每個個體都包含了一個網絡所有的權值和閾值，個體利用適應度函數計算個體適應度值，遺傳算法使用選擇、交叉和變異操作找到最優適應度值對應的個體[9]。BP神經網絡預測就是使用遺傳算法得到最優個體對網絡初始權值和閾值賦值，網絡經過訓練后對函數輸出預測。遺傳算法優化BP神經網絡流程如圖2所示。

圖2 遺傳算法優化BP神經網絡流程

2.2 模型建立

建立BP神經網絡模型，步驟如下。

⑴數據輸入及處理

針對前20 個對生物活性最具顯著影響的分子描述符及生物活性值，在1974 個化合物中預留出474 個化合物的數據進行之后的BP 神經網絡系統的測試和檢驗。

⑵網絡訓練

用歸一化的1500組數據進行網絡訓練，構建該網絡結構模型需要20個輸入參數、1個輸出參數，所以設置的BP 神經網絡結構為20-＞41-＞1，即輸入層有20個節點，隱含層有41 個節點，輸出層有1 個節點，共有20 × 41+1=821 個權值，41+1=42 個閾值，所以遺傳算法個體編碼長度為821+42=863，均方誤差為0.065；學習率為0.1；訓練結果與實際結果對比如圖3。

圖3 GA優化BP神經網絡輸出

⑶BP神經網絡的優化

本文采用遺傳算法對BP 神經網絡的初始權值和閾值進行優化尋求最優個體，在一定范圍內對遺傳算法的參數設定，進化代數為50，種群規模是50，交叉概率選擇為0.4，變異概率為0.3；適應度曲線和預測誤差如圖4所示。

圖4 適應度曲線和預測誤差

⑷網絡性能評估

均方誤差作為評價點估計最經典的標準，均方誤差越小，越能說明點估計精確性。所謂均方誤差，就是預測數據和原始數據對應點誤差的平方和的均值，計算公式為：

其中，n為樣本個數，SSE為和方差[8]。

優化前后BP 神經網絡誤差對比如圖5。首先，將未優化前的BP 神經網絡與優化后的神經網絡性能評估對比，對二者的最小均方誤差值作比較，優化前的MSE值為0.22097，優化后的MSE值為0.06923，經過遺傳算法優化后的BP 神經網絡的MSE 值更小。其次，觀察預測輸出與期望輸出對比圖，遺傳算法優化的輸出與期望輸出比較接近，擬合效果比較好。通過預測誤差得知，當樣本數在100左右時，預測誤差最小。因此，經過遺傳算法優化后的BP神經網絡對生物活性值的預測結果更加準確。

圖5 優化前后BP神經網絡誤差對比

⑸神經網絡預測

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

表2 生物活性預測值

3 結論

本文針對BP 神經網絡預測化合物生物活性值存在學習過程收斂速度比較慢，容易陷入局部最優解的缺陷，提出了基于遺傳算法優化的BP神經網絡預測模型，可以有效降低均方誤差，遺傳算法優化的BP 神經網絡能夠較為準確地預測化合物的生物活性值，預測結果符合預期值。