遺傳算法優化的BP神經網絡拷貝數變異檢測

黃體浩，李俊青，趙海勇

聊城大學計算機學院，山東 聊城 252000

拷貝數變異（copy number variation，CNV）是指DNA拷貝數偏離正常拷貝數時基因組序列發生擴增或缺失，是基因組結構變異的一種重要形式。生物醫學研究估計，人類基因組中包含1 000 多個CNV 區域，約占基因序列的1%。CNV 與腫瘤疾病的產生和發展有著極其重要的聯系，也是引發眾多復雜腫瘤疾病的主要原因之一，因此發現原癌基因和抑癌基因尤為關鍵的一步就是準確檢測目標基因組的拷貝數變異。

利用下一代測序（next-generation sequencing，NGS）技術檢測拷貝數變異的策略主要分為四種：基于讀深（read depth，RD）、基于雙端比對、基于分裂讀段、基于序列組裝。其中最常用的拷貝數檢測方法是基于RD的方法，Xie等[1]通過利用病例樣本與對照樣本之間的讀段數目（read count，RC）比例建立統計模型，并根據比例差異識別CNV；Yoon等[2]采用基于事件測試的方法在個體基因組中檢測到的缺失和擴增跨多個基因組進行檢查，以確定個體之間的多態性；Ivakhno 等[3]采用離散小波變換對QC 過濾的RC 進行去噪，并使用隱馬爾可夫模型（hidden markov model，HMM）進行分段。早期基于RD的拷貝數檢測方法過度依賴RD信號，Abyzov等[4]采用均值漂移算法對單個樣本的RD信號進行處理，并以此統計模型建立檢測CNV；Xi等[5]將最小化貝葉斯信息準則算法擴展應用到多個樣本中，以檢測重復出現的CNV；Miller等[6]建立服從負二項分布的RC統計模型以減少均方根誤差；Boeva 等[7]對GC 含量進行建模，并將模型用于歸一化RD 信號；Gusnanto 等[8]通過損失函數糾正RC比例以消除GC含量偏差的影響。然而上述方法只能檢測出顯著性的CNV，且不能準確定位CNV 的斷點位置。反向傳播（back propagation，BP）神經網絡擁有強大的非線性映射能力，可用來解決拷貝數變異檢測中多個特征的聯合作用，然而BP 神經網絡在實際應用中仍存在諸多缺點，如局部最優、收斂速度慢、過擬合等。

基于上述研究，本文提出CNV-GABP（copy number variation detection of BP neural network based on genetic algorithm）檢測方法。方法分為四個部分：特征值計算與歸一化、BP神經網絡結構確定、遺傳算法優化和BP神經網絡預測。該方法的優點在于：（1）充分考慮了多個特征之間的聯合作用；（2）在BP 神經網絡基礎上，使用遺傳算法（genetic algorithm，GA）優化BP 神經網絡模型，提高模型的泛化能力；（3）優化后的BP 神經網絡用以解決多個特征之間的聯合作用，進一步提高了檢測準確率和精度，可有效識別特征不明顯的CNV。

1 相關工作

近年來，研究者綜合考慮RD信號和其他因素，提出了很多CNV檢測算法，Wang等[9-11]采用基于HMM的方法，分別利用期望最大化（expectation-maximization，EM）算法和置信傳播算法估計模型的最優參數；Onsongo等[12]采用隨機森林算法結合機器學習識別CNV；Johansson 等[13]選擇與分析樣本最具相似度的對照樣本，歸一化后分析兩個樣本間的平均覆蓋率差異；Dharanipragada等[14]使用覆蓋深度方法檢測CNV；Chen等[15]采用最大懲罰似然估計的方法識別CNV和檢測CNV邊界；Melcher等[16]和Jugas 等[17]分別采用局部多項式回歸擬合算法和中位數方法進行GC偏差矯正，然后利用全變分模型和循環二元分割算法進行分割讀取；Yuan等[18]采用基于統計混合模型的貝葉斯方法來推斷拷貝數的擴增和缺失；Xiao 等[19]采用基于高斯混合模型的聚類策略消除假陽性，并在聚類步驟中引入EM算法優化模型參數；Li等[20]建立高斯混合模型，并用先大窗后小窗的分割方式進行分割；Yuan等[21]采用孤立森林算法將RD信號生成孤立樹集合，并使用全變分模型進行去噪；Li等[22]在RD信號的基礎上引入了連續位置的相關性作為統計量以提高對CNV 的識別能力；Yuan 等[23]為每個分段分配一個離群點因子，采用局部異常因子算法識別CNV；Yuan等[24]利用腫瘤純度和RD信號建立非線性模型，采用窮舉搜索策略計算最優解；Zhao 等[25]利用BP 神經網絡構建多特征訓練算法以預測CNV；Yuan 等[26]利用潛變量模型重建觀察到的信號，并對其進行分層處理。上述方法不僅僅依賴RD信號，減少了由于讀段數目分布不均和定位模糊產生的影響，提高了靈敏度和特異性，但是仍然忽略了多個特征之間聯合作用的影響，并且由于測序誤差、GC 含量偏差等影響，一些特征較不明顯的CNV 無法被識別。

