基于加權基因共表達網絡探索喉癌干性關鍵基因

昝雨欣，丁 妍

1. 湖北醫藥學院生物醫藥研究院（湖北十堰 442000） 2. 湖北醫藥學院附屬太和醫院生命科學研究所（湖北十堰 442000） 3. 胚胎干細胞研究湖北省重點實驗室（湖北十堰 442000）

喉癌是耳鼻喉科最常見的惡性腫瘤，好發于歐洲地區[1]。目前喉癌的治療手段包括手術、放化療及生物治療等，雖然有一定療效，但喉癌患者的5年生存率僅為30%～40%，且大多數患者在確診時往往處于晚期，多死于腫瘤復發或轉移[2]。因此，亟待找到一種有效可靠的方法對喉癌患者進行早期診治，以改善預后。

腫瘤干細胞是一類具有自我更新和分化能力的異質性腫瘤細胞，在腫瘤的增殖、轉移和復發中發揮重要作用[3]。通過人工智能和深度學習算法，建立一種描述腫瘤干細胞的新方法，稱為腫瘤干性指數，用于描述腫瘤細胞和干細胞之間的相似程度。研究表明，mRNAsi 可能是腫瘤患者生存、分類和疾病進展的有效指標[4]。目前對喉癌的認識主要依賴于頭頸部相關研究，其獨立研究相對缺乏，且基于生物信息學探索喉癌與mRNAsi 潛在聯系的研究較少。本研究采用加權基因共表達網絡分析（weighted gene co-expression network analysis, WGCNA），基于大規模數據集識別喉癌潛在相關基因模塊，旨在識別mRNAsi 基因潛在生物標志物或治療靶點，并驗證與喉癌免疫預后相關的新生物標志物[5]。

1 資料與方法

1.1 數據來源和處理

從癌癥基因組圖譜（The Cancer Genome Atlas,TCGA）數據庫（https://www.cancer.gov/ccg/research/genome-sequencing/tcga）中獲取11 個正常樣本和112 個喉癌樣本的FPKM 數據類型和相應臨床特征的RNA 測序數據。采用R 軟件的limma 軟件包篩選喉癌和正常組織之間的差異表達基因，錯誤發現率（false discovery rate, FDR）＜0.05 且|log（FC）|＞1 被認為是顯著的。

1.2 WGCNA網絡構建

WGCNA 是一種評估基因之間相關性并將高度相關基因分類到同一模塊的方法。首先通過差異表達分析進行處理，以過濾無關信息。拓撲重疊測度用于識別共表達基因模塊，以降低網絡對虛假連接的敏感性。通過將聚類樹切割成分支，將具有高度絕對相關性的基因聚類到相同模塊。只有超過50 個基因才會被定義為一個模塊。具有較高相關性的模塊將被合并（r＜0.25），每個模塊被分配不同顏色進行可視化。

1.3 關鍵模塊和基因篩選

為了評估模塊的顯著性，將基因與樣本性狀之間的相關性計算為基因顯著性（gene significance,GS）。模塊自身基因和基因表達譜之間的相關性被定義為模塊成員（module membership, MM）。選擇mRNAsi 和表觀遺傳調控的mRNAsi（EGERmRNAsi）作為臨床表型進一步分析。|r|＞0.3 被認為具有相關性。

1.4 風險模型構建

首先進行單因素Cox 回歸分析，以評估預后相關模塊的mRNA，再采用LASSO 回歸分析確定mRNA 系數，采用以下算法計算患者的風險得分：

風險評分=∑ni=coefi×mRNA 表達

采用Rtsne 和ggplot2 R 軟件包，在mRNAs 特征表達的前提下進行主成分分析（principal component analysis, PCA），以區分喉癌患者的高低風險隊列。

1.5 藥物敏感性分析

在CellMiner 數據庫（https://discover.nci.nih.gov/cellminer/home.do）[6]下載常用藥物敏感性和基因表達數據，與上述風險系數進行相關性分析，通過R 語言中impute、limma、ggpubr、ggplot2 包篩選具有高評分的基因和對應的藥物。

1.6 免疫相關性分析

