顧云婧, 朱 平

(江南大學 理學院， 無錫 214122)

卵巢癌是最致命的女性疾病之一，2012年全球共診斷出約238 700例新增卵巢癌患者，居于女性癌癥死因的第4位[1]。由于潛在癥狀不明顯，且缺乏可靠的早期篩查方法，約70%卵巢癌患者被診斷時已為晚期(III-IV期)，晚期患者的5年生存率低于30%[2]。目前，卵巢癌患者的治療主要通過手術(shù)切除和化療[3]，但大多數(shù)晚期患者會在18個月內(nèi)復(fù)發(fā)[4]，預(yù)后情況并不理想。因此，這就需要對卵巢癌的早期診斷和預(yù)后標志物進行更深入的研究，將其應(yīng)用于卵巢癌的治療指導(dǎo)和患者管理，以開發(fā)更有效的治療方法來改善卵巢癌的預(yù)后效果。

近30年，基因、蛋白質(zhì)、代謝等組學的迅速發(fā)展以及對腫瘤生物學的深入研究，目前已經(jīng)發(fā)現(xiàn)了VEGF、HE4等數(shù)百種具有預(yù)后價值的卵巢癌生物標志物[5]。因此，臨床上迫切需要對卵巢癌預(yù)后標志物進行系統(tǒng)分析，以便于指導(dǎo)臨床治療策略，基因相互作用網(wǎng)絡(luò)則為系統(tǒng)分析提供了可能。研究表明，基因相互作用網(wǎng)絡(luò)可以用于探究功能基因簇(功能模塊)和預(yù)后因素之間可能存在的內(nèi)在聯(lián)系[6]。通過加權(quán)基因共表達網(wǎng)絡(luò)分析方法(Weighted gene co-expression network analysis，WGCNA)構(gòu)建穩(wěn)健的基因共表達網(wǎng)絡(luò)，其中有意義的模塊可用于推斷腫瘤機制，預(yù)測患者存活率，以及建立新的診斷或治療目標。WGCNA方法為系統(tǒng)生物學提供了功能性解釋工具，目前已運用于乳腺癌、子宮內(nèi)膜癌等腫瘤中[7-10]。在卵巢癌中，研究人員應(yīng)用WGCNA研究了TP53錯義和無效突變[11]。

為了提高預(yù)后標志物的準確性，在構(gòu)建共表達網(wǎng)絡(luò)時融合基因的蛋白質(zhì)相互作用(Protein-protein interaction, PPI)數(shù)據(jù)，從而將基因的生物學意義納入網(wǎng)絡(luò)。此外，普遍認為拓撲上重要的(例如高度連接的)基因在疾病的發(fā)展中具有重要功能，并且傾向于在患者中表現(xiàn)出一致的表達變異[12]。因此，選擇拓撲重要基因作為預(yù)測因子構(gòu)建預(yù)后模型，能夠使得模型的性能在以下兩個方面得到改善：增強預(yù)后模型的穩(wěn)健性，因為拓撲上重要的基因具有一致的表達變化，從而降低了結(jié)果的偶然性；提高所鑒定的基因與給定疾病的生物學相關(guān)性。

本文運用篩選得到的預(yù)后相關(guān)基因構(gòu)建了卵巢癌預(yù)后加權(quán)基因共表達網(wǎng)絡(luò)，并對網(wǎng)絡(luò)進行系統(tǒng)分析。對于網(wǎng)絡(luò)中的卵巢癌預(yù)后模塊，整合模塊中基因的PPI數(shù)據(jù)對模塊網(wǎng)絡(luò)進行重新加權(quán)，得到預(yù)后模塊重加權(quán)共表達網(wǎng)絡(luò)(Reweighted co-expression network for module prognosis-related genes, RMPCEN)，并利用RMPCEN中基因的拓撲性質(zhì)對模塊基因進行排序。在同時考慮模型預(yù)測能力和模型基因個數(shù)的條件下，選擇合適數(shù)量的生物標志物構(gòu)建最終預(yù)后模型，預(yù)測卵巢癌患者的預(yù)后情況。

1 材料和方法

1.1 數(shù)據(jù)獲取