2 CNV-GABP檢測方法

本文提出的CNV-GABP方法是一種基于多個特征的檢測方法，用于從NGS 數據中準確檢測CNV。目前很多方法也逐漸考慮基因組位置之間的相關性，例如文獻[4]和文獻[18]都將基因組位置之間的相關性綜合到統計量中來提高RD信號的檢測效率，但是這些方法忽略了特征值之間的相互作用。CNV-GABP將GC含量、RD 信號、基因組相鄰位置之間的相關性以及比對質量這四個特征融入到基因組序列的分析中，并通過訓練BP 神經網絡使四個特征相互作用。為了克服BP 神經網絡全局搜索能力差、易陷入局部最優等缺點，本文使用遺傳算法優化BP 神經網絡模型，增強模型的泛化能力，提高方法的檢測性能。

2.1 算法流程

CNV-GABP 方法工作流程如圖1 所示。首先輸入一個參考基因組和一個測序樣本，然后對輸入數據進行預處理，最后經過四個步驟實現CNV 的預測：（1）定義與CNV相關的四個特征值，計算四個特征值并歸一化；（2）利用四個特征值構建基于遺傳算法優化的BP 神經網絡；（3）訓練BP神經網絡，其中標記為CNV的數據從仿真數據和真實數據中采樣獲取；（4）根據訓練完成的BP 神經網絡模型預測每個窗的CNV 狀態（正常、擴增或缺失）。

2.2 數據預處理與特征值量化

CNV-GABP使用比對算法BWA來處理格式為fasta的參考文件和格式為fastq 的序列樣本文件，以生成格式為bam 的比對文件。使用samtools 軟件從bam 文件中獲取RC 文件，RD 信號是根據RC 文件計算得來的，為方便計算RD 信號，CNV-GABP 將提取的RC 文件分成連續非重疊且大小相等的窗，可根據需求自由設置窗的大小，本文中窗的大小為1 000 bp，取每個窗內的平均RC值作為該窗的RD值。由于GC含量偏差對RD信號產生噪聲，通常如文獻[21]中的方法會預先對GC 含量偏差進行校正，以期降低GC含量偏差的影響，這種方法適用于一些特定的情況，但在大多數應用場景下并未取得十分顯著的效果。因此CNV-GABP并未同常規方法般預先對每個窗的GC 含量偏差進行校正，而是將GC含量作為該窗的特征值之一。基因組相鄰位置之間的相關性和比對質量在一定程度上反應出該區域的CNV狀態，因此CNV-GABP 也將這兩個因素作為特征值。其中基因組相鄰位置之間的相關性值參考文獻[25]，計算方式如式（1）：

其中Ci表示基因組第i個位置與相鄰位置的相關性值，Ri表示表示第i個窗的RD 值，n表示窗的總數，且j

每個窗的量化方式如式（2）：

其中，Gi表示GC 含量，Qi表示比對質量，Ci和Ri的定義由式（1）給出。

由于四個特征值的表示標準不同，為消除特征值之間的量綱影響，CNV-GABP 對特征值進行歸一化處理。以比對質量為例，歸一化方式如式（3）：

2.3 BP神經網絡構建

