許德華，陳曉琳，廖苑君，藍(lán)樹金，孫勝南，李 讓，饒紹奇

（廣東醫(yī)科大學(xué)公共衛(wèi)生學(xué)院1，醫(yī)學(xué)系統(tǒng)生物學(xué)研究所2，廣東 東莞 523808）

乳腺癌（breast cancer，BC）和甲狀腺癌（thyroid cancer，TC）是我國常見的女性惡性腫瘤。BC 年發(fā)病數(shù)約為30.4 萬，位于女性惡性腫瘤發(fā)生首位[1]。TC年發(fā)病數(shù)約為15.10 萬，位于女性惡性腫瘤發(fā)生第4位[1]。BC與TC 的發(fā)病率都呈逐年上升趨勢[2，3]，且根據(jù)以往研究報道[4]，TC 或BC 患者發(fā)生第2 原發(fā)BC或TC 的風(fēng)險相對更高。我國2001-2010 年調(diào)查研究發(fā)現(xiàn)[5，6]，患BC 的女性未來患TC 的可能性是普通人的2 倍，而患TC 的女性未來患BC 的幾率比普通人群高67%。研究顯示[7，8]，BC 患者一級親屬患TC的風(fēng)險更高，提示BC 和TC 之間可能存在某種家族聯(lián)系。有研究顯示[9]，BC 和TC 這兩種癌癥都受下丘腦-垂體-腺體軸調(diào)控，與性別、環(huán)境、藥物治療以及遺傳等因素有關(guān)，但根本原因尚未明確。因此，本研究基于基因多效性理論，利用轉(zhuǎn)錄組數(shù)據(jù)結(jié)合生物學(xué)網(wǎng)絡(luò)分析方法挖掘兩種疾病的共享功能模塊及核心致病基因，為揭示兩種癌癥發(fā)病的關(guān)聯(lián)機制提供新思路，現(xiàn)報道如下。

1 資料與方法

1.1 數(shù)據(jù)來源 從TCGA 數(shù)據(jù)庫（https://tcga-data.nci.nih.gov/tcga/）中分別下載BC 和TC 的RNA-seq數(shù)據(jù)和臨床資料。樣本剔除標(biāo)準(zhǔn)：①重復(fù)測量樣本；②臨床信息缺失樣本。異常表達(dá)基因剔除標(biāo)準(zhǔn)：表達(dá)值在80%的樣本中均少于0 的基因。生物分子網(wǎng)絡(luò)的構(gòu)建利用人類蛋白質(zhì)參考數(shù)據(jù)庫[10]（Human Protein Reference Database，HPRD）中的蛋白質(zhì)互作關(guān)系（protein-protein interaction，PPI）。

1.2 共享風(fēng)險基因的篩選 利用“edgeR”R 軟件包對BC 和TC 的RNA-seq 數(shù)據(jù)分別進(jìn)行差異分析，以FDR 校正后P＜0.05 且|logFC|＞1 為標(biāo)準(zhǔn)篩選差異表達(dá)基因（differentially expressed genes，DEGs）。最后，利用韋恩圖對BC 和TC 的DEGs 按共同上、下調(diào)取交集，得到兩種癌癥的共享風(fēng)險基因（common risk genes，CRGs）。

1.3 風(fēng)險基因網(wǎng)絡(luò)的構(gòu)建 以CRGs 為初始種子基因，當(dāng)任意兩個基因在HPRD 存在互作關(guān)系，那么視為它們在網(wǎng)絡(luò)中存在連接,網(wǎng)絡(luò)的節(jié)點代表著基因（包含種子基因及其一級鄰居基因），最終構(gòu)建出一個無向的共享風(fēng)險基因網(wǎng)絡(luò)。在Cytoscape 軟件（版本3.6.1）進(jìn)行網(wǎng)絡(luò)可視化。

1.4 BC與TC 關(guān)鍵功能模塊的挖掘 利用Newman網(wǎng)絡(luò)分割算法[11,12]進(jìn)行功能模塊的挖掘，此算法的核心是將網(wǎng)絡(luò)分解問題轉(zhuǎn)變成二次型優(yōu)化問題，即求使目標(biāo)函數(shù)Q 最大的列向量s 所對應(yīng)的數(shù)值：

