一種二進(jìn)制癌癥單驅(qū)動(dòng)通路識(shí)別模型和算法

摘 要：針對(duì)驅(qū)動(dòng)通路識(shí)別的相關(guān)研究依賴傳統(tǒng)生物實(shí)驗(yàn)方法，存在費(fèi)時(shí)費(fèi)力且經(jīng)濟(jì)成本高的問(wèn)題，提出一種新的二進(jìn)制癌癥驅(qū)動(dòng)通路識(shí)別方法PEA-BLMWS。首先，利用已有的基因表達(dá)數(shù)據(jù)，通過(guò)對(duì)比正常基因與突變基因表達(dá)量的差異，挖掘潛在的基因突變數(shù)據(jù)；其次，引入蛋白質(zhì)相互作用網(wǎng)絡(luò)數(shù)據(jù)，構(gòu)建出一個(gè)改進(jìn)的二進(jìn)制線性最大權(quán)重子矩陣模型；最后，提出一種雙親協(xié)同進(jìn)化算法求解該矩陣模型。在GBM（glioblastoma）和OVCA（ovarian cancer）數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明，相比于其他先進(jìn)的Dendrix、CCA-NMWS和CGP-NCM識(shí)別方法，PEA-BLMWS識(shí)別的基因集中有更多基因富集在已知的信號(hào)通路中，未富集在信號(hào)通路中的基因也與癌癥的發(fā)生密切相關(guān)，故該識(shí)別方法可作為一種驅(qū)動(dòng)通路識(shí)別的有效工具。

關(guān)鍵詞：驅(qū)動(dòng)通路； 基因突變； 基因表達(dá)； 進(jìn)化算法

中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-3695（2024）06-018-1728-07

doi：10.19734/j.issn.1001-3695.2023.10.0497

Binary cancer single-driver pathway identification model and algorithm

Abstract：The researches on driver pathway identification in cancer rely on traditional biological experiments， which have the drawbacks of being time-consuming， labor-intensive and costly. This paper proposed a novel binary cancer driver pathway identification method called PEA-BLMWS（parental evolutionary algorithm-binary linear maximum weight sub-matrix）. Firstly， it utilized the existing gene expression data to uncovered potential gene mutation data by comparing the differences in expression levels between normal and mutated genes. Secondly，it incorporated protein-protein interaction network data to construct an improved binary linear maximum weight sub-matrix model. Finally，it proposed a parental evolutionary algorithm to solve this matrix model. Experimental results on the GBM（glioblastoma） and OVCA（ovarian cancer） datasets show that compared to other advanced identification methods such as Dendrix， CCA-NMWS and CGP-NCM， the gene set identified by PEA-BLMWS has more genes enriched in known signaling pathways， and genes not enriched in signaling pathways are also closely related to the occurrence of cancer. Therefore， this identification method can serve as an effective tool for driving pathway identification.

Key words：driver pathway; gene mutation; gene expression; evolutionary algorithm

0 引言

癌癥的產(chǎn)生和發(fā)展與基因組的突變密切相關(guān)，且癌癥高度受體細(xì)胞基因突變驅(qū)動(dòng)，體細(xì)胞基因突變的類型主要包括單核苷酸變異、拷貝數(shù)變異以及核苷酸序列重復(fù)、插入以及缺失等［1］。并不是所有的基因突變都能導(dǎo)致癌癥的發(fā)生，只有少部分突變能夠使細(xì)胞具有生長(zhǎng)優(yōu)勢(shì)，進(jìn)而對(duì)癌癥的產(chǎn)生與發(fā)展起到促進(jìn)作用，這種突變稱之為驅(qū)動(dòng)突變［2］。對(duì)癌癥發(fā)展起促進(jìn)作用的基因稱為驅(qū)動(dòng)基因［3］。癌癥細(xì)胞中驅(qū)動(dòng)基因的集合稱作驅(qū)動(dòng)通路。不同的驅(qū)動(dòng)基因可能共同作用于相同的生物學(xué)通路，生物學(xué)通路是指在生物體細(xì)胞內(nèi)，一系列生物化學(xué)分子通過(guò)多個(gè)級(jí)聯(lián)反應(yīng)以實(shí)現(xiàn)一系列特定功能的路徑。驅(qū)動(dòng)通路中的一個(gè)基因發(fā)生突變時(shí)，可能干擾整個(gè)生物學(xué)通路的調(diào)控機(jī)制，從而促使癌癥的發(fā)生［4］。故將關(guān)注點(diǎn)轉(zhuǎn)移到生物學(xué)通路水平上，可以幫助科學(xué)家更加了解驅(qū)動(dòng)通路與癌癥的靶向關(guān)系，對(duì)于研究癌癥的發(fā)病機(jī)理以及開(kāi)發(fā)抗癌藥物具有重要意義［5］。然而，通過(guò)生物實(shí)驗(yàn)進(jìn)行驅(qū)動(dòng)通路的識(shí)別需要對(duì)大量基因進(jìn)行測(cè)序，花費(fèi)大量時(shí)間并且經(jīng)濟(jì)成本高昂。隨著高通量測(cè)序技術(shù)的發(fā)展，許多大型規(guī)模癌癥基因測(cè)序項(xiàng)目得以開(kāi)展，例如癌癥基因組圖譜計(jì)劃（The Cancer Genome Atlas，TCGA）［6］和國(guó)際腫瘤基因協(xié)作組（International Cancer Genome Consortium，ICGC）［7］。通過(guò)進(jìn)行基因測(cè)序項(xiàng)目，本文獲取了大量的多組學(xué)數(shù)據(jù)，這為利用計(jì)算方法來(lái)識(shí)別驅(qū)動(dòng)通路提供了可能。與傳統(tǒng)的生物實(shí)驗(yàn)相比，使用計(jì)算方法來(lái)識(shí)別驅(qū)動(dòng)通路具有高效率和低成本的優(yōu)勢(shì)［8］。