BP神經網絡是一種按照誤差逆向傳播算法訓練的多層前饋式神經網絡，以其優越的非線性映射能力和柔性的網絡結構被廣泛應用于分類問題，BP 神經網絡算法的學習過程包括信號的正向傳播與誤差的反向回傳。正向傳播時，數據從輸入層傳入，經隱含層處理，傳向輸出層，若輸出層輸出與期望輸出不一致，則將誤差作為調整信號反向回傳，調整神經元之間的權值矩陣，減小誤差。本文用歸一化之后的特征值組成的四元組Xi=(Ri,Gi,Qi,Ci)作為輸入構建BP 神經網絡。如圖2表示該網絡的拓撲結構。

圖2 BP神經網絡拓撲結構Fig.2 Topology of BP neural network

其中ω表示輸入層與隱含層、隱含層與輸出層之間的權重，θ表示偏差，X表示神經元輸入向量，即歸一化之后的特征值組成的四元組。除此之外，學習速率與收斂性密切相關，是神經網絡的重要參數之一，本文使用BP神經網絡的默認值0.1作為該網絡的學習速率。

2.4 基于遺傳算法的BP神經網絡

遺傳算法是一種啟發式隨機搜索算法，能夠在不需要誤差函數梯度信息的情況下學習近似最優解，具有良好的全局搜索能力和隱并行性，適合復雜的非線性問題求解。遺傳算法的主要步驟包括編碼、種群初始化、選擇、交叉和變異。使用遺傳算法優化BP 神經網絡主要是對神經網絡各項參數的調整，本文使用遺傳算法優化BP神經網絡的權值和閾值。

基于遺傳算法優化BP神經網絡算法的基本思想是將BP神經網絡算法與遺傳算法相結合，當BP算法訓練網絡的收斂速度較慢時，將BP 神經網絡每個隱含層節點的閾值和權重作為遺傳算法的輸入信息，并對它們進行編碼以生成染色體，然后使用遺傳算法的選擇算子、交叉算子和變異算子來生成新的后代作為BP算法的初始值，最后繼續使用BP 算法訓練網絡。基于遺傳算法優化的BP神經網絡算法流程如圖3所示。

圖3 基于遺傳算法的BP神經網絡程序流程圖Fig.3 Flow chart of BP neural network based on genetic algorithm

算法分為兩個部分：遺傳算法優化部分和BP 神經網絡部分。

（1）遺傳算法優化部分。首先對種群進行初始化，初始化種群由BP 神經網絡的初始權值和閾值組成，隨后執行選擇、交叉和變異操作，將種群中執行完操作的個體代入BP神經網絡進行迭代以尋找最優個體。具體內容如下：

步驟1 初始化種群。本文采用實數編碼，種群中的每個個體都是由輸入層與隱含層之間的權值、隱含層閾值、隱含層與輸出層之間的權值、輸出層閾值組成的實數串。

步驟2 計算適應度函數。本文將期望輸出與預測輸出之間的誤差絕對值作為個體的適應度值F，計算方式如下：

其中m表示神經網絡輸出節點總數，di表示神經網絡第i個節點的期望輸出，pi表示神經網絡第i個節點的預測輸出。

步驟3 選擇操作。本文采用輪盤賭法作為選擇策略，選擇適應度好的個體組成新的種群，選擇算子的計算方式如下：

其中Fi表示個體i的適應度值，N為種群中個體的總數。

步驟4 交叉操作。在種群中隨機選擇兩個個體按照一定概率進行交叉操作，兩個父代個體Pk和Pl在第j個點位的實數交叉方式如下：

其中α是介于（0，1）之間得隨機數。

步驟5 變異操作。在種群中隨機選擇一個個體按照一定概率進行變異操作，個體Pi在第j個點位變異方式如下：

其中Pmax表示個體的上界，Pmin表示個體的下界，v 表示當前迭代次數，Vmax表示最大進化次數，μ是介于（0，1）之間的隨機數。

步驟6 檢查新的個體是否符合最優個體標準，若不符合，則返回步驟2；若符合，則將最優個體分割為BP神經網絡的權值和閾值，調整網絡參數，重復學習訓練，直至達到所需精度或學習次數的上限為止。