采用R 包GSVA 對免疫細胞和功能的基因表達譜進行單樣本基因集富集分析（single sample gene set enrichment analysis，ssGSEA），采用ggpubr 包將結果可視化。同時通過腫瘤免疫功能障礙與排斥（Tumor Immune Dysfunction and Exclusion, TIDE）數據庫（http://tide.dfci.harvard.edu）下載腫瘤評分數據，采用limma、ggpubr 軟件包得到高低風險基因與腫瘤免疫逃逸結果。

1.7 功能注釋

利用clusterProfiler 包對高低風險表達組進行基因集富集分析（gene set enrichment analysis，GSEA），下載MSigDB 數據庫中的“c2.cp.kegg.v7.5.1.entrez.gmt”和“c5.go.v7.5.1.symbols.gmt”基因集，進行基因本體（gene ontology, GO）分析及京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes, KEGG）通路富集分析，將FDR ＜0.05 設置為閾值。

1.8 統計分析

采用R 4.0.3 軟件進行統計分析。分類變量采用χ2檢驗，連續變量采用秩檢驗。采用Pearson相關性分析評估關鍵基因的表達與藥物敏感性之間的相關性。以P＜0.05 為差異具有統計學意義。

2 結果

2.1 mRNAsi相關差異表達基因

為進一步探討與mRNAsi 相關的功能基因，計算正常組織和喉癌組織之間存在的顯著差異，紅色為高表達基因，綠色為低表達基因。共篩選得到5 302 個差異表達基因，并繪制火山圖，見圖1-A。其中，前14 個差異表達基因熱圖見圖1-B，KRT4、CRNIN、MAL等基因在正常組織中高表達（P＜0.05）。

圖1 喉癌的差異表達基因Figure 1. Differentially expressed genes in laryngeal cancer注：A. 火山圖；B. 前14個差異表達基因熱圖。

2.2 WGCNA篩選的重要模塊與基因

選擇上述差異表達基因作為WGCNA 分析的輸入端數據，構建一個基因共表達網絡并從中識別與mRNAsi / EREG-mRNAsi 顯著關聯的基因模塊，軟閾值（β=5，R2=0.9）用于保證無規模網絡，最終獲得30 個符合要求的模塊。然后，對相關模塊的全基因表達水平進行分析，以確定相應模塊與mRNAsi 之間的相關性，見圖2-A。在30 個模塊中，紫色模塊與喉癌樣本的mRNAsi分值呈顯著正相關（MS=0.53，P＜0.05），同時隨機選擇另一種正相關的模塊（淺黃色模塊）加以驗證。紫色模塊內共有7 個顯著差異表達基因（GS ＞0.5 且MM ＞0.7），見圖2-B。淺黃色模塊內共有4 個顯著差異表達基因（GS ＞0.3 且MM ＞0.8），見圖2-C。最后，選擇上述淺黃色模塊差異表達基因（MTHFD2、LINC01932、AL121761.1、TBX2）和紫色模塊差異表達基因（NEK2、KPNA2、ASF1B、MCM10、RFC5、FAM72A、DBF4）共11個基因作為喉癌的關鍵基因。

圖2 WGCNA篩選的重要模塊與基因Figure 2. The important modules and genes for WGCNA screening注：A. 模塊-mRNAsi（EREG-mRNAsi）關聯圖；B. 紫色模塊散點圖；C. 淺黃色模塊散點圖。

2.3 預后基因篩選與驗證