1 相關(guān)研究

根據(jù)識(shí)別的驅(qū)動(dòng)通路數(shù)量的不同，現(xiàn)有的驅(qū)動(dòng)通路識(shí)別問(wèn)題分為單驅(qū)動(dòng)通路識(shí)別［9］和協(xié)同驅(qū)動(dòng)通路識(shí)別［10］兩大類。協(xié)同驅(qū)動(dòng)通路識(shí)別融入了基因、miRNA等相關(guān)組學(xué)數(shù)據(jù)，能夠更加全面地分析癌癥發(fā)生的機(jī)制，但存在處理數(shù)據(jù)量大，模型計(jì)算難度大等問(wèn)題［11，12］。與協(xié)同驅(qū)動(dòng)通路識(shí)別相比，盡管單一驅(qū)動(dòng)通路識(shí)別方法僅考慮與基因相關(guān)的信息，會(huì)忽略一些在癌癥發(fā)生中具有重要作用的組學(xué)數(shù)據(jù)，如miRNA等，導(dǎo)致突變矩陣融入的數(shù)據(jù)類型受限，但這種方法可快速確定致癌基因，從而為癌癥治療提供可靠依據(jù)。兩類識(shí)別方法都存在各自的優(yōu)缺點(diǎn)，本文聚焦單驅(qū)動(dòng)通路識(shí)別問(wèn)題進(jìn)行研究。

2012年，Vandin等人［9］首次提出驅(qū)動(dòng)通路具備高覆蓋度和高互斥度兩個(gè)重要性質(zhì)。高覆蓋度是指大量患者樣本被驅(qū)動(dòng)基因覆蓋；高互斥度是指在一個(gè)癌癥患者樣本中，同一驅(qū)動(dòng)通路中的基因同時(shí)發(fā)生突變的概率很小。其提出的最大權(quán)重子矩陣（maximum weight submatrix，MWS）問(wèn)題模型是后續(xù)驅(qū)動(dòng)通路識(shí)別方法的基礎(chǔ)。對(duì)應(yīng)求解MWS問(wèn)題模型的基于馬爾可夫鏈的求解方法Dendrix，僅使用突變數(shù)據(jù)識(shí)別癌癥驅(qū)動(dòng)通路，雖然這種方法在大多數(shù)情況下表現(xiàn)出良好的效果，但在搜索過(guò)程中可能會(huì)受限于局部最優(yōu)解，無(wú)法跳出這種局限性。2012年，Zhao 等人［13］開(kāi)發(fā)了一個(gè)軟件包MDPFinder求解MWS問(wèn)題模型，MDPFinder包括二進(jìn)制線性規(guī)劃方法（binary linear programming，BLP）和基于遺傳算法的求解方法（genetic algorithm，GA）兩種識(shí)別方法。BLP和GA針對(duì)不同類型的權(quán)重函數(shù)及其他組學(xué)數(shù)據(jù)均適用。MDPFinder的局限性是沒(méi)有考慮覆蓋度和互斥度的平衡，如果識(shí)別的兩個(gè)基因集權(quán)重相同，其中覆蓋度和互斥度更加平衡的基因集會(huì)被忽略。2016年，Zheng等人［14］整合基因表達(dá)數(shù)據(jù)和突變數(shù)據(jù)，提出一種基于遺傳算法的多目標(biāo)優(yōu)化方法MOGA求解MWS問(wèn)題模型。MOGA將基因表達(dá)數(shù)據(jù)融入MWS問(wèn)題模型，通過(guò)調(diào)整模型中覆蓋度和互斥度的權(quán)重，改進(jìn)了MWS問(wèn)題模型。相比之前方法，MOGA在調(diào)整高覆蓋度和高互斥度之間的權(quán)衡方面更加可靠，但其局限性是未考慮到基因內(nèi)部之間的聯(lián)系。2021年，Wu等人［15］提出了一種新的識(shí)別模型和算法CGA-MWS，通過(guò)整合拷貝數(shù)變異、體細(xì)胞突變以及基因表達(dá)數(shù)據(jù)，構(gòu)建了一個(gè)加權(quán)的非二進(jìn)制突變矩陣，并重新定義了覆蓋度和互斥度的概念。該模型的局限性在于它僅考慮了與基因相關(guān)的數(shù)據(jù)，而沒(méi)有考慮其他可能的影響因素。2022年，Wu等人［16］提出了一種非線性最大權(quán)重子矩陣識(shí)別模型NMWS，基于覆蓋度、互斥度以及網(wǎng)絡(luò)連通性實(shí)現(xiàn)。NMWS 將癌癥樣本中的突變基因視為無(wú)向圖中彼此連接的頂點(diǎn)，并通過(guò)圖中頂點(diǎn)的聚集度來(lái)衡量驅(qū)動(dòng)通路的互斥性，然后提出一種競(jìng)爭(zhēng)協(xié)同演化算法CCA來(lái)求解NMWS模型。但NMWS模型只研究了基因集內(nèi)部相互之間的聯(lián)系，對(duì)其他網(wǎng)絡(luò)信息考慮不全面。2023年，Zhu等人［17］提出了一種新的癌癥驅(qū)動(dòng)通路識(shí)別方法CGP-NCM，該方法在不預(yù)先設(shè)置人工參數(shù)的情況下構(gòu)建了一個(gè)加權(quán)的非二進(jìn)制突變矩陣，并提出了一種無(wú)參數(shù)的蛋白質(zhì)相互作用網(wǎng)絡(luò)識(shí)別模型，最后，設(shè)計(jì)了一種基于粒子群優(yōu)化算法的協(xié)同進(jìn)化算法，用于求解該模型。