（2）BP 神經網絡部分。將經過遺傳算法優化得到的初始權值和閾值等個體，代入BP神經網絡，輸入訓練數據，進行BP 神經網絡訓練。利用遺傳算法檢驗個體的適應度，無法達到適應度標準的個體進行選擇、交叉和變異操作，重新進行BP神經網絡訓練，直至個體滿足適應度標準或滿足精度要求為止。具體內容如下：

步驟1 初始化網絡。X為神經元輸入向量，n為輸入層節點數；Y為神經元輸出向量，m為輸出層節點數。

其中輸入向量根據式（1）進行量化，輸出向量分別代表拷貝數變異的3種狀態。

步驟2 計算隱含層輸出。隱含層節點數為p，輸入層與隱含層之間的權值為ωij，閾值為λa，隱含層輸出值Ha的計算方式如下：

其中f為激勵函數，計算方式如式（4）所示。

步驟3 計算輸出層輸出。隱含層與輸出層之間的權值為ωjk，閾值為λb，輸出層輸出值Hb的計算方式如下：

步驟4 計算誤差并驗證。誤差e表示為預測輸出Hb與期望輸出Y之差，計算方式如下：

假定e0為允許的誤差最大值，若max(em)>e0，則調整BP神經網絡參數，重新進行學習訓練，直至全部訓練數據滿足誤差標準。

步驟5 更新權值和閾值。根據上述步驟計算出的誤差e更新ωij和ωjk，計算方式如下：

其中η為學習速率，介于（0，1）之間。閾值修正如下：

BP 神經網絡算法的結果還取決于網絡的原始狀態，在局部搜索中使用BP 算法更為有效，GA 則擅長全局搜索。因此，本文使用GA 優化初始權值和閾值，并在解空間中找到一些更好的搜索空間，然后使用BP 算法在這些小的解空間中搜索最優解。通常使用混合算法的效果要優于僅使用BP神經網絡算法的效果。

3 實驗結果與分析

為防止細胞污染等外部因素造成的誤差，本文在不同的腫瘤純度和覆蓋深度中選擇樣本作為訓練數據。利用經過上述步驟訓練完成的BP 神經網絡模型，對測試數據集進行預測。對于每個窗，算法輸出3個從tansig函數映射的0 到1 之間的值，分別表示CNV 這3 種狀態的概率，并根據最大的概率確定該窗的CNV 狀態。例如，代表擴增的概率大于代表其他兩種CNV 狀態的概率，則該窗的CNV狀態為擴增。

本文中BP神經網絡結構為3層，利用式（4）給出的激活函數，輸入層節點數為4，隱含層節點數為5，輸出層節點數為3，共4×5+5×3=35個權值，5+3=8個閾值，因此需要優化43 個參數。交叉概率pc=0.2，變異概率pm=0.1，初始種群規模S=10。

3.1 仿真數據研究

為了驗證CNV-GABP 方法的合理性和可靠性，首先使用模擬數據進行實驗。本文從基因組變異模擬方法中選擇了一種綜合模擬方法IntSIM[27]來生成測序數據。腫瘤純度設置為0.2、0.3和0.4，覆蓋深度設置為4x和6x，以生成不同配置的數據集。從標準參考基因組中選擇一條染色體Chr21作為IntSIM輸入序列，并針對每種配置模擬50個樣本，共計300個樣本。將CNV-GABP與4 種現有方法（CNAnator[4]、FREEC[7]、MFCNV[25]、GROM_RD[28]）在靈敏度、精度和F1評分上進行比較，指標的計算方式參考文獻[16]。每個指標均取50 個樣本的測量均值，5種方法的性能比較結果如圖4所示。

圖4 仿真數據中5種方法的性能比較Fig.4 Performance comparison of five methods in simulation data

由圖4可以看出，腫瘤純度增加的同時，5種方法的檢測性能都在提高。例如腫瘤純度0.2、覆蓋深度均為4x 時GROM_RD 的靈敏度約為0.45，而腫瘤純度0.4 時約為0.56；5種方法F1 評分均值在腫瘤純度0.2、覆蓋深度均為4x 時約為0.59，而在腫瘤純度0.4 時約為0.71。就精度而言，在6類樣本中GROM_RD獲得2個最高值，CNVnator 獲得1 個最高值，CNV-GABP 獲得3 個最高值；就靈敏度和F1 評分而言，在6類樣本中CNV-GABP均獲得最高值。