進一步繪制11 個基因在正常組織和腫瘤組織中的表達趨勢圖，發現淺黃色、紫色模塊中的候選基因在喉癌中大多上調（P＜0.05），見圖3-A。模塊內基因相關性分析結果顯示，紫色模塊協同表達關系較強的基因對共有兩對，包括FAM72A和NEK2（r=0.77，P＜0.05）、MCM10和RFC5（r=0.71，P＜0.05）。淺黃色模塊協同表達關系較強的基因對共有三對，包括MTHFD2和TBX2（r=0.74，P＜0.05）、LINC01932和TBX2（r=0.71，P＜0.05）、LINC01932和MTHFD2（r=0.70，P＜0.05），見圖3-B。對11 個基因進行Cox 回歸和LASSO 回歸分析，構建預后基因風險模型：風險分數=TBX2×0.021+MTHFD2×0.085，見圖3-C。PCA 分析結果表明，所得的風險評分有效區分了高低風險組，見圖3-D。通過化療藥物數據庫篩選出前12 個結果進行可視化，結果顯示TBX2基因與維莫非尼（Vemurafenib）（r=0.587，P＜0.001）、達拉菲尼（Dabrafenib）（r=0.439，P＜0.001）均呈正相關；MTHFD2基因與奧沙利鉑（Oxaliplatin）（r=0.368，P=0.004）呈正相關，與凡德他尼（Vandetanib）（r=-0.354，P=0.005）呈負相關，見圖3-E。

圖3 喉癌與mRNAsi相關的關鍵基因表達Figure 3. Expression of key genes related to mRNAsi in laryngeal cancer注：A. 模塊基因表達水平；B. 模塊基因的相關性；C. 單變量Cox回歸分析；D. 基于風險評分的PCA圖和t-SNE圖；E. 化療藥物篩選；*P＜0.05，***P＜0.001。

2.4 預后基因免疫微環境