此外，單一組學(xué)數(shù)據(jù)存在噪聲且包含的信息過(guò)于獨(dú)立，導(dǎo)致基于單一組學(xué)數(shù)據(jù)的驅(qū)動(dòng)通路識(shí)別得到的普遍不是最優(yōu)解。因此，針對(duì)現(xiàn)有單驅(qū)動(dòng)通路識(shí)別方法的局限性，本文提出了一種新的二進(jìn)制癌癥驅(qū)動(dòng)通路識(shí)別方法PEA-BLMWS。首先，整合體細(xì)胞突變數(shù)據(jù)、拷貝數(shù)變異數(shù)據(jù)和基因表達(dá)數(shù)據(jù)，生成一個(gè)二進(jìn)制突變矩陣；其次，引入蛋白質(zhì)相互作用網(wǎng)絡(luò)數(shù)據(jù)，構(gòu)建一個(gè)改進(jìn)的二進(jìn)制線性最大權(quán)重子矩陣模型（binary linear maximum weight submatrix，BLMWS），BLMWS模型充分考慮了驅(qū)動(dòng)通路的高覆蓋度和高互斥度的特性、突變基因?qū)ζ渌赐蛔兓虻挠绊憽⒌鞍踪|(zhì)網(wǎng)絡(luò)連通性等因素，提高了基因在樣本中的突變可信度；最后，提出一種改進(jìn)的雙親進(jìn)化算法PEA（parental evolutionary algorithm）求解BLMWS模型。實(shí)驗(yàn)結(jié)果表明，本文提出的PEA-BLMWS可快速高效地識(shí)別出多種具有生物學(xué)意義的單驅(qū)動(dòng)通路。

2 BLMWS模型構(gòu)建

2.1 矩陣定義

使用計(jì)算的方法識(shí)別驅(qū)動(dòng)通路，需要定義一組形式化表述的矩陣符號(hào)，對(duì)應(yīng)表示多種組學(xué)數(shù)據(jù)。一組癌癥樣本集合表示為P，數(shù)量為|P|。一組備選基因集合表示為G，數(shù)量為|G|。

體細(xì)胞突變矩陣S|P|×|G|，矩陣元素的值為sij∈{0，1}（i=1，2，…，|P|，j=1，2，…，|G|），如果第j個(gè)基因在第i個(gè)癌癥樣本中發(fā)生突變，則sij=1，若無(wú)突變，則sij=0。

拷貝數(shù)變異矩陣C|P|×|G|，矩陣元素的取值范圍為cij∈{-2，-1，0，1，2}（i=1，2，…，|P|，j=1，2，…，|G|），cij表示第j個(gè)基因在第i個(gè)癌癥樣本中的拷貝數(shù)變異值。

基因表達(dá)矩陣E|P|×|G|，矩陣元素值為實(shí)數(shù)，矩陣元素eij表示第j個(gè)基因在第i個(gè)癌癥樣本中的表達(dá)量，即eij∈R（i=1，2，…，|P|，j=1，2，…，|G|）。

蛋白質(zhì)互作矩陣Q|G|×|G|，矩陣的行和列均代表基因，其中矩陣元素qij∈{0，1}（i=1，2，…，|G|，j=1，2，…，|G|，i≠j），qij=1表示基因i與j有相互作用，qij=0則表示兩個(gè)基因無(wú)相互作用。

突變矩陣D|P|×|G|，體細(xì)胞突變矩陣S|P|×|G|與拷貝數(shù)變異矩陣C|P|×|G|的交集稱為突變矩陣，即矩陣D|P|×|G|=S|P|×|G|∩C|P|×|G|。其中矩陣元素dij∈{0，1}（i=1，2，…，|P|，j=1，2，…，|G|），若sij和cij取值都為1，則dij=1，反之，則dij=0。

倍數(shù)矩陣T|P|×|G|，矩陣元素tij（i=1，2，…，|P|，j=1，2，…，|G|）的值與基因表達(dá)矩陣E|P|×|G|的矩陣元素eij有關(guān)。設(shè)基因表達(dá)矩陣E|P|×|G|中的第j個(gè)基因的所有正常樣本的基因表達(dá)量的平均值為，則矩陣元素tij=eij/。

加權(quán)突變矩陣A|P|×|G|，為了提高未突變基因的潛在突變可能性，對(duì)矩陣D|P|×|G|中的未突變基因進(jìn)行突變，提升了潛在突變基因的受關(guān)注度，得到一個(gè)加權(quán)突變矩陣A|P|×|G|。矩陣元素aij的值如式（1）所示。

其中：λ是用于確定是否進(jìn)行突變操作的倍數(shù)閾值，經(jīng)大量實(shí)驗(yàn)驗(yàn)證，λ的最優(yōu)值為3。加權(quán)突變矩陣A|P|×|G|的具體構(gòu)建過(guò)程，如圖1所示。

2.2 通路特性定義

根據(jù)先前的研究，癌癥驅(qū)動(dòng)通路最大子矩陣模型具備高覆蓋度和高互斥度兩個(gè)重要特征。本文在基因相關(guān)數(shù)據(jù)和蛋白質(zhì)相互作用的數(shù)據(jù)基礎(chǔ)上，引入相關(guān)度和連通性的概念或指標(biāo)，將癌癥驅(qū)動(dòng)通路的識(shí)別模型轉(zhuǎn)換為一個(gè)多目標(biāo)優(yōu)化問(wèn)題。在這個(gè)問(wèn)題中，通過(guò)參數(shù)K2來(lái)減弱基因相關(guān)度和蛋白質(zhì)連通性對(duì)驅(qū)動(dòng)通路的影響，以優(yōu)化模型的性能。