其中B 是一個對稱矩陣，表示為

s 為一個列向量，可表示為

si或sj是指節(jié)點i 或j 所屬的子網(wǎng)；sT表示s 的倒置，表示行向量；A 為s 內(nèi)所有節(jié)點的鄰接矩陣，反映基因是否存在互作關(guān)系，當(dāng)Aij=1 表示，i 基因和j 基因互作，Aij=0 則無。m 為網(wǎng)絡(luò)中邊的總數(shù)，ki或kj為節(jié)點i 或j 的連通度。根據(jù)柯西（Cauchy）不等式，s 取B 最大特征值對應(yīng)的特征向量，使目標(biāo)函數(shù)的使Q 達(dá)到最大。根據(jù)s 取值符號方向?qū)⒕W(wǎng)絡(luò)切割成兩部分，重復(fù)這一步驟，至網(wǎng)絡(luò)無法分解為止。最后將節(jié)點數(shù)大于50 的子網(wǎng)作為模塊。以上計算過程在R 語言（version3.6.1）中進(jìn)行。

1.5 拓?fù)鋵W(xué)分析及模塊核心基因識別 利用網(wǎng)絡(luò)直徑、特征路徑長度、連通度、無標(biāo)度網(wǎng)絡(luò)特性等指標(biāo)對生物分子網(wǎng)絡(luò)進(jìn)行拓?fù)鋵W(xué)分析。基于泊松分布，通過識別某個基因的連通度是否顯著大于一般基因進(jìn)行連通度（k）的檢驗來篩選核心基因：

其中，λ=NP，λ 表示節(jié)點的期望連通度；N 為節(jié)點數(shù)目；P為節(jié)點間發(fā)生連接的概率。以Bonferroni校正后P＜0.05 為差異有統(tǒng)計學(xué)意義。

1.6 功能富集分析 利用“ClusterProfiler”R 軟件包[13]對關(guān)鍵功能模塊進(jìn)行KEGG 通路富集分析。設(shè)定FDR 校正后P＜0.01 為統(tǒng)計學(xué)意義顯著。對各模塊富集到的通路經(jīng)KEGG 數(shù)據(jù)庫進(jìn)行分類，探討關(guān)鍵功能模塊之間的關(guān)聯(lián)。

1.7 生存分析 利用“survival”R 軟件包中的Kaplan-Meier 法繪制生存曲線對核心基因進(jìn)行生存分析，以P＜0.05 為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 共享風(fēng)險基因的篩選 對TCGA 數(shù)據(jù)集進(jìn)行預(yù)處理后，得到BC 樣本1036 個（病例924 個+正常112 個）；TC 樣本500 個（病例442 個+正常58 個）。差異分析后，在BC 中有5289 個DEGs（3200 個上調(diào)，2089 個下調(diào)），TC 中2972 個DEGs（1926 上調(diào)，1046 下調(diào)）。將BC 和TC 的DEGs 按上下調(diào)取交集后得到1136 個CRGs（708 個上調(diào)，428 個下調(diào)），見圖1。

圖1 共享風(fēng)險基因韋恩圖

2.2 BC 和TC 共享風(fēng)險基因網(wǎng)絡(luò)的構(gòu)建 基于HRPD 數(shù)據(jù)庫的PPI，構(gòu)建出一個包含2856 個節(jié)點和4408 個互作對的共享風(fēng)險基因網(wǎng)絡(luò)。剔除游離節(jié)點（CYB5A、TSKS、ACACB 等），得到了包含2654個節(jié)點和4282 個互作對的最大子網(wǎng)。共享風(fēng)險基因網(wǎng)絡(luò)的網(wǎng)絡(luò)直徑為12，聚類系數(shù)是0.034，網(wǎng)絡(luò)平均路徑長度為4.972。