卵巢癌患者的研究數(shù)據(jù)從癌癥基因組數(shù)據(jù)庫(The cancer genome atlas, TCGA)(https://portal.gdc.cancer.gov/)中下載，TCGA數(shù)據(jù)庫收錄了眾多癌癥的實驗數(shù)據(jù)和臨床數(shù)據(jù)。本文從數(shù)據(jù)庫中下載了卵巢癌患者的全基因高通量轉(zhuǎn)錄組測序數(shù)據(jù)和臨床信息，并對編碼基因表達量進行歸一化處理。最終選取了320個同時具有基因表達數(shù)據(jù)和臨床信息的卵巢癌樣本，且生存時間均大于30 d。

1.2 篩選預(yù)后相關(guān)基因

對于患者的臨床數(shù)據(jù)，由于每個病例隨訪時間不同，發(fā)生終點事件的可能性也不相等。因此，本文選用Cox單因素回歸分析[13]篩選與總體生存(Overall survival, OS)相關(guān)的基因，作為卵巢癌預(yù)后相關(guān)基因。對于每個基因，根據(jù)其表達值進行中值分割，將患者分為高表達和低表達兩組進行單因素回歸分析?；貧w分析采用R軟件進行，P<0.05認為有統(tǒng)計學意義。

1.3 預(yù)后基因共表達網(wǎng)絡(luò)的構(gòu)建方法

本文采用加權(quán)的基因共表達網(wǎng)絡(luò)分析方法構(gòu)建卵巢癌預(yù)后網(wǎng)絡(luò)，并對網(wǎng)絡(luò)進行模塊挖掘。

運用皮爾遜相關(guān)系數(shù)衡量基因之間的相關(guān)性，得到基因表達的相關(guān)性矩陣：

S=[Smn=|cor(m,n)|]

(1)

其中，Smn表示基因m和n的相關(guān)性系數(shù)。

選擇適當?shù)募訖?quán)系數(shù)β作為軟閾值，運用冪指數(shù)函數(shù)計算基因之間的鄰接系數(shù)，將相關(guān)矩陣轉(zhuǎn)換為鄰接矩陣，使得網(wǎng)絡(luò)節(jié)點之間的連接符合無尺度網(wǎng)絡(luò)分布，即符合冪律分布，鄰接矩陣:

A=[amn=power(smn,β)]=[|smn|β]

(2)

其中，amn表示基因m和n的鄰接系數(shù)。

基于選定的軟閾值，進一步地將鄰接矩陣轉(zhuǎn)化為拓撲重疊矩陣(Topological overlap matrix, TOM)[14]：

(3)

根據(jù)計算得到的基因拓撲重疊矩陣，運用層次聚類對網(wǎng)絡(luò)進行劃分，并識別網(wǎng)絡(luò)中的模塊。同時，運用Benjamini-Hochberg法校正模塊中基因?qū)ο嚓P(guān)系數(shù)P值，選擇P<0.01的基因?qū)?gòu)建模塊網(wǎng)絡(luò)。

1.4 模塊網(wǎng)絡(luò)的重加權(quán)及網(wǎng)絡(luò)分析方法