2.2.1 覆蓋度

設(shè)矩陣A|P|×|G|的任意一個(gè)子矩陣為M|P|×K（1<K<|G|），當(dāng)aij=1（i=1，2，…，|P|）時(shí)，令Γ（j）表示矩陣M|P|×K中第j個(gè)基因發(fā)生突變的樣本集合。令CO（M）表示矩陣M|P|×K的覆蓋度，定義如式（2）所示。

2.2.2 互斥度

假設(shè)p和q代表矩陣M|P|×K中任意兩個(gè)基因（p≠q），當(dāng)aip=1且aiq=1時(shí)，令Γ（pq）表示p和q同時(shí)發(fā)生突變的樣本集合。令MU（M）表示矩陣M|P|×K的互斥度，定義如下：

2.2.3 相關(guān)度

2.2.4 連通性

令PE（M）表示蛋白質(zhì)互作矩陣Q|G|×|G|中基因之間的連通性，定義如式（5）所示。

2.3 構(gòu)建過(guò)程

基于上述定義，定義一個(gè)函數(shù)W（M），在給定的加權(quán)突變矩陣A|P|×|G|中尋找使W（M）函數(shù)值最大的子矩陣M|P|×K，W（M）的定義如式（6）所示。

W（M）=CO（M）+MU（M）+RE（M）+PE（M）（6）

因此，對(duì)應(yīng)使函數(shù)W（M）取得最大值的矩陣M|P|×K中所包含的基因組就是驅(qū)動(dòng)通路識(shí)別問(wèn)題的解。識(shí)別癌癥驅(qū)動(dòng)通路問(wèn)題轉(zhuǎn)換為求解使函數(shù)W（M）取得最大值的子矩陣M|P|×K。

3 PEA算法構(gòu)建

3.1 染色體編碼方案

在遺傳算法中，種群是指用遺傳算法解決問(wèn)題時(shí)，初始時(shí)給定的多個(gè)解的集合。初始時(shí)設(shè)置兩個(gè)種群，每個(gè)種群包含n條染色體，每條染色體代表問(wèn)題的一個(gè)解。染色體編碼一般有十進(jìn)制編碼和二進(jìn)制編碼兩種方式。本文種群中的染色體采用十進(jìn)制編碼方式，在備選基因集合G中隨機(jī)挑選K個(gè)基因，生成初始染色體X={x1，x2，…，xK}（xi∈G，i=1，2，3，…，K）。

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

適應(yīng)度函數(shù)是一種用來(lái)對(duì)種群中任意個(gè)體（染色體）的環(huán)境適應(yīng)性進(jìn)行度量的函數(shù)，其函數(shù)值是遺傳算法實(shí)現(xiàn)優(yōu)勝劣汰的依據(jù)。設(shè)每條染色體X代表子矩陣M|P|×K中的K個(gè)基因，對(duì)應(yīng)的適應(yīng)度函數(shù)F（X）定義如下：

F（X）=W（M）（7）

3.3 選擇算子

選擇算子的作用是對(duì)種群中的個(gè)體進(jìn)行優(yōu)勝劣汰操作，根據(jù)個(gè)體適應(yīng)度函數(shù)值從種群中選出優(yōu)秀的個(gè)體，傳給下一代種群。為了避免算法陷入過(guò)早收斂，選擇算子使用輪盤賭選擇法、錦標(biāo)賽選擇法以及精英策略三種方法結(jié)合的方式實(shí)現(xiàn)以下的流程：

a）對(duì)于每一個(gè)種群，精英策略直接將種群中最優(yōu)個(gè)體從父代傳遞到子代；

b）用輪盤賭和錦標(biāo)賽選擇法從兩個(gè)種群中各選擇一半個(gè)體進(jìn)入子代。

以上三種方法相結(jié)合實(shí)現(xiàn)的選擇算子能夠明顯增加搜索的廣度，雖然算法的收斂速度略有降低，但總體效果較好。

3.4 交叉算子

交叉算子的目的是增強(qiáng)算法的局部搜索能力，對(duì)應(yīng)的實(shí)現(xiàn)流程如下：

3.5 變異算子

變異算子可以增加種群的多樣性，避免算法過(guò)早陷入局部最優(yōu)解，本文設(shè)計(jì)了單點(diǎn)變異算子和概率變異自適應(yīng)算子來(lái)增加種群的多樣性。

3.5.1 單點(diǎn)變異算子

對(duì)于染色體中的基因，執(zhí)行基于貪心策略的單點(diǎn)變異算子，具體流程如下：

a）給定一個(gè)染色體X，隨機(jī)刪除X中一個(gè)基因得到染色體X′，對(duì)應(yīng)的候選基因集為G′={g|g∈G，gX}；

b）從候選基因集G′中隨機(jī)挑選一個(gè)基因加入到X′中；

c）比較F（X）與F（X′），如果執(zhí)行單點(diǎn)變異后的染色體的適應(yīng)度小于變異前染色體適應(yīng)度，則重復(fù)執(zhí)行步驟b）。若候選基因集G′中的基因都不滿足F（X）<F（X′），則染色體X保持不變，否則得到變異后的新染色體X′。

3.5.2 概率變異自適應(yīng)算子

對(duì)于種群中的染色體，采用概率變異算子實(shí)現(xiàn)自適應(yīng)的改變，具體流程如下：

