引用本文： 張元欣, 杜正貴, 李宏江. 加權基因共表達網絡分析法挖掘乳腺癌發生發展相關Hub基因. 華西醫學, 2020, 35(9): 1074-1081. doi: 10.7507/1002-0179.201904162 復制
目前，乳腺癌發病率已經與肺癌發病率并列第一，占比均為 11.6%[1]，成為威脅當代女性健康的第一殺手。然而我們對包括乳腺癌在內的腫瘤發生發展認識尚不清楚，目前被廣泛認可的理論為乳腺癌是多因素導致的一種異質性極強的疾病，例如女性年齡的增加[2]、絕經后體重增加[3]、體內激素水平高[4]、初潮較早或絕經較晚及不良的生活習慣等因素均是女性罹患乳腺癌的高危風險因素。然而這些危險因素都是與基因相互作用，最后引起腫瘤的發生發展[5-7]。因此基因的改變是機體發出的第一信號，研究與乳腺癌發生發展相關基因應該作為乳腺癌病因學、發病學等研究的起點。基因的高通量檢測方法使我們獲得大量的生物遺傳信息，加權基因共表達網絡分析法（weighted gene co-expression network analysis，WGCNA）[8]便是基因網絡分析方法的一種，特指利用基因間的表達相關性預測基因間調控關系的方法。與傳統的差異基因分析方法不同，WGCNA 首先定義一群功能相似且相關性較強的基因歸納為一個模塊，模塊可以將復雜的基因網絡分解為有限的幾個功能區域，從而簡化網絡結構[9-10]。使模塊與感興趣的臨床特征構建關聯，并從中挖掘出關聯度最高的基因，我們稱之為“Hub 基因”。有研究結果顯示，Hub 基因不僅是在模塊中連接度最高的基因，也是與臨床特征關系最為密切的基因[11-12]。因此本研究主要應用 WGCNA 篩選出乳腺癌發生發展過程中的關鍵模塊，并從中挖掘出 Hub 基因，分析 Hub 基因中是否具有乳腺癌發生發展相關的特異性，從而為乳腺癌發病機制的探索提供理論基礎。現報告如下。
1 資料與方法
1.1 研究對象
在美國國家生物技術信息中心 Gene Expression Omnibus（GEO）數據庫（
圖1
GSE5364 樣本分布情況
a. 6 種實體腫瘤的數量分布柱形圖；b. 惡性腫瘤與正常組織分布的柱形圖
1.2 研究方法
1.2.1 WGCNA 聚類形成模塊
使用 AffymetrixU133A 平臺進行組織基因微陣列分析。微陣列注釋信息用于將探針與相應的基因進行匹配，去除擁有多個基因的探針后，計算出對應多個探針的基因平均值。WGCNA 包含的 Affy 包可依賴最近鄰法分析對寡核苷酸芯片和雙色 cDNA 微陣列產生的實驗數據進行預處理，如背景校正、歸一化、估計缺失值[14]，使用對數表達量方差篩選的方法選出方差較大的前 25% 基因[15]進入下一步分析。采用相關系數[8]對腫瘤組織樣本和正常樣本分別進行層次聚類分析，剔除聚類樹狀圖上無法歸類的樣本—離群樣本。研究者們觀察到基因網絡是由少數高度連接的集散節點（關鍵節點）串聯起來所控制的系統，絕大部分的節點之間連接稀疏，其分布規律實際上是服從冪律分布的規律，即出現連接數為 k 的概率 P（k），反比于 kn［P（k）～k-n］。因此應用冪指數鄰接函數Axy= power（Sxy，β）=|Sxy|β將基因間的相似矩陣轉換為鄰接矩陣，其中參數β的選擇可以使構建的網絡更加接近基因網絡真實的分布[16-17]，利用 WGCNA 包中的函數 pick Soft Threshold 即可計算權重值，自動選擇出最佳的β值。本研究設定篩選標準R2>0.85。然后找出最佳β構建無尺度分布。應用拓撲重疊測量方法（topological overlap measure，TOM）來計算基因與基因之間的關聯程度[9，18]，將鄰接矩陣轉化為拓撲矩陣。其中具有高 TOM 的一組基因被定義為“模塊”[10-19]。本研究使用的是動態剪切樹的方法，使用 R 語言中 dynamic Tree Cut 包修剪樹枝后，合并相似的共表達基因，進一步篩選聚類即可確定模塊。拓撲重疊矩陣圖（TOM 圖）是基因與基因之間相關性的熱圖。TOM 圖使用顏色編碼來描述 TOM 相異度，其中行和列均使用 TOM 相異度進行分層聚類分析。
1.2.2 挖掘與乳腺癌發生發展相關關鍵模塊及 Hub 基因
模塊內基因的整體水平可以用特征向量基因代表。計算特征向量基因與關注表型的相關系數，便得到了模塊與臨床特征的相關度[20]。使用 DAVID（The Database for Annotation，Visualization and Integrated Discovery）（
每個基因與臨床特征的相關性稱為基因顯著性（gene significance，GS），GS 越高，代表這個基因與臨床特征的關系也越強。而模塊身份（module membership，MM）表明基因屬于哪個模塊，可以反映基因在模塊內部的連接度，連接度越高，代表與這個基因關聯的其他基因節點越多。一般來說 Hub 基因需要滿足以下條件：GS>0.2，MM>0.8，P≤0.05[11]。
1.2.3 挖掘與腫瘤合集發生發展相關 Hub 基因
將 GSE5364 微陣列中所有組織類型腫瘤歸納為腫瘤合集，使用同樣的方法，挖掘出與腫瘤合集發生發展相關 Hub 基因。
1.3 統計學方法
本研究所有分析在 R 2.3.12 軟件環境中應用 WGCNA 軟件包進行分析。模塊與各種實體腫瘤發生發展的相關性采用 Pearson 相關性分析，基因網絡采用 Cytoscape 軟件可視化展示。檢驗水準為α=0.05。
2 結果
2.1 離群值篩選及軟閾值確定
對所有樣本進行聚類樹分析發現 3 個離群樣本（圖 2a）。去除離群樣本后再次進行聚類，得到的樣本聚類圖顯示所有樣本合格（圖 2b）。對 log（k）與 log［P（k）］按照軟閾值的不同進行作圖，當軟閾值β=4 時，R2=0.85，且鄰接系數最低，此時為構建無尺度網絡的最優軟閾值，因此本研究選擇的軟閾值為 4（圖 3）。
圖2
樣本聚類樹狀圖和臨床特征的熱圖
a. 發現 3 個離群樣本；b. 去除離群樣本后再次進行樣本聚類，提示未見明顯離群值
圖3
軟閾值確定
a. 構建無尺度拓撲模型，相關系數平方需要至少大于 0.8，滿足無尺度網絡構建要求，此相關系數越高，表明該網絡越符合無尺度網絡的分布；b. 平均連接度，表示不同的軟閾值對應的基因網絡中基因鄰接系數的均值，反映了網絡的平均連接水平，平均連接度越低，代表構建的網絡更加接近真實世界的網絡分布情況；c. 連接度與頻率的相關柱狀圖，可見基因網絡中隨著基因之間連接度較高的基因出現的頻率較低，符合文中提及的冪律分布；d. 確定最佳軟閾值，當β=4 時，其
2.2 分層聚類鑒定模塊
以β=4 進行網絡構建，WGCNA 最終聚類得到 10 個模塊。每個模塊用不同的顏色加以區分，條帶的寬度與基因的數目相匹配，灰色條帶代表未分類到有意義模塊的基因（圖 4）。每個模塊中的基因具體數目如表 1 所示。
表1
每個模塊的顏色及對應分布的基因數量（個）
圖4
聚類樹狀圖
上方樹狀圖代表聚類生成模塊的過程，下垂的“樹枝”即代表聚類的模塊，中間彩色條帶不同的顏色代表不同的基因模塊；下方彩色條帶提示每個模塊中基因與乳腺癌發生發展的關系，顏色越紅代表與乳腺癌發生的關系越呈正相關，顏色越淺則越為負相關
2.3 模塊與臨床信息的聯系
本研究納入數據中包含乳腺癌在內的 6 種實體惡性腫瘤組織，WGCNA 可分別探索模塊與每一種癌癥發生發展的相關性，也可以將 6 種實體腫瘤合并為一個合集，探討模塊與惡性腫瘤發生的相關程度。首先可以看到與乳腺癌發生發展具有正相關性的是藍色模塊（rp=0.6，P=3×10?34）、粉色模塊（rp=0.44，P=8×10?18）和黃色模塊（rp=0.18，P=8×10?4），其中藍色模塊與乳腺癌發生發展相關性最高；另一方面，棕色模塊呈現出與乳腺癌的發生發展負相關的趨勢（rp=?0.48，P=3×10?21）。其余模塊與乳腺癌發生發展相關性較弱，未見明顯聯系（圖 5a）。
圖5
篩選關鍵模塊及模塊可視化展示
a. 基因模塊與臨床信息的相關性分析，每個格子的顏色代表與臨床特征的相關性強弱，越紅代表與該種腫瘤的發生正相關越強，越綠則是與腫瘤的發生負相關越強；b. Catoscape 軟件將藍色模塊中基因可視化呈現
而與腫瘤合集發生發展相關性最高的關鍵模塊包括：藍色模塊（rp=0.57，P=4×10?30），粉色模塊（rp=0.42，P=7×10?16）和紅色模塊（rp=0.35，P=6×10?11），其中也是藍色模塊中的基因與腫瘤合集發生發展關系最為密切（圖 5a）。
2.4 GO/KEGG 功能富集分析
進一步探索候選模塊的功能意義，對上述 4 個與乳腺癌發生發展相關模塊的基因進行了 KEGG/GO 的功能富集分析，結果顯示：藍色模塊中的基因在細胞周期（KEGG ID: cfa04110）和病毒致癌通路（KEGG ID: hsa05203）、癌癥通路（KEGG ID: hsa05200）和系統性紅斑狼瘡（KEGG ID: hsa05322）等通路中最顯著地富集，參與了細胞周期（GO: 0022402）、有絲分裂核分裂（GO: 1903047）、有絲分裂細胞質分裂（GO: 0000910）和微管運動（GO: 0007017）等生物學過程。粉色模塊主要富集在細胞外基質（extracellular matrix，ECM）受體相互作用（KEGG ID: hsa04512）、病灶粘連（KEGG ID: hsa04510）、蛋白的溶解與吸收（KEGG ID: hsa04974）、磷酯酰肌醇 3 激酶-絲氨酸/蘇氨酸激酶信號通路（KEGG ID: hsa04151）等通路。黃色模塊主要參與了免疫反應（GO: 0006955）、調節白細胞活動（GO: 0002694）等生物過程，在細菌感染（KEGG ID: hsa05150）和風濕性關節炎（KEGG ID: hsa05323）相關通路富集。而與腫瘤發生發展負相關的棕色模塊，主要富集在甲狀腺激素合成通路（KEGG ID: hsa04918）、脂肪酸代謝通路（KEGG ID: hsa01212）、調控干細胞多能性的信號通路（KEGG ID: hsa04550）和酪氨酸代謝通路（KEGG ID: hsa00350）等。
由模塊功能富集結果可知，藍色模塊中基因的功能及參與通路與癌癥相關性最高，甚至直接參與了腫瘤形成的信號通路、病毒致腫瘤形成通路，在細胞周期、增殖、分裂等生物過程中扮演重要功能，因此將藍色模塊作為我們重點研究的目標關鍵模塊。最后使用 Catoscape 軟件將藍色模塊中所有基因可視化呈現出來（圖 5b）。
2.5 挖掘與乳腺癌發生發展相關 Hub 基因
如前所述，4 個模塊中藍色模塊為與乳腺癌發生發展相關度最高的模塊，因此從藍色模塊中篩選出 32 個 Hub 基因，為了鎖定最有意義的基因，按照 GS 降序，篩選出前 8 個與乳腺癌發生發展相關 Hub 基因，分別是NUSAP1、FOXM1、KIF20A、BIRC5、TOP2A、RRM2、CEP55、ASPM。
采用同樣的方法對腫瘤合集發生發展相關 Hub 基因進行篩選，結果提示，與腫瘤發生相關 Hub 基因包括TOP2A、MCM4、KIF20A、CDK1、NEK2、PRC1、NUSAP1、ASPM、CEP55。其中 5 個基因（NUSAP1、KIF20A、TOP2A、CEP55、ASPM）也是乳腺癌發生發展最關鍵的基因。
2.6 Hub 基因的數據庫驗證
為了驗證篩選出的 Hub 基因的準確性，單獨在數據庫中對 8 個 Hub 基因的表達情況進行驗證。熱圖顯示，Hub 基因在乳腺癌組織與正常組織間表達存在顯著的差異，8 個基因在乳腺癌腫瘤組織中均呈現高表達（紅色），而在正常組織中呈低表達或者不表達狀態（綠色）（圖 6a）；同樣的，在腫瘤合集中觀察這 8 個 Hub 基因的表達情況，仍然能對絕大部分腫瘤組織和正常組織進行聚類區分（圖 6b），但較乳腺癌組織和正常組織對比特異性稍低。
圖6
8 個 Hub 基因表達熱圖
a. 8 個 Hub 基因（關鍵基因）在乳腺癌組織和正常組織中的表達對比情況；b. 8 個 Hub 基因（關鍵基因）在腫瘤組織合集和正常組織中的表達對比情況。
3 討論
本研究使用了 WGCNA 對乳腺癌發生發展相關模塊進行了篩選，結果發現與乳腺癌發生發展正相關模塊包括藍色、粉色、黃色，負相關模塊有棕色模塊。對這 4 個模塊進行 KEGG/GO 功能富集分析的結果提示，粉色模塊主要參與 ECM 的重構、蛋白質的溶解與吸收、磷酯酰肌醇 3 激酶-絲氨酸/蘇氨酸激酶信號通路等，與 ECM 相互作用、病灶粘連等反應密切相關。而 ECM 影響腫瘤微環境的主要成分，大量研究發現 ECM 與腫瘤細胞的增殖、侵襲等能力密切相關[23]。完整的 ECM 可以發揮限制腫瘤細胞增殖、分化和轉移等功能；但隨著腫瘤等進展，腫瘤釋放大量的細胞因子進入 ECM，同時 ECM 自身也分泌多種細胞因子，導致 ECM 重塑，這時的 ECM 往往為腫瘤細胞提供一個舒適的生長環境，不僅促進腫瘤的分化、增殖，并且利于腫瘤的侵襲及遠處轉移[24]。
正常情況下，機體的監控機制發現異常突變的細胞時會第一時間將其進行清除，然而，惡性腫瘤細胞在面對免疫反應增強時，不僅可以繼續頑固生長，而且發展出抗腫瘤免疫的功能—免疫逃逸[25]。既往研究提示腫瘤患者體內往往伴隨著免疫相關物質的高表達[26]，結合本研究結果，黃色模塊中基因主要參與炎癥反應、免疫相關通路并且與乳腺癌發生發展呈現正相關趨勢，提示腫瘤的發生發展過程中伴隨著機體持續抵抗的免疫反應。可推測，腫瘤細胞可能“綁架”了機體某些免疫細胞或者免疫反應過程，使其不僅不能對腫瘤產生殺傷作用，并且為腫瘤所用，營造出適合腫瘤物質吸收代謝、促進增殖等的環境。
乳腺與甲狀腺都是由下丘腦、垂體軸調控的內分泌器官或內分泌作用靶器官。臨床上隱約觀察到內分泌功能的變化與乳腺癌發生及預后間有一定關聯，然而直至今日，乳腺癌與甲狀腺疾病的關系仍存在巨大爭議。許多臨床型研究報道了高水平的甲狀腺激素水平與乳腺癌發生相關，并且激素水平越高，預后越差[27-28]。有研究者推測其可能的機制是甲狀腺激素可與乳腺癌細胞表面 S1、S2 位點結合，結合部位是質膜蛋白整合素 αVβ3，通過與 S1 結合可以激活磷酯酰肌醇 3-激酶通路刺激缺氧誘導因子 1α 的表達，從而增強乳腺癌侵襲轉移能力[29]；或者通過與 S2 結合可以激活絲裂原活化蛋白激酶，誘導雌激素受體 α 刺激雌激素反應元件表達，在轉錄水平改變促進乳腺細胞增殖和存活[30]。相反的是，Ferreira 等[31]在乳腺癌小鼠模型中構建甲狀腺功能亢進（甲亢）組和甲狀腺功能減退組以及健康對照組，結果表明，甲亢組小鼠的乳腺癌細胞增殖明顯減少。Gago-Dominguez 等[32]也發現，甲亢或體內高含量碘可降低乳腺癌相關風險，甲亢患者機體內通過產生氧化應激反應，從而促進腫瘤細胞凋亡。本研究發現棕色模塊中基因主要參與甲狀腺激素合成、脂肪酸的代謝等通路，并且與乳腺癌發生發展呈現負相關，因此從基因網絡的層面推測體內高甲狀腺激素可能對乳腺癌發生發展是一種保護性因素，可能可以降低乳腺癌的相關風險。
藍色模塊內的基因在核的有絲分裂、細胞質的分裂及 DNA 的復制等生物學過程中明顯富集，參與了細胞周期的大部分過程。在與臨床特征相關聯的分析當中我們發現，藍色模塊是與乳腺癌發生發展最為相關的模塊，同時藍色模塊也與腫瘤合集的發生發展呈現較強相關性，尤其在肝癌、肺癌、食管癌的發生發展中，相關性指數較高。為了驗證挖掘出的 Hub 基因是否具有乳腺癌特異性，本研究同時也關注了 6 種常見的實體腫瘤構成的腫瘤合集發生發展的 Hub 基因，結果發現兩者之間存在高度的重復性，由此可推測，乳腺癌發生發展的 Hub 基因并不存在乳腺癌特異性，不同組織類型的腫瘤發生發展最為關鍵的過程可能就是細胞分裂的過程，并且這個過程在多種實體腫瘤中可能存在相似的促發通路。
綜上，WGCNA 能夠將繁復的基因聚類成為模塊從而簡化網絡的結構，篩選出的 Hub 基因通過數據庫的驗證，在正常組織和乳腺癌組織中表達確實存在明顯的差異，在功能富集中可看到，Hub 參與了有絲分裂中核分裂、細胞質分裂及 DNA 的復制等生物學過程，通過促進細胞周期進程從而促進腫瘤的發生發展。本研究為后續乳腺癌及其他惡性腫瘤的發生發展機制的探索提供了豐富的基礎證據。后期我們將設計實驗對 8 個關鍵基因促進乳腺癌發生發展進行驗證。
目前，乳腺癌發病率已經與肺癌發病率并列第一，占比均為 11.6%[1]，成為威脅當代女性健康的第一殺手。然而我們對包括乳腺癌在內的腫瘤發生發展認識尚不清楚，目前被廣泛認可的理論為乳腺癌是多因素導致的一種異質性極強的疾病，例如女性年齡的增加[2]、絕經后體重增加[3]、體內激素水平高[4]、初潮較早或絕經較晚及不良的生活習慣等因素均是女性罹患乳腺癌的高危風險因素。然而這些危險因素都是與基因相互作用，最后引起腫瘤的發生發展[5-7]。因此基因的改變是機體發出的第一信號，研究與乳腺癌發生發展相關基因應該作為乳腺癌病因學、發病學等研究的起點。基因的高通量檢測方法使我們獲得大量的生物遺傳信息，加權基因共表達網絡分析法（weighted gene co-expression network analysis，WGCNA）[8]便是基因網絡分析方法的一種，特指利用基因間的表達相關性預測基因間調控關系的方法。與傳統的差異基因分析方法不同，WGCNA 首先定義一群功能相似且相關性較強的基因歸納為一個模塊，模塊可以將復雜的基因網絡分解為有限的幾個功能區域，從而簡化網絡結構[9-10]。使模塊與感興趣的臨床特征構建關聯，并從中挖掘出關聯度最高的基因，我們稱之為“Hub 基因”。有研究結果顯示，Hub 基因不僅是在模塊中連接度最高的基因，也是與臨床特征關系最為密切的基因[11-12]。因此本研究主要應用 WGCNA 篩選出乳腺癌發生發展過程中的關鍵模塊，并從中挖掘出 Hub 基因，分析 Hub 基因中是否具有乳腺癌發生發展相關的特異性，從而為乳腺癌發病機制的探索提供理論基礎。現報告如下。
1 資料與方法
1.1 研究對象
在美國國家生物技術信息中心 Gene Expression Omnibus（GEO）數據庫（
圖1
GSE5364 樣本分布情況
a. 6 種實體腫瘤的數量分布柱形圖；b. 惡性腫瘤與正常組織分布的柱形圖
1.2 研究方法
1.2.1 WGCNA 聚類形成模塊
使用 AffymetrixU133A 平臺進行組織基因微陣列分析。微陣列注釋信息用于將探針與相應的基因進行匹配，去除擁有多個基因的探針后，計算出對應多個探針的基因平均值。WGCNA 包含的 Affy 包可依賴最近鄰法分析對寡核苷酸芯片和雙色 cDNA 微陣列產生的實驗數據進行預處理，如背景校正、歸一化、估計缺失值[14]，使用對數表達量方差篩選的方法選出方差較大的前 25% 基因[15]進入下一步分析。采用相關系數[8]對腫瘤組織樣本和正常樣本分別進行層次聚類分析，剔除聚類樹狀圖上無法歸類的樣本—離群樣本。研究者們觀察到基因網絡是由少數高度連接的集散節點（關鍵節點）串聯起來所控制的系統，絕大部分的節點之間連接稀疏，其分布規律實際上是服從冪律分布的規律，即出現連接數為 k 的概率 P（k），反比于 kn［P（k）～k-n］。因此應用冪指數鄰接函數Axy= power（Sxy，β）=|Sxy|β將基因間的相似矩陣轉換為鄰接矩陣，其中參數β的選擇可以使構建的網絡更加接近基因網絡真實的分布[16-17]，利用 WGCNA 包中的函數 pick Soft Threshold 即可計算權重值，自動選擇出最佳的β值。本研究設定篩選標準R2>0.85。然后找出最佳β構建無尺度分布。應用拓撲重疊測量方法（topological overlap measure，TOM）來計算基因與基因之間的關聯程度[9，18]，將鄰接矩陣轉化為拓撲矩陣。其中具有高 TOM 的一組基因被定義為“模塊”[10-19]。本研究使用的是動態剪切樹的方法，使用 R 語言中 dynamic Tree Cut 包修剪樹枝后，合并相似的共表達基因，進一步篩選聚類即可確定模塊。拓撲重疊矩陣圖（TOM 圖）是基因與基因之間相關性的熱圖。TOM 圖使用顏色編碼來描述 TOM 相異度，其中行和列均使用 TOM 相異度進行分層聚類分析。
1.2.2 挖掘與乳腺癌發生發展相關關鍵模塊及 Hub 基因
模塊內基因的整體水平可以用特征向量基因代表。計算特征向量基因與關注表型的相關系數，便得到了模塊與臨床特征的相關度[20]。使用 DAVID（The Database for Annotation，Visualization and Integrated Discovery）（
每個基因與臨床特征的相關性稱為基因顯著性（gene significance，GS），GS 越高，代表這個基因與臨床特征的關系也越強。而模塊身份（module membership，MM）表明基因屬于哪個模塊，可以反映基因在模塊內部的連接度，連接度越高，代表與這個基因關聯的其他基因節點越多。一般來說 Hub 基因需要滿足以下條件：GS>0.2，MM>0.8，P≤0.05[11]。
1.2.3 挖掘與腫瘤合集發生發展相關 Hub 基因
將 GSE5364 微陣列中所有組織類型腫瘤歸納為腫瘤合集，使用同樣的方法，挖掘出與腫瘤合集發生發展相關 Hub 基因。
1.3 統計學方法
本研究所有分析在 R 2.3.12 軟件環境中應用 WGCNA 軟件包進行分析。模塊與各種實體腫瘤發生發展的相關性采用 Pearson 相關性分析，基因網絡采用 Cytoscape 軟件可視化展示。檢驗水準為α=0.05。
2 結果
2.1 離群值篩選及軟閾值確定
對所有樣本進行聚類樹分析發現 3 個離群樣本（圖 2a）。去除離群樣本后再次進行聚類，得到的樣本聚類圖顯示所有樣本合格（圖 2b）。對 log（k）與 log［P（k）］按照軟閾值的不同進行作圖，當軟閾值β=4 時，R2=0.85，且鄰接系數最低，此時為構建無尺度網絡的最優軟閾值，因此本研究選擇的軟閾值為 4（圖 3）。
圖2
樣本聚類樹狀圖和臨床特征的熱圖
a. 發現 3 個離群樣本；b. 去除離群樣本后再次進行樣本聚類，提示未見明顯離群值
圖3
軟閾值確定
a. 構建無尺度拓撲模型，相關系數平方需要至少大于 0.8，滿足無尺度網絡構建要求，此相關系數越高，表明該網絡越符合無尺度網絡的分布；b. 平均連接度，表示不同的軟閾值對應的基因網絡中基因鄰接系數的均值，反映了網絡的平均連接水平，平均連接度越低，代表構建的網絡更加接近真實世界的網絡分布情況；c. 連接度與頻率的相關柱狀圖，可見基因網絡中隨著基因之間連接度較高的基因出現的頻率較低，符合文中提及的冪律分布；d. 確定最佳軟閾值，當β=4 時，其
2.2 分層聚類鑒定模塊
以β=4 進行網絡構建，WGCNA 最終聚類得到 10 個模塊。每個模塊用不同的顏色加以區分，條帶的寬度與基因的數目相匹配，灰色條帶代表未分類到有意義模塊的基因（圖 4）。每個模塊中的基因具體數目如表 1 所示。
表1
每個模塊的顏色及對應分布的基因數量（個）
圖4
聚類樹狀圖
上方樹狀圖代表聚類生成模塊的過程，下垂的“樹枝”即代表聚類的模塊，中間彩色條帶不同的顏色代表不同的基因模塊；下方彩色條帶提示每個模塊中基因與乳腺癌發生發展的關系，顏色越紅代表與乳腺癌發生的關系越呈正相關，顏色越淺則越為負相關
2.3 模塊與臨床信息的聯系
本研究納入數據中包含乳腺癌在內的 6 種實體惡性腫瘤組織，WGCNA 可分別探索模塊與每一種癌癥發生發展的相關性，也可以將 6 種實體腫瘤合并為一個合集，探討模塊與惡性腫瘤發生的相關程度。首先可以看到與乳腺癌發生發展具有正相關性的是藍色模塊（rp=0.6，P=3×10?34）、粉色模塊（rp=0.44，P=8×10?18）和黃色模塊（rp=0.18，P=8×10?4），其中藍色模塊與乳腺癌發生發展相關性最高；另一方面，棕色模塊呈現出與乳腺癌的發生發展負相關的趨勢（rp=?0.48，P=3×10?21）。其余模塊與乳腺癌發生發展相關性較弱，未見明顯聯系（圖 5a）。
圖5
篩選關鍵模塊及模塊可視化展示
a. 基因模塊與臨床信息的相關性分析，每個格子的顏色代表與臨床特征的相關性強弱，越紅代表與該種腫瘤的發生正相關越強，越綠則是與腫瘤的發生負相關越強；b. Catoscape 軟件將藍色模塊中基因可視化呈現
而與腫瘤合集發生發展相關性最高的關鍵模塊包括：藍色模塊（rp=0.57，P=4×10?30），粉色模塊（rp=0.42，P=7×10?16）和紅色模塊（rp=0.35，P=6×10?11），其中也是藍色模塊中的基因與腫瘤合集發生發展關系最為密切（圖 5a）。
2.4 GO/KEGG 功能富集分析
進一步探索候選模塊的功能意義，對上述 4 個與乳腺癌發生發展相關模塊的基因進行了 KEGG/GO 的功能富集分析，結果顯示：藍色模塊中的基因在細胞周期（KEGG ID: cfa04110）和病毒致癌通路（KEGG ID: hsa05203）、癌癥通路（KEGG ID: hsa05200）和系統性紅斑狼瘡（KEGG ID: hsa05322）等通路中最顯著地富集，參與了細胞周期（GO: 0022402）、有絲分裂核分裂（GO: 1903047）、有絲分裂細胞質分裂（GO: 0000910）和微管運動（GO: 0007017）等生物學過程。粉色模塊主要富集在細胞外基質（extracellular matrix，ECM）受體相互作用（KEGG ID: hsa04512）、病灶粘連（KEGG ID: hsa04510）、蛋白的溶解與吸收（KEGG ID: hsa04974）、磷酯酰肌醇 3 激酶-絲氨酸/蘇氨酸激酶信號通路（KEGG ID: hsa04151）等通路。黃色模塊主要參與了免疫反應（GO: 0006955）、調節白細胞活動（GO: 0002694）等生物過程，在細菌感染（KEGG ID: hsa05150）和風濕性關節炎（KEGG ID: hsa05323）相關通路富集。而與腫瘤發生發展負相關的棕色模塊，主要富集在甲狀腺激素合成通路（KEGG ID: hsa04918）、脂肪酸代謝通路（KEGG ID: hsa01212）、調控干細胞多能性的信號通路（KEGG ID: hsa04550）和酪氨酸代謝通路（KEGG ID: hsa00350）等。
由模塊功能富集結果可知，藍色模塊中基因的功能及參與通路與癌癥相關性最高，甚至直接參與了腫瘤形成的信號通路、病毒致腫瘤形成通路，在細胞周期、增殖、分裂等生物過程中扮演重要功能，因此將藍色模塊作為我們重點研究的目標關鍵模塊。最后使用 Catoscape 軟件將藍色模塊中所有基因可視化呈現出來（圖 5b）。
2.5 挖掘與乳腺癌發生發展相關 Hub 基因
如前所述，4 個模塊中藍色模塊為與乳腺癌發生發展相關度最高的模塊，因此從藍色模塊中篩選出 32 個 Hub 基因，為了鎖定最有意義的基因，按照 GS 降序，篩選出前 8 個與乳腺癌發生發展相關 Hub 基因，分別是NUSAP1、FOXM1、KIF20A、BIRC5、TOP2A、RRM2、CEP55、ASPM。
采用同樣的方法對腫瘤合集發生發展相關 Hub 基因進行篩選，結果提示，與腫瘤發生相關 Hub 基因包括TOP2A、MCM4、KIF20A、CDK1、NEK2、PRC1、NUSAP1、ASPM、CEP55。其中 5 個基因（NUSAP1、KIF20A、TOP2A、CEP55、ASPM）也是乳腺癌發生發展最關鍵的基因。
2.6 Hub 基因的數據庫驗證
為了驗證篩選出的 Hub 基因的準確性，單獨在數據庫中對 8 個 Hub 基因的表達情況進行驗證。熱圖顯示，Hub 基因在乳腺癌組織與正常組織間表達存在顯著的差異，8 個基因在乳腺癌腫瘤組織中均呈現高表達（紅色），而在正常組織中呈低表達或者不表達狀態（綠色）（圖 6a）；同樣的，在腫瘤合集中觀察這 8 個 Hub 基因的表達情況，仍然能對絕大部分腫瘤組織和正常組織進行聚類區分（圖 6b），但較乳腺癌組織和正常組織對比特異性稍低。
圖6
8 個 Hub 基因表達熱圖
a. 8 個 Hub 基因（關鍵基因）在乳腺癌組織和正常組織中的表達對比情況；b. 8 個 Hub 基因（關鍵基因）在腫瘤組織合集和正常組織中的表達對比情況。
3 討論
本研究使用了 WGCNA 對乳腺癌發生發展相關模塊進行了篩選，結果發現與乳腺癌發生發展正相關模塊包括藍色、粉色、黃色，負相關模塊有棕色模塊。對這 4 個模塊進行 KEGG/GO 功能富集分析的結果提示，粉色模塊主要參與 ECM 的重構、蛋白質的溶解與吸收、磷酯酰肌醇 3 激酶-絲氨酸/蘇氨酸激酶信號通路等，與 ECM 相互作用、病灶粘連等反應密切相關。而 ECM 影響腫瘤微環境的主要成分，大量研究發現 ECM 與腫瘤細胞的增殖、侵襲等能力密切相關[23]。完整的 ECM 可以發揮限制腫瘤細胞增殖、分化和轉移等功能；但隨著腫瘤等進展，腫瘤釋放大量的細胞因子進入 ECM，同時 ECM 自身也分泌多種細胞因子，導致 ECM 重塑，這時的 ECM 往往為腫瘤細胞提供一個舒適的生長環境，不僅促進腫瘤的分化、增殖，并且利于腫瘤的侵襲及遠處轉移[24]。
正常情況下，機體的監控機制發現異常突變的細胞時會第一時間將其進行清除，然而，惡性腫瘤細胞在面對免疫反應增強時，不僅可以繼續頑固生長，而且發展出抗腫瘤免疫的功能—免疫逃逸[25]。既往研究提示腫瘤患者體內往往伴隨著免疫相關物質的高表達[26]，結合本研究結果，黃色模塊中基因主要參與炎癥反應、免疫相關通路并且與乳腺癌發生發展呈現正相關趨勢，提示腫瘤的發生發展過程中伴隨著機體持續抵抗的免疫反應。可推測，腫瘤細胞可能“綁架”了機體某些免疫細胞或者免疫反應過程，使其不僅不能對腫瘤產生殺傷作用，并且為腫瘤所用，營造出適合腫瘤物質吸收代謝、促進增殖等的環境。
乳腺與甲狀腺都是由下丘腦、垂體軸調控的內分泌器官或內分泌作用靶器官。臨床上隱約觀察到內分泌功能的變化與乳腺癌發生及預后間有一定關聯，然而直至今日，乳腺癌與甲狀腺疾病的關系仍存在巨大爭議。許多臨床型研究報道了高水平的甲狀腺激素水平與乳腺癌發生相關，并且激素水平越高，預后越差[27-28]。有研究者推測其可能的機制是甲狀腺激素可與乳腺癌細胞表面 S1、S2 位點結合，結合部位是質膜蛋白整合素 αVβ3，通過與 S1 結合可以激活磷酯酰肌醇 3-激酶通路刺激缺氧誘導因子 1α 的表達，從而增強乳腺癌侵襲轉移能力[29]；或者通過與 S2 結合可以激活絲裂原活化蛋白激酶，誘導雌激素受體 α 刺激雌激素反應元件表達，在轉錄水平改變促進乳腺細胞增殖和存活[30]。相反的是，Ferreira 等[31]在乳腺癌小鼠模型中構建甲狀腺功能亢進（甲亢）組和甲狀腺功能減退組以及健康對照組，結果表明，甲亢組小鼠的乳腺癌細胞增殖明顯減少。Gago-Dominguez 等[32]也發現，甲亢或體內高含量碘可降低乳腺癌相關風險，甲亢患者機體內通過產生氧化應激反應，從而促進腫瘤細胞凋亡。本研究發現棕色模塊中基因主要參與甲狀腺激素合成、脂肪酸的代謝等通路，并且與乳腺癌發生發展呈現負相關，因此從基因網絡的層面推測體內高甲狀腺激素可能對乳腺癌發生發展是一種保護性因素，可能可以降低乳腺癌的相關風險。
藍色模塊內的基因在核的有絲分裂、細胞質的分裂及 DNA 的復制等生物學過程中明顯富集，參與了細胞周期的大部分過程。在與臨床特征相關聯的分析當中我們發現，藍色模塊是與乳腺癌發生發展最為相關的模塊，同時藍色模塊也與腫瘤合集的發生發展呈現較強相關性，尤其在肝癌、肺癌、食管癌的發生發展中，相關性指數較高。為了驗證挖掘出的 Hub 基因是否具有乳腺癌特異性，本研究同時也關注了 6 種常見的實體腫瘤構成的腫瘤合集發生發展的 Hub 基因，結果發現兩者之間存在高度的重復性，由此可推測，乳腺癌發生發展的 Hub 基因并不存在乳腺癌特異性，不同組織類型的腫瘤發生發展最為關鍵的過程可能就是細胞分裂的過程，并且這個過程在多種實體腫瘤中可能存在相似的促發通路。
綜上，WGCNA 能夠將繁復的基因聚類成為模塊從而簡化網絡的結構，篩選出的 Hub 基因通過數據庫的驗證，在正常組織和乳腺癌組織中表達確實存在明顯的差異，在功能富集中可看到，Hub 參與了有絲分裂中核分裂、細胞質分裂及 DNA 的復制等生物學過程，通過促進細胞周期進程從而促進腫瘤的發生發展。本研究為后續乳腺癌及其他惡性腫瘤的發生發展機制的探索提供了豐富的基礎證據。后期我們將設計實驗對 8 個關鍵基因促進乳腺癌發生發展進行驗證。