2.3 共享致病網(wǎng)絡(luò)模塊的挖掘和拓?fù)鋵傩苑治?共享風(fēng)險基因網(wǎng)絡(luò)被Newman 算法分解識別出17 個有統(tǒng)計學(xué)意義的關(guān)鍵功能模塊。各模塊的特征路徑長度均小于6，無標(biāo)度檢驗的指數(shù)參數(shù)為2～3，且連通度均符合冪律分布（KS 檢驗，P＞0.05），反映所構(gòu)建的網(wǎng)絡(luò)符合無標(biāo)度性質(zhì)，見表1。

表1 各功能模塊基本特征和拓?fù)鋵傩?/p>

2.4 模塊核心基因的識別 本研究識別了105 個核心基因，其中M3 模塊中的核心基因數(shù)量最多，達(dá)到24 個。通過查閱GeneCards 和PudMed 數(shù)據(jù)庫發(fā)現(xiàn)69 個基因已被證實與BC 或TC 相關(guān)，其中ESR1、RXRG、S100B 等25 個基因被證實與兩種癌癥均相關(guān)，見表2。

表2 各功能模塊的核心基因

2.5 關(guān)鍵模塊的功能富集分析 KEGG 富集分析顯示，各模塊主要富集于癌癥通路、p53 信號通路、雌激素信號途徑，見表3。進(jìn)一步通過KEGG 網(wǎng)站查詢挖掘各功能模塊間的相互關(guān)系發(fā)現(xiàn)，17 個關(guān)鍵功能模塊被劃分為免疫系統(tǒng)（immune system）、癌癥（cancer）、傳染?。╥nfectious disease）功能類別，見圖2。

圖2 功能模塊與通路分類的關(guān)系

表3 各功能模塊的KEGG 通路分析結(jié)果

2.6 生存分析 共得到16 個基因與BC 預(yù)后相關(guān)，5個與TC 預(yù)后相關(guān)。預(yù)后相關(guān)基因中有8 個未有文獻(xiàn)報道與這兩種疾病相關(guān)，見表4。

表4 各模塊兩種癌癥預(yù)后相關(guān)的核心基因

3 討論

生物學(xué)功能的實現(xiàn)并不是由單一基因?qū)崿F(xiàn)的，往往需要通過多個基因的協(xié)同作用。在PPI 網(wǎng)絡(luò)中，基因通過協(xié)同作用實現(xiàn)不同的生物學(xué)功能，表現(xiàn)為形成多個成簇聚合的高度緊密連接的節(jié)點，從而把復(fù)雜網(wǎng)絡(luò)按照不同的生物學(xué)功能劃分為不同的模塊，又稱為網(wǎng)絡(luò)模塊化。根據(jù)這一特征，可結(jié)合目前存在的多種網(wǎng)絡(luò)分析方法[14]將復(fù)雜的生物分子網(wǎng)絡(luò)分解為具有統(tǒng)計學(xué)意義的關(guān)鍵功能模塊，從而深入的挖掘疾病的分子機制。

本研究把共享風(fēng)險基因網(wǎng)絡(luò)分解為17 個功能模塊，各模塊在直徑、特征路徑長度、連通度等拓?fù)鋮?shù)接近，且擬合優(yōu)度檢驗均符合無標(biāo)度網(wǎng)絡(luò)屬性。各模塊富集到的通路與兩種癌癥密切相關(guān)，如M17中的雌激素信號通路（hsa04915），其產(chǎn)物雌激素作為BC 的良好預(yù)測和預(yù)后的標(biāo)記物，在腫瘤微環(huán)境中具有潛在的免疫調(diào)節(jié)作用[15]。在TC 中，雌激素促進(jìn)了腫瘤細(xì)胞的增長并參與調(diào)節(jié)與TC 預(yù)后最為密切的血管生成和癌細(xì)胞轉(zhuǎn)移[16]。M16 模塊中的Hippo信號通路（hsa04390）是促進(jìn)BC 細(xì)胞增殖和遷移以及腫瘤生長所必需途徑[17]。在TC 中，多個致病因子（如SNHG15、ST6GAL2）的異常表達(dá)是通過作用于Hippo 信號通路來促進(jìn)腫瘤的發(fā)生發(fā)展[18,19]，反映出所識別出的每一個模塊均與兩種疾病的發(fā)病機制密切相關(guān)。對各模塊富集到的KEGG 通路進(jìn)行類別劃分發(fā)現(xiàn)，M4 模塊參與多個功能類別，該模塊及其核心基因?qū)τ诮沂緝煞N疾病的共同發(fā)病機制具有最重要的價值。