a）種群平均適應(yīng)度是指該種群內(nèi)所有個(gè)體適應(yīng)度值的平均值，設(shè)任一種群的種群平均適應(yīng)度為average_fit，設(shè)有該種群任一染色體X，則其適應(yīng)度值為F（X）。

b）若F（X）≥average_fit，則該染色體突變基因數(shù)量為大于等于0.1×K的最小整數(shù)；若F（X）<average_fit，則該染色體突變基因數(shù)量為小于等于0.5×K的最小整數(shù)；設(shè)最終步驟b）染色體突變數(shù)量為temp。

d）從候選基因集G′中隨機(jī)挑選temp個(gè)基因加入到X′中，若X=X′，則重復(fù)執(zhí)行步驟d），直到得到多點(diǎn)變異后的新染色體X′。

3.6 算法流程

算法：PEA

4 實(shí)驗(yàn)與分析

4.1 實(shí)驗(yàn)設(shè)置

4.1.1 數(shù)據(jù)集預(yù)處理

對(duì)膠質(zhì)母細(xì)胞瘤（GBM）和卵巢癌（OVCA）兩種原始的真實(shí)癌癥數(shù)據(jù)集進(jìn)行相應(yīng)的預(yù)處理，蛋白質(zhì)互作網(wǎng)絡(luò)數(shù)據(jù)取自文獻(xiàn)［19］。

a）原始膠質(zhì)母細(xì)胞瘤數(shù)據(jù)集包括91個(gè)患者樣本的體細(xì)胞突變數(shù)據(jù)，206個(gè)患者樣本的拷貝數(shù)變異數(shù)據(jù)以及259個(gè)患者樣本的基因表達(dá)數(shù)據(jù)。對(duì)冗余數(shù)據(jù)處理后，得到覆蓋90個(gè)患者樣本的440個(gè)基因的基因突變矩陣和基因表達(dá)矩陣，組成實(shí)驗(yàn)所用的新膠質(zhì)母細(xì)胞瘤（GBM）數(shù)據(jù)集。

b）原始卵巢癌數(shù)據(jù)集包括313個(gè)患者樣本的突變數(shù)據(jù)，以及489個(gè)患者樣本的基因表達(dá)數(shù)據(jù)。對(duì)冗余數(shù)據(jù)處理后，得到覆蓋313個(gè)患者樣本的2 547個(gè)基因的突變矩陣和基因表達(dá)矩陣，組成實(shí)驗(yàn)所用的新卵巢癌（OVCA）數(shù)據(jù)集。

4.1.2 參數(shù)設(shè)置

經(jīng)大量實(shí)驗(yàn)驗(yàn)證后，識(shí)別方法PEA-BLMWS對(duì)應(yīng)的超參數(shù)值設(shè)置如表1所示。

4.1.3 實(shí)驗(yàn)環(huán)境

實(shí)驗(yàn)在一臺(tái)電腦（AMD Ryzen 7 5800H，3.20 GHz，內(nèi)存16 GB）上進(jìn)行，編譯運(yùn)行工具為R 4.2.1。

4.2 實(shí)驗(yàn)結(jié)果

設(shè)置參數(shù)K從2變化到8，將三種基準(zhǔn)方法（Dendrix、CCA-NMWS、CGP-NCM）和PEA-BLMWS方法識(shí)別到的基因集作比較，GBM數(shù)據(jù)集識(shí)別的結(jié)果見(jiàn)表2，OVCA數(shù)據(jù)集識(shí)別的結(jié)果見(jiàn)表3。圖2和4為基因集對(duì)應(yīng)信號(hào)通路圖，圖中實(shí)線表示抑制作用，實(shí)線箭頭表示直接激活作用，虛線箭頭表示間接激活作用，有圖案填充的橢圓代表識(shí)別的基因，無(wú)填充的橢圓表示與識(shí)別到的基因相關(guān)聯(lián)的基因。

4.2.1 GBM數(shù)據(jù)集實(shí)驗(yàn)結(jié)果

GBM識(shí)別結(jié)果如表2所示。結(jié)果分析如下：

a）基因集分析

①當(dāng)K=2時(shí)，如圖2（a）所示，三種基準(zhǔn)方法識(shí)別到的基因集（CDKN2B CDK4）都富集在RB信號(hào)通路中， PEA-BLMWS可以識(shí)別出其他方法未識(shí)別到的基因集（CDKN2A CDK4）， （CDKN2A CDK4） 是P53信號(hào)通路的一部分。

②當(dāng)K=3時(shí)，如圖2（a）所示，Dendrix識(shí)別的基因集（CDKN2B CDK4 RB1）富集在RB信號(hào)通路中。如圖2（b）所示，CCA-NMWS、CGP-NCM、PEA-BLMWS識(shí)別的基因集（CDKN2A MDM2 TP53） 同樣富集在P53信號(hào)通路中。相關(guān)研究表明，MDM2基因在P53信號(hào)通路中失調(diào)會(huì)抑制P53信號(hào)通路在抑制癌癥方面的作用［20］。

③當(dāng)K=4時(shí)，如圖2（b）所示，CCA-NMWS和CGP-NCM識(shí)別到的基因集（CDKN2A MDM2 MDM4 TP53）是P53信號(hào)通路的一部分。PEA-BLMWS識(shí)別到不同的基因集（CDKN2A CDK4 MDM2 TP53），該基因集同樣富集在P53信號(hào)通路中。P53信號(hào)通路的失調(diào)與膠質(zhì)母細(xì)胞瘤細(xì)胞的增殖、侵襲、轉(zhuǎn)移、凋亡以及癌細(xì)胞干性有關(guān)，P53蛋白的致癌變體對(duì)膠質(zhì)母細(xì)胞瘤有促進(jìn)作用［21］。

