基于加權基因共表達網絡分析建立結腸腺癌預后風險模型以及關鍵基因藥物敏感性分析

黎蘊琪，蔣利和，陳闖*（. 廣西醫科大學附屬腫瘤醫院，南寧 53002；2. 右江民族醫學院，廣西 百色533000；3. 廣西大學，南寧 530004）

結直腸癌（colorectal cancer，CRC）是世界上最常見的消化系統惡性腫瘤之一。據統計，在全球范圍內，CRC發病率排第三，病死率排第二[1]。結腸腺癌（colorectal adenocarcinoma，COAD）是CRC中最常見的病理分型。目前，根治性手術治療是早期COAD最主要的治療方法，而對中晚期患者來說，化療是最常用且有效的選擇，可限制腫瘤細胞的擴散和轉移，進而控制病情。但化療藥物耐藥性的產生和腫瘤的高復發率仍是COAD臨床實踐中的主要問題[2-3]。

近年來，生物信息學在挖掘癌癥生物標志物[4]和抗腫瘤藥物分子靶點[5]等研究方面發揮重要作用。加權基因共表達網絡分析（WGCNA）主要運用于分析大樣本基因表達數據，將基因根據表達模式相似性分為不同模塊[6]，這種方法跟普通聚類方法的不同在于其將基因表達值的相關系數取了N次冪，從而使得相關系數分布更加符合無尺度網絡分析，這更符合生物學規律[7]。相比于只關注差異表達的基因，WGCNA充分利用信息，把基因與表型的關聯轉換為基因集與表型的關聯，免去了多重假設檢驗校正的問題。本研究基于WGCNA結合差異表達基因（DEGs）的方法，構建高性能的預后預測模型，同時對這些靶點進行化療藥物敏感性分析，旨在篩選出高風險腫瘤患者并對其進行精準的靶向治療，以達到個體COAD治療療效最大化。

1 材料和方法

1.1 數據采集

微陣列數據集GSE39582和GSE41258以及這些數據集的相應臨床數據從GEO數據庫（https：//www.ncbi.nlm.nih.gov/geo）中下載，通過R軟件進行數據歸一化處理，利用limma包去除批次效應，通過ggord包繪制PCA圖。GEO數據用于篩選DEGs，構建預后模型。從癌癥基因組圖譜（TCGA）數據庫（http：//cancergenome.nih.gov/）下載COAD的RNA-seq數據和相應臨床信息。該數據集包括424個COAD樣本，用于驗證預后模型和生存分析。

1.2 差異表達基因分析

采用limma軟件包尋找mRNA的差異表達，并通過pheatmap包繪制熱圖，“AdjustedP＜0.05且|log2FC|＞1”被定義為閾值，得到的DEGs被用于后續分析。

1.3 WGCNA

在評估GSE39582和GSE41258的基因表達譜后，使用WGCNA包構建無標度的共表達網絡，計算基因的皮爾遜相關系數以獲得相關矩陣，選擇適當的軟閾值來測量基因之間的連接性，將鄰接矩陣轉換為拓撲重疊矩陣（TOM），并計算相應的差異（1-TOM）。具有相似表達模式的基因被分為同一個共表達模塊。此外，我們計算了模塊特征基因的差異性，對模塊進行了層次聚類，并合并了相似的模塊（abline＝0.25）。

1.4 預后風險模型的構建

將DEGs和關鍵模塊基因取交集，得到共同基因用于進一步分析。將這些基因的表達譜和臨床信息合并，再隨機分為訓練集（70%）和測試集（30%），使用訓練集建立模型并在測試集進行驗證。在訓練集中，通過單變量Cox比例風險回歸分析篩選與預后相關的基因，納入Lasso回歸分析，防止模型過度擬合，然后納入多因素Cox回歸分析，得出COAD預后風險模型，計算每個患者的風險評分，依據訓練集風險評分的中位值，分別將訓練集和測試集分為高風險組和低風險組。運用Kaplan-Meier生存分析和ROC曲線評估該預后風險模型的準確性。利用生存時間和基因風險模型分別繪制散點圖和高低風險熱圖，利用風險評分結合臨床信息繪制列線圖。

1.5 基因集富集分析（GSEA）

為了進一步確認預后關鍵基因的潛在功能和識別高、低風險組的相關通路，使用GSEA軟件（4.0.1）進行基因組富集分析，以P＜0.05為閾值，展示KEGG相關通路的富集結果。

1.6 藥物敏感相關性分析