探索風險評分與mRNAsi 的相關性，其與RNAss 呈正相關（r=0.35，P＜0.05），與DNAss 的相關性差異無統計學意義（r=0.073，P=0.44），見圖4-A。同時檢測風險評分與腫瘤免疫逃逸的相關性，發現高低風險組存在顯著差異。TIDE 是一種模擬腫瘤免疫逃逸機制的計算方法，以預測樣本對免疫檢查點抑制劑的反應。兩基因構建的低風險組TIDE 明顯高于高風險組（P＜0.001），見圖4-B。在低風險隊列中，漿細胞樣樹突狀細胞（plasmacytoid dendritic cells，pDCs）和輔助性T 細胞2（T helper 2 cell，Th2）的比例增加，見圖4-C。在低風險隊列中，趨化因子受體（chemokine receptor，CCR）、人類白細胞抗原（human leukocyte antigen，HLA）及T 細胞共刺激(T cell co-stimulation） 3 種免疫功能均有差異表達，見圖4-D。

圖4 喉癌特異性腫瘤微環境Figure 4. Laryngeal carcinoma specific tumor microenvironment注：A. 風險曲線與mRNAsi相關性；B. 風險評分與TIDE的關系；C、D. 低風險和高風險之間免疫相關細胞與途徑的ssGSEA富集分數比較隊列；*P＜0.05，**P＜0.01，***P＜0.001。

2.5 GO和KEGG富集分析

GO 功能富集分析結果顯示，MTHFD2主要富集于表皮細胞分化，TBX2主要富集于結締組織替代的正向調節，且二者與角化作用、角化細胞分化密切相關，見圖5-A。KEGG 信號通路富集分析結果顯示，MTHFD2主要富集于嗅覺信號轉導途徑，TBX2主要富集于同種異體移植排斥反應和蛋白酶體，見圖5-B。

圖5 MTHFD2與TBX2基因功能富集分析Figure 5. Functional enrichment analysis of MTHFD2 and TBX2 genes注：A. GO富集分析；B. KEGG富集分析。

3 討論

喉癌的發病率較低，每10 萬人中僅有1.5 例發病，臨床較少見[1]。隨著人口老齡化，喉癌的發病率逐漸上升，喉癌患者在生物學行為、治療反應和預后方面的表現具有較大異質性。年齡、腫瘤分級和TNM 分類等常見的預后因素無法準確預測喉癌患者的預后，因而聚焦該疾病具有重要的臨床實際意義。

腫瘤干細胞是惡性腫瘤無限增殖和復發的來源。然而，早期研究缺乏對癌癥干細胞的統一描述，喉癌的個體化治療進展緩慢。而mRNAsi 作為量化腫瘤患者腫瘤干性的指標被開發出來，可以更加高效便捷地探索與腫瘤干性相關的基因。腫瘤干細胞在喉癌中的研究仍處于初級階段。Prince 等曾報道在原發性喉癌中檢測到CD44+癌細胞，雖然其占比不到10%，但在體外具有很強的成瘤能力[7]。研究表明，CD133 是喉癌干細胞的標志物之一，可以為靶向治療提供證據[8]。這表明腫瘤干細胞標記物與喉癌治療相結合，可能對喉癌患者的長期預后產生積極影響，并為癌癥治療帶來新思路。本研究構建了一個基于WGCNA 的mRNAsi 的喉癌預測模型，得到了關鍵基因MTHFD2和TBX2，為喉癌的機制研究提供了新見解。

mRNAsi 相關基因的藥敏數據將有助于喉癌患者的治療。與MTHFD2相關的前兩種化療藥物維莫非尼和達拉菲尼均為有效的絲氨酸/蘇氨酸蛋白激酶選擇性抑制劑，且前者在甲狀腺乳頭狀癌患者中顯示出抗腫瘤活性，后者可抑制成釉細胞癌與甲狀腺癌進展[9-10]。奧沙利鉑和凡德他尼是兩種與TBX2相關的化療藥物，奧沙利鉑是一種DNA 合成抑制劑，可以誘導免疫原細胞死亡，并對喉癌細胞具有抑制作用[11]；凡德他尼是一種多通道腫瘤信號傳導抑制劑，通過調節細胞依賴性DNA 修復和腫瘤微環境，使頭頸鱗癌對光動力學治療過敏，從而導致體內腫瘤生長的顯著增加[12]。雙基因共相關的達沙替尼是一種具有抗實體腫瘤雙重SRC/ABL 激酶抑制劑，抑制SRC 磷酸化反應誘導細胞周期停止和凋亡，從而抑制腫瘤進展[13]。這些藥物是喉癌患者的潛在新治療選擇。

腫瘤微環境的免疫細胞數量和比例隨著腫瘤進展而變化，在某些情況下可以對抗癌治療產生反應。豐富的抗腫瘤浸潤免疫細胞（pDCs 和Th2）富集在低風險隊列中，pDCs 和Th2 細胞均是參與免疫反應并在抗腫瘤反應中發揮積極作用的主要細胞。在扁桃體上皮中，pDCs 不僅可以通過識別病毒并釋放干擾素幫助抵御病毒感染，而且能夠在不受微生物過度刺激的情況下啟動免疫。因此，pDCs可以被認為是扁桃體上皮的“守護者”，幫助維護免疫系統的穩定和健康[14]。此外，臨床證據顯示，使用pDCs 的藥物分類法識別惡性腫瘤治療敏感的患者，可有效提高治療反應率[15]。以上均提示pDCs 在喉癌進展中發揮抑制作用。在喉鱗狀細胞癌患者中，Th2 型轉錄因子和細胞因子的表達占優勢[16]。經批準的免疫檢查點抑制劑可減少免疫逃逸，從而改善癌癥患者的臨床預后。臨床研究發現一例特殊的轉移性喉癌病例，患者在PD-1 抑制劑治療失敗后，出現了縱隔腹膜后淋巴結和骨質進行性疾病，但在經過免疫檢查點抑制劑藥物治療后，其骨痛、腹膜后淋巴結和胸腔淋巴結病癥狀得到緩解，臨床狀況也得到改善，這與TIDE 評分結果一致[17]，有助于開發個性化和精確的免疫治療方案。在功能富集研究中，發現關鍵基因與角化細胞密切相關，而喉癌前病變喉角化癥是一種常見的粘膜病變，起源于造血干細胞的細胞缺陷會導致喉角化狀癌發生，這提示未來治療可通過免疫學對惡性腫瘤進行早期監測[18]。

本研究存在一定局限性，盡管采用了生物信息學分析，但實際臨床研究的樣本量較少，降低了研究的可行性。綜上所述，TBX2和MTHFD2基因可成為喉癌早期檢測、診斷和預后評估的潛在生物標志物。將WGCNA模塊與mRNAsi相結合，通過精準定位這些模塊中的關鍵基因，可以確定在癌癥進展中起關鍵作用的潛在治療靶點，這些特定基因或途徑的靶向可能為開發喉癌新的治療策略和確定藥物靶點提供參考。