④當(dāng)K=5時(shí)，如圖2（b）所示，CCA-NMWS、CGP-NCM和PEA-BLMWS識(shí)別的基因集中，只有基因RB1與其他基因沒(méi)有富集在同一個(gè)信號(hào)通路中。RB1是一個(gè)重要的癌基因，當(dāng)其表達(dá)產(chǎn)物超過(guò)正常水平時(shí)會(huì)加快細(xì)胞增殖速度并降低細(xì)胞活性［22］。

⑤當(dāng)K=6、7時(shí)，如圖2（b）所示，除了基因RB1，CCA-NMWS、CGP-NCM與PEA-BLMWS識(shí)別到的基因集均富集在P53信號(hào)通路中。基因MDM4與其他三種基因CDKN2ATP53和MDM2共同參與P53信號(hào)通路。基因MDM4與MDM2具有相似的結(jié)構(gòu)和功能，它們都在P53信號(hào)通路中充當(dāng)致癌抑制劑以抑制腫瘤活性［21］。

⑥當(dāng)K=8時(shí)，如圖2（c）所示， PEA-BLMWS識(shí)別到一個(gè)不同的基因集（CDK4 EGFR ERBB2 PIK3R1 PTEN TP53），該基因集富集在PI3K-AKT信號(hào)通路中。PTEN基因在抑制腫瘤的產(chǎn)生和發(fā)展中起到重要作用，當(dāng)PTEN基因發(fā)生缺失或突變時(shí)，會(huì)導(dǎo)致細(xì)胞增殖加速，并減緩細(xì)胞的凋亡速度。這些變化可以促進(jìn)膠質(zhì)母細(xì)胞瘤的發(fā)展，PTEN和TP53的失活會(huì)對(duì)膠質(zhì)母細(xì)胞瘤的發(fā)展起到促進(jìn)作用［23］。

b）通路特性分析

在膠質(zhì)母細(xì)胞瘤癌癥樣本中，本文PEA-BLMWS方法識(shí)別的基因集在不同K值的覆蓋和互斥下的情況如圖3所示。PEA-BLMWS方法在基因集的覆蓋度方面表現(xiàn)出色。它能夠覆蓋絕大多數(shù)樣本中的突變基因，并且絕大部分突變基因之間呈現(xiàn)互斥的特點(diǎn)。驅(qū)動(dòng)通路具有高覆蓋度和高互斥度這兩個(gè)重要性質(zhì)，因此，根據(jù)覆蓋度和互斥度的表現(xiàn)，可以認(rèn)為PEA-BLMWS方法在GBM數(shù)據(jù)集上表現(xiàn)良好。

4.2.2 OVCA數(shù)據(jù)集實(shí)驗(yàn)結(jié)果

OVCA識(shí)別結(jié)果如表3所示。結(jié)果分析如下：

a）基因集分析

由表3可知，當(dāng)K=2，3，4時(shí)，Dendrix識(shí)別的基因集只有2個(gè)基因富集在同一信號(hào)通路，當(dāng)K=5，6，7，8時(shí)，Dendrix識(shí)別的基因集只有3個(gè)基因富集在同一信號(hào)通路。與CCA-NMWS、CGP-NCM和PEA-BLMWS方法相比， Dendrix方法識(shí)別的基因集中富集在同一信號(hào)通路中的基因較少。

①當(dāng)K=3時(shí)，如圖4（b）所示，PEA-BLMWS識(shí)別到的基因集中，有2個(gè)基因（BRCA1 MYC）富集在PI3K-AKT信號(hào)通路。PI3K-AKT信號(hào)通路在各種癌細(xì)胞的遷移中起著重要作用，研究表明，至少15%的侵襲性卵巢癌是由基因BRCA1和BRCA2突變引起的［24］。SMARCA4基因是一個(gè)重要的基因，SMARCA4表達(dá)缺失會(huì)造成卵巢高鈣血癥小細(xì)胞的癌變［25］。

②當(dāng)K=4時(shí)，如圖4（b）所示，在PEA-BLMWS識(shí)別的基因集中，有3個(gè)基因（BRCA1 CCNE1 MYC）富集在PI3K-AKT信號(hào)通路，而其他三種方法識(shí)別的基因集只有2個(gè)基因富集在PI3K-AKT信號(hào)通路。

③當(dāng)K=5時(shí)，如圖4（b）所示，CCA-NMWS、CGP-NCM與PEA-BLMWS識(shí)別的基因集中，有3個(gè)基因（BRCA1 CCNE1 MYC）富集在PI3K-AKT信號(hào)通路中。其中新識(shí)別的基因FANCA具有腫瘤抑制作用，其等位基因突變可能增加卵巢癌的患病幾率［26］。當(dāng)K>4時(shí)，Dendrix識(shí)別的基因集有且只有3個(gè)基因富集在同一信號(hào)通路，富集在同一信號(hào)通路中基因的數(shù)量并沒(méi)有隨著K的增大而增加。

④當(dāng)K=6時(shí)，如圖4（b）所示，CCA-NMWS方法識(shí)別的基因中只有3個(gè)基因富集在同一個(gè)信號(hào)通路。如圖4（a）（c）所示，CGP-NCM與PEA-BLMWS識(shí)別的基因中有4個(gè)基因（BRCA1 CCNE1 MYC KRAS）富集在PI3K-AKT信號(hào)通路中，其中基因（MYC KRAS）也在MAPK以及ERBB信號(hào)通路中富集。MAPK信號(hào)通路在卵巢癌發(fā)生中起重要作用，并參與卵巢癌細(xì)胞的遷移［27］。具有YAP的ERBB信號(hào)通路形成自分泌環(huán)以誘導(dǎo)卵巢癌的發(fā)生和發(fā)展［28］。