此外，本研究所構(gòu)建的網(wǎng)絡(luò)中每一模塊符合無標(biāo)度性，即存在著少數(shù)連通度較高的核心基因，它們位于中樞位置，維系模塊網(wǎng)絡(luò)的穩(wěn)定和模塊內(nèi)基因的相互調(diào)控作用，從而維持模塊的生物學(xué)功能。當(dāng)它們發(fā)生突變或者受到干擾時，模塊對應(yīng)的生物學(xué)功能也會發(fā)生異常，從而導(dǎo)致疾病的發(fā)生發(fā)展。由此，本研究對各模塊節(jié)點進(jìn)行泊松分布檢驗，識別出105 個核心基因。經(jīng)GeneCards 和PudMed 檢索發(fā)現(xiàn)，S100B、ESR1、CDKN2A 等基因被報道與BC 和TC 的癌細(xì)胞的生長、分化、侵襲和轉(zhuǎn)移密切相關(guān)。如ESR1，其編碼雌激素受體，還作為轉(zhuǎn)錄因子去調(diào)控DNA 結(jié)合、轉(zhuǎn)錄激活以及激素結(jié)合。作為BC 研究的焦點，ESR1 的激活突變是BC 治療中獲得性內(nèi)分泌耐藥的關(guān)鍵機制[20]。同時，ESR1 的高表達(dá)與侵襲性TC 行為密切相關(guān)，還可作為預(yù)測甲狀腺外延伸、淋巴結(jié)轉(zhuǎn)移、遠(yuǎn)處轉(zhuǎn)移和高TNM分期的指標(biāo)[21]。且ESR1 在本研究中的BC 和TC 疾病組中均呈高表達(dá)，反映出該基因的多效性。同時，未見DLG2、KCNJ2、GHR 等與BC 或TC 相關(guān)基因的報道，因此對其研究有助于挖掘兩種癌發(fā)病關(guān)聯(lián)的潛在機制。本研究BC 和TC 疾病組中均低表達(dá)的DLG2 是潛在的腫瘤抑制因子，其異常表達(dá)會導(dǎo)致多種癌癥上皮細(xì)胞惡性增殖和贅生性轉(zhuǎn)化[22,23]。此外，參與了細(xì)胞中鉀通道的形成的KCNJ2，研究發(fā)現(xiàn)[24]，KCNJ2 在侵襲性胃癌細(xì)胞中呈高表達(dá)水平，而沉默KCNJ2 會顯著降低胃癌細(xì)胞的侵襲和轉(zhuǎn)移能力。Liu H 等[25]在小細(xì)胞型肺癌的研究中發(fā)現(xiàn)，KCNJ2 高表達(dá)會促進(jìn)癌細(xì)胞生長并使其對化學(xué)治療藥物不敏感，表明DLG2 在多種癌癥中異常表達(dá)，表現(xiàn)出其基因多效性，可作為癌癥新的預(yù)后標(biāo)志物和治療靶標(biāo)。此外，本研究對各模塊生存分析顯示，模塊M3 中的預(yù)后相關(guān)基因數(shù)量最多，且部分基因與兩種癌癥預(yù)后均相關(guān)，提示M3 模塊對兩種癌癥的預(yù)后有著重要的研究價值；且其預(yù)后相關(guān)的基因中CD19、ACAN、RORB 等8 個基因未曾有文獻(xiàn)報道，有望成為新的癌癥預(yù)測靶點。

綜上所述，本研究基于基因多效性理論對兩種癌癥的構(gòu)建共享功能模塊網(wǎng)絡(luò)進(jìn)行分析，挖掘它們致病的關(guān)鍵功能模塊以及核心基因，為兩種癌癥發(fā)病關(guān)聯(lián)機制以及預(yù)后提供了新的思路。