為了提高預(yù)測效率，本文考慮模塊網(wǎng)絡(luò)的生物學意義，構(gòu)建模塊網(wǎng)絡(luò)時通過融合基因的PPI數(shù)據(jù)對網(wǎng)絡(luò)進行重新加權(quán)。從STRING[15]數(shù)據(jù)庫(https://string-db.org)中下載了模塊基因的蛋白質(zhì)相互作用數(shù)據(jù)，利用PPI的置信度對網(wǎng)絡(luò)重新加權(quán)：

(4)

對于重加權(quán)的模塊網(wǎng)絡(luò)進行拓撲特征分析，驗證其是否屬于小世界網(wǎng)絡(luò)。網(wǎng)絡(luò)特征分析包含網(wǎng)絡(luò)的特征路徑長度(Characteristic path length, CPL)、聚類系數(shù)(Clustering coefficient, CC)和小世界指數(shù)(Small-world index, SW)。網(wǎng)絡(luò)的CPL為網(wǎng)絡(luò)中所有節(jié)點對之間最短路徑長度的平均值，節(jié)點的CC為該節(jié)點直接連接的節(jié)點中實際存在的邊與所有可能存在的邊的比例，網(wǎng)絡(luò)的CC為網(wǎng)絡(luò)中所有單個節(jié)點的CC平均值。通過比較模塊網(wǎng)絡(luò)和相同規(guī)模的隨機網(wǎng)絡(luò)，網(wǎng)絡(luò)參數(shù)CPL和CC可以定義小世界網(wǎng)絡(luò)[16]，當其滿足以下條件時：

CCsmall-world?CCrandom

(5)

CPLsmall-world≥CPLrandom

(6)

其中，CCrandom和CPLrandom是1000個與目標網(wǎng)絡(luò)相同規(guī)模的Erdos-Rényi(ER)隨機網(wǎng)絡(luò)的平均值。

進一步地，網(wǎng)絡(luò)小世界指數(shù)被定義為[17]：

σSW=(CC/CCrandom)/(CPL/CPLrandom)

(7)

當網(wǎng)絡(luò)的小世界指數(shù)σSW大于1時，該網(wǎng)絡(luò)為小世界網(wǎng)絡(luò)。

同時，基于小世界網(wǎng)絡(luò)的生物學特性，對網(wǎng)絡(luò)進行生物功能分析。通過DAVID(https://david.ncifcrf.gov/home.jsp)在線工具[18]對模塊網(wǎng)絡(luò)包含的基因進行GO(Gene ontology)功能富集的生物過程(Biological process, BP)分析。

1.5 生物標志物的選取及預(yù)后模型的構(gòu)建方法

為了精確識別卵巢癌預(yù)后生物標志物，需要選擇能夠準確預(yù)測生存結(jié)局的數(shù)量最少的基因。根據(jù)網(wǎng)絡(luò)節(jié)點的度(degree)、緊密中心度(Closeness centrality, CCL)、中介中心度(Betweenness centrality, BC)3個拓撲性質(zhì)來衡量節(jié)點在網(wǎng)絡(luò)中的重要性，對節(jié)點進行重要性排序。節(jié)點的度為該節(jié)點與其他節(jié)點直接連接的邊數(shù)，CCL衡量了節(jié)點到其他所有節(jié)點之間的距離，BC用來衡量節(jié)點在網(wǎng)絡(luò)信息傳播中的作用。按以下公式計算節(jié)點的得分，對節(jié)點進行排序：

score=BC×CCL×log2degree

(8)

根據(jù)score選出得分最高的10個基因作為候選預(yù)后標志物，根據(jù)排名依次組合構(gòu)建預(yù)測模型。預(yù)測模型的構(gòu)建采用Cox比例風險模型進行，同時考慮各模型預(yù)測能力和模型基因個數(shù)，選擇最優(yōu)模型作為最終預(yù)測模型。利用Cox回歸構(gòu)建風險評分模型[19]如下：

(9)

其中N是候選預(yù)后標志物的數(shù)量，xi是預(yù)后標志物的表達水平，wi是Cox回歸分析中候選預(yù)后標志物的回歸系數(shù)。

根據(jù)風險評分Riskscore的中值將樣本分為高低兩類，運用KM(Kaplan-meier)圖表估計生存函數(shù)，分析兩組之間的預(yù)后差異。ROC(Receiver operating characteristic)曲線[20]和AUC(Area under curve)值被用來評價預(yù)后預(yù)測模型的效率，AUC值越大的分類器，預(yù)測模型的正確率越高。

2 結(jié)果及分析

2.1 卵巢癌預(yù)后相關(guān)基因

利用Cox回歸對所有卵巢癌基因進行單因素回歸分析，從17 779個基因中共篩選出了747個與OS顯著相關(guān)的基因，稱為卵巢癌預(yù)后相關(guān)基因。在這些預(yù)后相關(guān)基因中，一部分基因在先前研究中已被證實與卵巢癌預(yù)后相關(guān)，如：KLK4[21](HR=1.385;95%CI=1.04～1.85;logrankP=0.027)、CXCL9[22](HR=0.897;95%CI=0.83～0.97;logrankP=0.0034)。

2.2 構(gòu)建的共表達網(wǎng)絡(luò)及網(wǎng)絡(luò)模塊

為使得構(gòu)建的共表達網(wǎng)絡(luò)呈現(xiàn)無尺度網(wǎng)絡(luò)分布，基于不同β條件下的log(k)與log(P(k))(k表示節(jié)點的連接度)的相關(guān)系數(shù)進行加權(quán)系數(shù)的選擇。根據(jù)無尺度網(wǎng)絡(luò)規(guī)則，加權(quán)系數(shù)β需滿足log(k)與log(P(k))呈負相關(guān)，且相關(guān)系數(shù)越大，網(wǎng)絡(luò)的無尺度特征越顯著。最終選擇相關(guān)系數(shù)大于0.9時，將共表達網(wǎng)絡(luò)軟閾值β設(shè)置為3(圖1)，此時網(wǎng)絡(luò)節(jié)點度符合冪律分布，R2=0.89(圖2)。在選定軟閾值β=3后，將747個基因的相關(guān)性系數(shù)矩陣轉(zhuǎn)換為鄰接矩陣，進而轉(zhuǎn)換為拓撲重疊矩陣，并運用層次聚類從中識別了一個預(yù)后模塊(圖3)。

圖1 確定加權(quán)系數(shù)β

圖2 無尺度網(wǎng)絡(luò)特性檢驗

圖3 基因聚類樹狀圖

2.3 重加權(quán)網(wǎng)絡(luò)拓撲分析及富集分析

對于基因共表達網(wǎng)絡(luò)篩選得到的預(yù)后模塊，用模塊基因的蛋白質(zhì)相互作用數(shù)據(jù)進行重新加權(quán)。從STRING數(shù)據(jù)庫中下載得到了1154個模塊基因的PPI，融合PPI信息后得到預(yù)后模塊重加權(quán)共表達網(wǎng)絡(luò)，該網(wǎng)絡(luò)共包含96個節(jié)點和1768條邊。進一步對網(wǎng)絡(luò)進行拓撲分析，驗證RMPCEN的小世界網(wǎng)絡(luò)屬性。拓撲分析表明，RMPCEN聚類系數(shù)為0.848，特征路徑長度為1.614。同樣節(jié)點數(shù)和邊數(shù)的ER隨機網(wǎng)絡(luò)的平均聚類系數(shù)為0.389，平均特征路徑長度為1.610。與ER隨機網(wǎng)絡(luò)相比，RMPCEN滿足公式(5)和公式(6)，因此RMPCEN符合小世界網(wǎng)絡(luò)特性，其小世界指數(shù)為2.173。同時，輸出RMPCEN中權(quán)重大于0.3的基因?qū)Γ肅ytoscape軟件[23]繪制了RMPCEN的權(quán)重網(wǎng)絡(luò)圖(圖4)。

圖4 預(yù)后模塊重加權(quán)共表達網(wǎng)絡(luò)

為了進一步驗證RMPCEN在卵巢癌中的潛在功能意義，對RMPCEN中的基因進行了生物學過程富集分析。富集分析表明，RMPCEN中包含的262個基因在387條GO term中顯著富集(FDR<0.05)，F(xiàn)DR值小于1×10-10的生物學過程見表1。可以看出，RMPCEN中的基因集中參與了免疫系統(tǒng)過程、免疫反應(yīng)、調(diào)節(jié)免疫系統(tǒng)過程等與免疫相關(guān)的生物學過程。

2.4 預(yù)后標志物的選取及預(yù)后模型的建立

計算RMPCEN中每個基因的score得分，并將基因按照score得分排序?？紤]預(yù)后模型中基因的個數(shù)，選擇了得分最高的10個基因測試其預(yù)測性能。排名前10的基因中，3個基因已被證實在卵巢癌的生長、侵襲、轉(zhuǎn)移等過程中起重要作用，可用作卵巢癌治療的良好靶點。這3個基因分別為：TBX21[24]、CXCR6[25-26]和TIGIT[27]。對于這10個候選預(yù)后標志物，根據(jù)score排名分別對于前n(1

表1 RMPCEN基因的生物學過程富集分析

表2 候選預(yù)后標志物及其構(gòu)建的預(yù)后預(yù)測模型

同時考慮模型預(yù)測能力和模型基因個數(shù)，最終選擇了前3個基因構(gòu)建的預(yù)后模型，這3個基因分別為SLAMF6(HR=0.84;95%CI=0.71～0.99;logrankP=0.0416),SLAMF1(HR=0.63;95%CI=0.46～0.87;logrankP=0.005),CD2(HR=0.89; 95%CI=0.80～0.98; logrankP=0.0226)。預(yù)后模型由這3個預(yù)后標志物的表達值和Cox回歸分析得出的回歸系數(shù)線性組合構(gòu)建：Riskscore=(-0.176 739 9×SLAMF6表達量)+(-0.459 172×SLAMF1表達量)+(-0.119 112 9×CD2表達量)。生存分析表明，高風險組(n=170)的患者生存率顯著低于低風險組(n=169，logrankP=0.0055，見圖5)，這3個標志物能夠較好地區(qū)分患者的預(yù)后情況。此外，模型的 ROC曲線顯示，預(yù)后模型的AUC為0.689，表明了較好的模型預(yù)測性能。

圖5 預(yù)后預(yù)測模型的KM生存曲線和ROC曲線

3 討論與結(jié)論

根據(jù)TCGA下載的數(shù)據(jù)構(gòu)建了卵巢癌預(yù)后基因共表達網(wǎng)絡(luò)，識別了網(wǎng)絡(luò)的模塊并利用模塊基因的PPI數(shù)據(jù)對預(yù)后模塊進行重新加權(quán)。得到的RMPCEN的小世界指數(shù)為2.173，符合小世界網(wǎng)絡(luò)特性。這意味著，一方面，RMPCEN中的節(jié)點高度聚集：當一個節(jié)點連接到另外兩個節(jié)點時，后兩個節(jié)點也傾向于彼此直接連接；另一方面，網(wǎng)絡(luò)中的平均最短路徑長度幾乎與隨機網(wǎng)絡(luò)一樣低，從而說明RMPCEN具有高效的信息傳遞能力。因此，對RMPCEN中高度聚集的節(jié)點進行挖掘，通過score得分衡量節(jié)點的拓撲重要性。在score得分前10的基因中，TBX21[24]、CXCR6[25-26]和TIGIT[27]在卵巢癌細胞增殖和腫瘤生長中具有重要功能，表明拓撲學上重要的基因往往在疾病中發(fā)揮關(guān)鍵作用。目前，其余7個基因在卵巢癌中的功能尚未得到證實，但已被證明與癌癥相關(guān)。ITK在大多數(shù)轉(zhuǎn)移性黑色素瘤中異常表達，表明ITK抑制劑可能對黑素瘤治療有效[28]。IL21R在乳腺癌細胞的增殖、遷移和侵襲中發(fā)揮作用[29]。因此，這些基因可能是卵巢癌的潛在疾病基因，值得進一步研究。

同時，對RMPCEN中的基因進行了GO富集分析，預(yù)測了與卵巢癌預(yù)后相關(guān)的潛在生物過程。生物學過程富集分析表明，RMPCEN中的基因集中參與了免疫系統(tǒng)過程、免疫反應(yīng)、調(diào)節(jié)免疫系統(tǒng)過程等與免疫相關(guān)的生物學過程，這與先前的研究成果相符。目前，越來越多的證據(jù)表明，卵巢癌本質(zhì)上是免疫原性腫瘤。流行病學和臨床數(shù)據(jù)證明，卵巢癌患者的生存期與自發(fā)抗腫瘤免疫反應(yīng)、腫瘤免疫逃逸機制相關(guān)。試驗數(shù)據(jù)也證實了免疫療法的功效，這將為卵巢癌免疫療法提供新思路[30]。

本文在全基因組水平上對卵巢癌基因表達數(shù)據(jù)進行了系統(tǒng)的分析，識別了具有預(yù)后價值的生物標志物。以上討論表明，通過整合PPI數(shù)據(jù)進行基因共表達網(wǎng)絡(luò)分析有助于選擇具有生物學意義的預(yù)后標志物。此外，拓撲上重要的基因在疾病的發(fā)展中表現(xiàn)出一致的表達變異，在疾病的發(fā)展中發(fā)揮關(guān)鍵作用。通過此方法挑選的預(yù)后生物標志物具有一定的可信度，可作為卵巢癌治療的潛在靶點。