⑤當(dāng)K=7時(shí)，如圖4（b）所示，CCA-NMWS、CGP-NCM與PEA-BLMWS方法識(shí)別的基因集有4個(gè)基因富集在同一個(gè)信號(hào)通路中。當(dāng)K=8時(shí)， PEA-BLMWS識(shí)別的基因集在信號(hào)通路富集情況優(yōu)于其他方法。CCA-NMWS與CGP-NCM方法識(shí)別到的基因集只有4個(gè)基因富集在同一信號(hào)通路，如圖4（b）所示。PEA-BLMWS識(shí)別到的基因集有5個(gè)基因（CCNE1 EGFR GRB2 MYC PIK3R1）富集在PI3K-AKT信號(hào)通路。如圖4（c）～（e）所示，基因（EGFR GRB2 MYC）富集在MAPK信號(hào)通路中；基因（EGFR GRB2 GAB2 PIK3R1）富集在RAS信號(hào)通路中；基因（EGFR PIK3R1）富集在RAP1信號(hào)通路中。RAS信號(hào)通路在驅(qū)動(dòng)細(xì)胞增殖中起著重要作用，該信號(hào)通路的失調(diào)通常發(fā)生在腫瘤發(fā)生過(guò)程中［29］。RAP1通路通過(guò)促進(jìn)MMP-2和MMP-9的分泌來(lái)調(diào)節(jié)卵巢癌的侵襲和轉(zhuǎn)移［30］。

b）通路特性分析

在卵巢癌樣本中，PEA-BLMWS方法識(shí)別的基因集在不同K值的覆蓋和互斥下的情況如圖5所示。它能夠覆蓋大多數(shù)樣本中的突變基因，且絕大部分突變基因之間呈現(xiàn)互斥的特點(diǎn)。因此，可以認(rèn)為PEA-BLMWS方法在OVCA數(shù)據(jù)集上表現(xiàn)良好。

4.3 顯著性檢驗(yàn)

為了評(píng)估實(shí)驗(yàn)結(jié)果的顯著性并衡量其與理論最優(yōu)值之間的偏離程度，本文采用了隨機(jī)檢驗(yàn)方法。設(shè)識(shí)別到的K個(gè)基因的函數(shù)值為W（M），W（Mi）為隨機(jī)選取的K個(gè)基因的函數(shù)值，重復(fù)隨機(jī)挑選1 000次，統(tǒng)計(jì)W（Mi）>W（M）的次數(shù)，顯著性檢驗(yàn)的公式如下：

經(jīng)大量實(shí)驗(yàn)，p=0.04時(shí)，實(shí)驗(yàn)結(jié)果是顯著的。

5 結(jié)束語(yǔ)

本文提出了一種單驅(qū)動(dòng)通路識(shí)別方法，不依賴于先驗(yàn)知識(shí)，而是采用從頭識(shí)別的方法。該方法基于體細(xì)胞突變數(shù)據(jù)、拷貝數(shù)變異數(shù)據(jù)和基因表達(dá)數(shù)據(jù)構(gòu)建了一個(gè)加權(quán)的二進(jìn)制突變矩陣。在確保驅(qū)動(dòng)通路具有高覆蓋度和高互斥度的基礎(chǔ)上，引入了基因之間的相互作用和蛋白質(zhì)互作網(wǎng)絡(luò)，構(gòu)建了一個(gè)新的識(shí)別模型。為了解決該模型，提出了一種雙親遺傳算法PEA。該算法使用錦標(biāo)賽選擇法和輪盤賭選擇的組合作為選擇算子，并基于貪心策略設(shè)計(jì)了交叉算子，從而能夠高效地找到最優(yōu)解。最后，通過(guò)在膠質(zhì)母細(xì)胞瘤和卵巢癌數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果表明，與其他基準(zhǔn)方法相比，本文PEA-BLMWS在識(shí)別基因集時(shí)能夠更好地捕捉到富集在已知信號(hào)通路中的基因，而那些在信號(hào)通路中沒(méi)有富集的基因已經(jīng)被證實(shí)與癌癥的發(fā)生和發(fā)展密切相關(guān)。為了進(jìn)一步提升模型的性能和可靠性，本文計(jì)劃在更大規(guī)模的數(shù)據(jù)集上展開(kāi)實(shí)驗(yàn)。同時(shí)，本文也將結(jié)合當(dāng)前機(jī)器學(xué)習(xí)和深度學(xué)習(xí)的知識(shí)，探索可能存在于驅(qū)動(dòng)通路中的其他特征，以進(jìn)一步完善驅(qū)動(dòng)通路識(shí)別模型。

參考文獻(xiàn)：

［1］Greenman C， Stephens P， Smith R，et al. Patterns of somatic mutation in human cancer genomes［J］. Nature， 2007，446（7132）：153-158.

［2］Hou J P， Ma Jian. DawnRank： discovering personalized driver genes in cancer［J］. Genome Medicine， 2014，6： article No.56.

［3］Ding Li， Getz G， Wheeler D A， et al. Somatic mutations affect key pathways in lung adenocarcinoma［J］. Nature， 2008，455（7216）： 1069-1075.

［4］Hahn W C， Weinberg R A. Modelling the molecular circuitry of can-cer［J］. Nature Reviews Cancer， 2002，2（5）： 331-341.