CellMiner是以美國國家癌癥研究所癌癥研究中心（NCI）所列出的60種癌細胞為基礎建立的數據庫[8]。從CellMiner數據庫下載NCI-60中基因mRNA表達數據和藥物活性相關數據，探究5個關鍵基因的表達與藥物敏感性之間的關系。此外，通過數據庫GSE81005分析5-氟尿嘧啶（5-FU）耐藥的結直腸癌細胞系和野生型結直腸癌細胞系中HSD17B2mRNA的表達高低。

1.7 統計學分析

使用Wilcoxon檢驗分析連續變量。使用Kaplan-Meier 生存分析和 Cox 單因素分析進行預后分析。Pearson相關性分析用于評估5個關鍵基因的表達與藥物敏感性之間的相關性。使用R軟件（3.6.3版本）進行統計分析，P＜0.05為差異有統計學意義。

2 結果

2.1 數據處理

GSE39582和GSE41258數據集共包含73個正常樣本和752個COAD組織樣本。從這些樣本中鑒定出490個DEGs，其中上調基因313個，下調基因177個（見圖1A、B）。數據批次效應情況通過批次去除前后的PCA圖進行評估（見圖1C、D）。

圖1 基因表達差異分析結果Fig 1 Gene expression differential analysis

2.2 WGCNA和關鍵模塊的識別

為了找出與COAD相關的關鍵基因，筆者使用正常樣本和COAD樣本的基因表達矩陣進行了WGCNA分析。將軟閾值β設置為5，平均連通性接近0（見圖2C、D）。將具有相似表達模式的DEGs聚集到相同的模塊中，最終得到11個共表達模塊（見圖2B），灰色模塊被認為是無法被分配給任何模塊的基因集合。綠松石色模塊與正常樣本表型相關系數為0.64，綠色模塊與腫瘤樣本表型相關系數為0.53，通過分析基因與模塊的相關性（MM）和基因與臨床特征（正常和腫瘤）的相關性（GS）之間的關聯，發現綠色和綠松石色模塊中MM和GS成正相關（見圖2E、F），說明這些與性狀高度相關的基因，在關鍵模塊中也扮演著舉足輕重的角色。因此選擇綠色和綠松色為目的模塊，并將DEGs和目的模塊基因取交集篩選出247個基因用于進一步分析。

圖2 WGCNA分析Fig 2 WGCNA analysis

2.3 預后模型的構建

將表達和生存數據合并后的740個GEO數據集樣本分為訓練集（70%）和測試集（30%）。首先使用訓練集的生存數據對247個DEGs進行單因素Cox比例風險回歸分析，鑒定出7個對預后有顯著影響的基因（P＜0.01），分別為GZMB、HSD17B2、PDK4、PIGR、PXMP2、GINS1、EMP1。將上述基因納入Lasso回歸分析降維（見圖3A、B），進一步使用多變量Cox比例風險回歸分析，最終得到5個與COAD預后相關的風險基因（見圖3C）。基于這5個關鍵基因進行風險預后模型的構建：Risk score＝-0.162GZMB+0.210HSD17B2+0.184PDK4-0.141PIGR-0.516PXMP2。生存分析結果均顯示低風險組的患者生存狀況明顯優于高風險組（P＜0.001，見圖4A、B）。訓練集1年、2年和3年的隨訪預測AUC分別為0.690、0.709、0.699（見圖4C），測試集的隨訪預測AUC，見圖4D。從生存時間和風險評分繪制的熱圖和散點圖中看出，隨著風險得分的增加，死亡的患者也增加，存活時間相對減少，由此可見模型有相對較好的預測能力（見圖5）。

圖3 LASSO及多變量Cox回歸分析進一步篩選預后基因Fig 3 Further screening of prognostic genes by the LASSO and multivariate Cox regression analysis

圖4 訓練集（A、C）、測試集（B、D）Kaplan-Meier曲線和時間依賴性ROC曲線Fig 4 Distribution of time-dependent ROC curves and Kaplan-Meier curves of the training set（A and C）and the test set（B and D）

圖5 訓練集（A）、測試集（B）風險評分的分布、患者的生存狀態和5個基因的表達熱圖Fig 5 Distribution of risk scores，patient survival status，and the five-gene expression heat map of the training set（A）and the test set（B）

2.4 預后風險模型的臨床應用

進一步通過單因素和多因素Cox回歸分析評估COAD患者風險評分和臨床病理信息與總生存期的關系。訓練集和測試集的結果均顯示風險評分可以作為獨立危險因素來預測COAD患者的生存預后（P＜0.01），多因素Cox回歸分析結果顯示，風險評分在訓練集（見圖6B）和測試集（見圖6D）中的HR分別為1.581和1.616。