檢測CNV邊界的準確性同樣是衡量算法可靠性的指標之一，邊界偏差越小，CNV 檢測方法性能越高。MFCNV在邊界偏差上取得了良好的效果，本文將CNVGABP與之比較，結果如圖5所示。

圖5 兩種方法的邊界偏差比較Fig.5 Boundary bias comparison of two methods

由圖5 可以看出，隨著腫瘤純度的增加，兩種方法的邊界偏差都隨之減小。實驗結果表明，CNV-GABP在六類樣本中的邊界偏差均小于MFCNV。因此，CNVGABP在模擬數據中的性能表現優于其他4種算法。

CNV-GABP 的性能優勢主要由三方面原因導致：一方面，CNV-GABP在不同的腫瘤純度和測序深度中選擇大量訓練數據進行訓練，提高了容錯率；另一方面，其他方法只考慮到多個特征的邊際效應，而CNV-GABP通過BP神經網絡解決多個特征之間的聯合作用；最后，CNV-GABP 使用遺傳算法優化BP 神經網絡模型，增強了模型的魯棒性。

3.2 真實數據研究

為驗證CNV-GABP 的有效性，應用該算法分析3個實際測序樣本：NA19238、NA19239 和NA19240，數據可以在1 000 基因組項目網站http：//www.internationalgome.org/上獲取。本文將5 種檢測方法分別用于這3個樣本，檢測結果如表1 所示。基因組變異數據庫（http：//DGV.tcag.ca/）提供已確認的CNV，可量化這5種方法的性能。

表1 5種方法檢測的CNV個數比較Table 1 Number comparison of detected copy number variations by 5 methods

由表1 展示的結果可以看出，在NA19238 樣本中，CNV-GABP 檢測的CNV 個數少于CNAnator，剩余的兩個樣本中，CNV-GABP 檢測的CNV 個數均多于其他方法。這表明針對特征不明顯的CNV 時，該算法具備較強的檢測能力。5種方法在真實數據上的性能比較結果如圖6所示。

由圖6 可以看出，3 個樣本中CNV-GABP 的F1 評分值均最高，平均值為0.85；就精度而言，MFCNV 表現優于CNV-GABP，其次是CNVnator；就靈敏度而言，CNV-GABP 獲得兩個最高值，其次是MFCNV。因此，CNV-GABP 在實際測序樣本數據中的整體性能表現優于其他4種算法，這說明該算法在實際應用中是有效且可靠的。

圖6 真實數據中5種方法的性能比較Fig.6 Performance comparison of five methods in real data

針對GC 含量偏差校正過程產生的誤差對現有CNV 檢測方法造成影響的問題，本文提出了一種遺傳算法優化的BP神經網絡拷貝數變異檢測算法。該算法無需進行GC含量偏差校正，采用BP神經網絡解決不同特征之間的聯合作用，然后使用遺傳算法優化BP 神經網絡模型獲取最優參數，最后使用優化的BP 神經網絡進行CNV 預測。本文基于仿真數據和1 000 基因組計劃的真實數據，將CNV-GABP與其他4種現有算法進行比較。實驗結果表明，CNV-GABP算法在靈敏度、精度、F1 評分和邊界偏差上均獲得了明顯的優勢，尤其提高了對特征不明顯的CNV的檢測能力。

對于將來的工作，將從以下幾方面著手對CNVGABP 方法進行改進和擴展：（1）將算法擴展至腫瘤-正常配對樣本的研究，對腫瘤樣本和正常樣本分別提取特征，并將提取后的特征作為CNV-GABP算法的輸入，以進行差異基因分析；（2）單細胞測序技術蓬勃發展，下一步可收集此類測序數據的基因組特征，將該算法擴展至單細胞測序數據中，以分析單個細胞的基因組突變；（3）同時可根據混合優化算法[29-31]優化BP神經網絡的結構，提高算法運行效率。

4 結束語