［5］Boca S M， Kinzler K W， Velculescu V E， et al. Patient-oriented gene set analysis for cancer mutation data［J］. Genome Biology， 2010，11（11）： article No.R112.

［6］Hudson T J， Anderson W， Aretz A， et al. International network of cancer genome projects［J］. Nature， 2010，464（7291）： 993-998.

［7］Chin L， Andersen J N， Futreal P A. Cancer genomics： from disco-very science to personalized medicine［J］. Nature Medicine， 2011，17（3）： 297-303.

［8］Zhang Junhua， Zhang Shihua. The discovery of mutated driver pathways in cancer： models and algorithms［J］. IEEE/ACM Trans on Computational Biology and Bioinformatics， 2016，15（3）： 988-998.

［9］Vandin F， Upfal E， Raphael B J. De novo discovery of mutated dri-ver pathways in cancer［J］. Genome research， 2012，22（2）： 375-385.

［10］Zhang Junhua， Wu Lingyun， Zhang Xiangsun， et al. Discovery of co-occurring driver pathways in cancer［J］. BMC Bioinformatics， 2014，15（1）： 271.

［11］Wang Jun， Yang Ziying， Domeniconi C， et al. Cooperative driver pathway discovery via fusion of multi-relational data of genes， mi-RNAs and pathways［J］. Briefings in Bioinformatics， 2021，22（2）： 1984-1999.

［12］Yang Ziying， Yu Guoxian， Guo Maozu， et al. CDPath： cooperative driver pathways discovery using integer linear programming and Markov clustering［J］. IEEE/ACM Trans on Computational Bio-logy and Bioinformatics， 2019，18（4）： 1384-1395.

［13］Zhao Junfei， Zhang Shihua， Wu Lingyun， et al. Efficient methods for identifying mutated driver pathways in cancer［J］. Bioinformatics， 2012，28（22）： 2940-2947.

［14］Zheng Chunhou， Yang Wu， Chong Yanwen， et al. Identification of mutated driver pathways in cancer using a multi-objective optimization model［J］. Computers in Biology and Medicine， 2016，72： 22-29.

［15］Wu Jingli， Zhu Kai， Li Gaoshi， et al. A model and algorithm for identifying driver pathways based on weighted non-binary mutation matrix［J］. Applied Intelligence， 2022，52（1）： 127-140.

［16］Wu Jingli， Chen Xiaorong， Li Gaoshi， et al. A nonlinear model and an algorithm for identifying cancer driver pathways［J］. Applied Soft Computing， 2022， 129： 109578.

［17］Zhu Kai， Wu Jingli， Li Gaoshi， et al. A model and cooperative co-evolution algorithm for identifying driver pathways based on the integrated data and PPI network［J］. Expert Systems with Applications， 2023， 212： 118753.

［18］Cohen J， Chen Jingdong， Huang Yiteng， et al. Pearson correlation coefficient［M］//Noise Reduction in Speech Processing. Berlin： Springer， 2009： 1-4.

［19］Leiserson M D M， Vandin F， Wu H T， et al. Pan-cancer network analysis identifies combinations of rare somatic mutations across pathways and protein complexes［J］. Nature Genetics， 2015， 47（2）： 106-114.

［20］Iwakuma T， Lozano G. MDM2， an introduction［J］. Molecular Cancer Research， 2003，1（14）： 993-1000.

［21］Zhang Ying， Dube C， Gibert Jr M， et al. The p53 pathway in glioblastoma［J］. Cancers， 2018，10（9）： 297.

［22］Goldhoff P， Clarke J， Smirnov I， et al. Clinical stratification of glioblastoma based on alterations in retinoblastoma tumor suppressor protein（RB1） and association with the proneural subtype［J］. Journal of Neuropathology & Experimental Neurology， 2012，71（1）： 83-89.

［23］McLendon R， Friedman A， Bigner D，et al. Comprehensive genomic characterization defines human glioblastoma genes and core pathways［J］. Nature， 2008， 455（7216）： 1061-1068.

［24］Pal T， Permuth-Wey J， Betts J A， et al. BRCA1 and BRCA2 mutations account for a large proportion of ovarian carcinoma cases［J］. Cancer： Interdisciplinary International Journal of the American Cancer Society， 2005，104（12）： 2807-2816.

［25］Moes-Sosnowska J，Szafron L，Nowakowska D，et al.Germline SMARCA4 mutations in patients with ovarian small cell carcinoma of hypercalcemic type［J］. Orphanet Journal of Rare Diseases， 2015，10： article No.32.

［26］Thompson E， Dragovic R L， Stephenson S A， et al. A novel duplication polymorphism in the FANCA promoter and its association with breast and ovarian cancer［J］. BMC Cancer， 2005，5： article No. 43.

［27］Ptak A， Hoffmann M， Gruca I， et al. Bisphenol a induce ovarian cancer cell migration via the MAPK and PI3K/Akt signalling pathways［J］. Toxicology Letters， 2014，229（2）： 357-365.

［28］He Chunbo， Lyu Xiangmin， Hua Guohua， et al. YAP forms autocrine loops with the ERBB pathway to regulate ovarian cancer initiation and progression［J］. Oncogene， 2015，34（50）： 6040-6054.

［29］Punekar S R， Velcheti V， Neel B G， et al. The current state of the art and future trends in RAS-targeted cancer therapies［J］. Nature Reviews Clinical Oncology， 2022，19（10）： 637-655.

［30］Che Yaling， Luo Shujuan， Li Gang， et al. The C3G/Rap1 pathway promotes secretion of MMP-2 and MMP-9 and is involved in serous ovarian cancer metastasis［J］. Cancer Letters， 2015， 359（2）： 241-249.