圖6 訓練集（A、B）、測試集（C、D）關于臨床特征因素和風險評分的獨立預后分析Fig 6 Independent prognostic analysis of clinical characteristics and risk score of both the training set（A and B）and the test set（C and D）

2.5 模型驗證

接下來，在TCGA驗證集（424例樣本）中使用上述模型公式計算每位患者的風險評分，將患者分為高風險和低風險組，并通過KM曲線比較預后。結果表明，該模型可以區分患者的預后（P＝0.013，見圖7）。

圖7 TCGA驗證集的生存分析Fig 7 Survival on TCGA validation set

2.6 預后列線圖的建立

筆者使用臨床病理特征（年齡、性別、腫瘤分級、T分期、N分期）和5個基因的表達在GEO測試集具有完整臨床資料的503例樣本中建立了建立列線圖來預測患者的生存率。結果顯示列線圖可以預測1年、2年和3年的生存率，接近理想模型（見圖8A）。校準圖顯示出較佳的預測準確性（見圖8B），并且預測的生存率與實際生存率基本重合。

圖8 預后風險模型列線圖（A）及校準曲線圖（B）Fig 8 Nomogram（A）and calibration（B）curve of the prognostic risk model

2.7 高風險和低風險組的DEGs功能富集分析

為了探索5個關鍵基因與COAD的潛在生物學機制，筆者使用GSEA在訓練集和測試集中對高、低風險組之間的DEGs進行了功能富集分析。訓練集和測試集結果均顯示低風險組的基因集在代謝以及細胞周期等相關途徑有明顯的富集作用（見圖9），提示低風險組可能通過細胞周期通路影響腫瘤細胞凋亡進而影響患者生存。

圖9 訓練集（A）、測試集（B）GSEA富集分析Fig 9 GSEA enrichment analysis of the training set（A）and the test set（B）

2.8 藥物敏感性分析

最后，筆者研究了5個關鍵基因在60種人類癌細胞系（NCI-60）中的表達水平及其與化療藥物之間的相關性。這些基因和多種腫瘤的化療藥物的耐藥性增加有關（P＜0.05），其中HSD17B2mRNA的表達增加與多種治療COAD的化療藥物的耐藥性增加有關（見圖10A～E），包括伊立替康（cor＝-0.262，P＝0.043，中晚期CRC的一線用藥[9]）、5-氟脫氧尿苷（cor＝-0.318，P＝0.013）、氟尿苷（cor＝-0.334，P＝0.009，治療結腸癌的肝轉移[10]）、依托泊苷（cor＝-0.308，P＝0.017）、雷替曲塞（cor＝-0.298，P＝0.021）。另外，通過數據庫GSE81005對5-FU耐藥的CRC細胞系和野生型CRC細胞系進行分析，發現5-FU耐藥的CRC細胞系中HSD17B2mRNA的表達明顯高于野生型CRC細胞系（P＜0.001，見圖10F）。這些結果表明HSD17B2在CRC的發生發展和耐藥性中發揮著重要生物學作用。

圖10 HSD17B2的表達與伊立替康（A）、5-氟脫氧尿苷（B）、氟尿苷（C）、依托泊苷（D）和雷替曲塞（E）的藥物敏感性的相關性分析以及HSD17B2在CRC野生型細胞及其5-FU誘導耐藥細胞系中的表達（F）Fig 10 Correlation between HSD17B2 expression and drug sensitivity of irinotecan（A），5-fluorodeoxyuridine（B），fluorouridine（C），etoposide（D）and raltitrexed（E）； expression of HSD17B2 in CRC wild type cells and its 5-FU-induced resistant cell lines（F）

3 討論

COAD是全球范圍內最常見的癌癥類型之一，惡性程度和病死率較高[1，11]。目前COAD的治療主要采取手術切除及化療等傳統治療方法，對于早期腫瘤效果較好，但對于中、晚期腫瘤只能延長患者生命，緩解臨床癥狀，不能從根本上提高患者的生活質量，患者預后仍較差。因此，尋找COAD中用于生存評估和靶向治療的生物標志物仍至關重要。

