基于知識學習挖掘乳腺癌與甲狀腺癌的共享功能模塊和核心基因

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

（廣東醫科大學公共衛生學院1，醫學系統生物學研究所2，廣東 東莞 523808）

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

1 資料與方法

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

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

1.3 風險基因網絡的構建 以CRGs 為初始種子基因，當任意兩個基因在HPRD 存在互作關系，那么視為它們在網絡中存在連接,網絡的節點代表著基因（包含種子基因及其一級鄰居基因），最終構建出一個無向的共享風險基因網絡。在Cytoscape 軟件（版本3.6.1）進行網絡可視化。

1.4 BC與TC 關鍵功能模塊的挖掘 利用Newman網絡分割算法[11,12]進行功能模塊的挖掘，此算法的核心是將網絡分解問題轉變成二次型優化問題，即求使目標函數Q 最大的列向量s 所對應的數值：

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

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

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

1.5 拓撲學分析及模塊核心基因識別 利用網絡直徑、特征路徑長度、連通度、無標度網絡特性等指標對生物分子網絡進行拓撲學分析。基于泊松分布，通過識別某個基因的連通度是否顯著大于一般基因進行連通度（k）的檢驗來篩選核心基因：

其中，λ=NP，λ 表示節點的期望連通度；N 為節點數目；P為節點間發生連接的概率。以Bonferroni校正后P＜0.05 為差異有統計學意義。

1.6 功能富集分析 利用“ClusterProfiler”R 軟件包[13]對關鍵功能模塊進行KEGG 通路富集分析。設定FDR 校正后P＜0.01 為統計學意義顯著。對各模塊富集到的通路經KEGG 數據庫進行分類，探討關鍵功能模塊之間的關聯。

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

2 結果

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

圖1 共享風險基因韋恩圖

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

2.3 共享致病網絡模塊的挖掘和拓撲屬性分析 共享風險基因網絡被Newman 算法分解識別出17 個有統計學意義的關鍵功能模塊。各模塊的特征路徑長度均小于6，無標度檢驗的指數參數為2～3，且連通度均符合冪律分布（KS 檢驗，P＞0.05），反映所構建的網絡符合無標度性質，見表1。

表1 各功能模塊基本特征和拓撲屬性

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

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

2.5 關鍵模塊的功能富集分析 KEGG 富集分析顯示，各模塊主要富集于癌癥通路、p53 信號通路、雌激素信號途徑，見表3。進一步通過KEGG 網站查詢挖掘各功能模塊間的相互關系發現，17 個關鍵功能模塊被劃分為免疫系統（immune system）、癌癥（cancer）、傳染病（infectious disease）功能類別，見圖2。

圖2 功能模塊與通路分類的關系

表3 各功能模塊的KEGG 通路分析結果

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

表4 各模塊兩種癌癥預后相關的核心基因

3 討論

生物學功能的實現并不是由單一基因實現的，往往需要通過多個基因的協同作用。在PPI 網絡中，基因通過協同作用實現不同的生物學功能，表現為形成多個成簇聚合的高度緊密連接的節點，從而把復雜網絡按照不同的生物學功能劃分為不同的模塊，又稱為網絡模塊化。根據這一特征，可結合目前存在的多種網絡分析方法[14]將復雜的生物分子網絡分解為具有統計學意義的關鍵功能模塊，從而深入的挖掘疾病的分子機制。

本研究把共享風險基因網絡分解為17 個功能模塊，各模塊在直徑、特征路徑長度、連通度等拓撲參數接近，且擬合優度檢驗均符合無標度網絡屬性。各模塊富集到的通路與兩種癌癥密切相關，如M17中的雌激素信號通路（hsa04915），其產物雌激素作為BC 的良好預測和預后的標記物，在腫瘤微環境中具有潛在的免疫調節作用[15]。在TC 中，雌激素促進了腫瘤細胞的增長并參與調節與TC 預后最為密切的血管生成和癌細胞轉移[16]。M16 模塊中的Hippo信號通路（hsa04390）是促進BC 細胞增殖和遷移以及腫瘤生長所必需途徑[17]。在TC 中，多個致病因子（如SNHG15、ST6GAL2）的異常表達是通過作用于Hippo 信號通路來促進腫瘤的發生發展[18,19]，反映出所識別出的每一個模塊均與兩種疾病的發病機制密切相關。對各模塊富集到的KEGG 通路進行類別劃分發現，M4 模塊參與多個功能類別，該模塊及其核心基因對于揭示兩種疾病的共同發病機制具有最重要的價值。

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

綜上所述，本研究基于基因多效性理論對兩種癌癥的構建共享功能模塊網絡進行分析，挖掘它們致病的關鍵功能模塊以及核心基因，為兩種癌癥發病關聯機制以及預后提供了新的思路。