本研究通過分析兩個GEO數據集來鑒定DEGs，構建了WGCNA共表達網絡，從中找出5個預后關鍵基因（GZMB、HSD17B2、PDK4、PIGR和PXMP2），構建COAD預后預測模型，并對其進行了驗證，列線圖和ROC曲線顯示該模型具有較好的臨床實用潛力。藥物敏感性的分析顯示，這些關鍵基因與多種化療藥物的耐藥性相關。先前的研究表明，顆粒酶（GZMs）是細胞毒性T淋巴細胞和自然殺傷細胞的關鍵效應分子，在控制細胞內病原體和癌癥方面發揮著關鍵作用[12-13]。顆粒酶B（GZMB）作為抗腫瘤和抗感染因子，在多種腫瘤中顯著表達并與良好預后相關，如非小細胞肺癌、CRC、骨肉瘤等[14-16]。17-β羥基類固醇脫氫酶2型（HSD17B2）是一種參與雌激素和雄激素調節的酶[17]，其可催化睪酮和雌二醇（E2）分別轉化為雄烯二酮和雌酮（E1）[18]。而雌二醇被報道可減少CRC細胞的增殖和促進凋亡[19]。本研究構建的風險模型也提示，隨著HSD17B2的表達增加，COAD患者的預后更差，這可能與雌二醇的失活相關。Lee等[20]發現HSD17B2可以作為獨立預測因子，其過表達與直腸癌患者的不良預后相關。丙酮酸脫氫酶激酶4（PDK4）是細胞能量代謝的關鍵調節劑，可以參與m6A調節的糖酵解和ATP生成[21]。研究表明PDK4的高表達可以促進腫瘤細胞的增殖、遷移和侵襲，例如CRC、宮頸癌、肝癌、膀胱癌細胞等[21-23]。Woolbright等[23]發現，抑制PDK4表達可能會加劇順鉑誘導的細胞死亡。聚合免疫球蛋白受體（PIGR）是黏膜免疫系統的關鍵組成部分，在非小細胞肺癌、鼻咽癌、上皮性卵巢癌等的表達降低[24-27]。PIGR被證實了和消化道腫瘤密切相關。研究發現，PIGR可以啟動Yes-Dap12-Syk-Rac1/CDC42-MEK/ERK信號級聯以促進腫瘤細胞的增殖和轉化[28]，富含PIGR的細胞外囊泡可以促進肝癌細胞的侵襲性[29]。Ohkuma等[30]發現，PIGR的mRNA和蛋白水平是胰腺癌的獨立預后因素，其高表達與患者預后不良有關。PXMP2是高等真核生物中最豐富的一種同源三聚體過氧化物酶體膜蛋白（PMP），具有通道形成活性[31-32]，參與廣泛的代謝途徑，包括脂質、氨基酸和羥基酸、嘌呤和活性氧物質的轉化[33]。但PXMP2基因與人類疾病的關系尚未被充分了解。

對高風險組和低風險組進行GSEA功能富集分析，發現低風險組在細胞周期、DNA復制、RNA聚合酶和代謝相關通路顯著富集。細胞周期失調和細胞周期蛋白依賴性激酶（CDK）激活導致的持續細胞增殖是癌癥的標志[34]。目前，已經有3種選擇性CDK4/6抑制劑獲得美國食品藥品監督管理局的批準[35]，多種CDK4/6抑制劑正處于治療癌癥的臨床試驗中[36]。此外，通過藥物敏感性分析發現這5個關鍵基因與多種化療藥物相關，特別是HSD17B2與結腸癌化療藥物（伊立替康、5-氟脫氧尿苷、依托泊苷、雷替曲塞等）存在顯著負相關。此外，5-FU耐藥的細胞系中的HSD17B2表達顯著高于野生組。5-FU是CRC治療中使用最廣泛的化療藥物之一。盡管近年越來越多的化療藥物出現，但以5-FU為基礎的化療（如 XELOX 和 FOLFOX）仍然是目前中晚期 CRC 的一線化療藥物。但由于化療耐藥性，其總體反應率僅為 10%～15%[37]，化療耐藥和嚴重的不良反應已成為影響化療治療效果的主要因素[38]。為了克服化學耐藥性和減少不良反應，許多藥物在臨床試驗中進行了測試，但迄今為止尚未取得重大進展，需要更多努力尋找能夠克服5-FU耐藥性且毒性較小的新藥以改善CRC患者。本研究提示HSD17B2可能是5-FU化療耐藥的潛在生物標志物，針對HSD17B2的靶向治療或許能在提高化療敏感性的同時改善患者的預后。總之，預后模型中的5個靶基因可能通過影響細胞周期和細胞對化療藥物的敏感性，進而影響COAD的進展和預后。

綜上，本研究篩選的5個靶基因是潛在COAD的預后生物標志物，HSD17B2更是有望作為降低COAD化療藥物耐藥性的靶點。然而，本研究是基于多個公共數據庫的數據分析，具有局限性，這些關鍵基因在結腸腺癌預后以及化療耐藥中的分子機制仍需要進一步的基礎實驗來證實。

4 結論

本研究篩選了5個樞紐基因可能為COAD的潛在預后標志物，并發現HSD17B2可能是化療耐藥的潛在生物標志物，為進一步研究COAD的發病機制和治療提供潛在靶點。