WGCNA聯合LASSO-COX方法篩選甲狀腺癌預后關鍵基因及其臨床價值分析

張澍漾，郭松雪，項承，支飛虎，謝立江，趙萍

1.紹興市中醫院 浙江中醫藥大學附屬紹興中醫院普外科，浙江紹興 312000；2.浙江大學醫學院附屬第二醫院整形外科，浙江杭州 310009；3.浙江大學醫學院附屬第二醫院甲狀腺外科，浙江杭州 310009

甲狀腺癌（thyroid cancer，THCA）發病率逐年增長，是全球增長速度最快的癌癥之一，已經成為內分泌系統中最常見的惡性腫瘤。據統計，全球THCA患者已經超過5.86萬例[1]，女性的患病率是男性的3倍[2]。近30年，女性發病率上升了22倍，男性發病率上升了15倍，已經成為了威脅人民健康的十大癌癥之一[3]。

近年來隨著生物信息學的發展與應用，基因靶向治療成為精準治療腫瘤的熱點。既往研究利用基因表達綜合數據庫（gene expression omnibus，GEO）和癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數據庫確定了肺鱗癌的鐵死亡差異表達基因（differentially expressed genes，DEGs），為癌癥的治療靶點和預后預測提供指導[4]。有研究通過免疫微環境、免疫治療反應和基因突變分析等方法，構建了穩健的腎透明細胞癌癥預后預測模型[5]。6-DNA甲基化特征是THCA患者無進展生存率的獨立預后標志物[6]。Hsa-miR-139-5p是一種THCA預后標志物，參與HNRNPF介導的可變剪接，是一種與調節THCA主要信號通路和腫瘤毒性相關的新型調控機制[7]。GPX4表達的增加通過抑制鐵死亡促進THCA的發生，且可以預測較差的臨床結果[8]。SERINC2在乳頭狀THCA中可作為潛在的腫瘤驅動生物標志物[9]。目前我國已經常規開展了BRAF基因、TERT基因、RET/PTC重排和RAS基因等突變檢測，用于THCA患者的輔助診斷[10]。

本研究通過TCGA-THCA數據集篩選出THCA和正常樣本間的DEGs，隨后進行權重基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）和套索聯系COX回歸分析（least absolute shrinkage and selection operator regression COX analysis，LASSO-COX）獲得5個與預后相關的關鍵基因，然后由關鍵基因構建預后預測模型，并采用免疫組化實驗進行驗證基因的表達，再進行受試者工作特征（receiver operating characteristic，ROC）曲線分析和生存分析，最后基于基因表達譜和風險評分進行基因集富集分析（gene set enrichment analysis，GSEA）用以評估相關途徑和分子機制，為患者預后預測提供一個潛在選擇。

1 資料與方法

1.1 數據收集與處理

在TCGA官方網站上收集2011年至2015年共511例THCA患者的臨床數據和mRNA測序數據，共569份組織樣本，其中THCA組織樣本511份，甲狀腺組織正常組織58份。隨訪時間至2015年，其中死亡人數為17例。使用Perl5.24.3軟件將原始mRNA測序數據轉換成mRNA表達矩陣。

1.2 差異基因的篩選

應用Limma R軟件包對RNA-seq表達數據行歸一化處理，并以P<0.01和|log2FC|>2為差異有統計學意義的閾值，進行差異基因表達分析篩選出THCA樣本組和正常組之間的DEGs。

1.3 WGCNA

利用基因表達譜分別計算每個基因的絕對中位差（median absolute deviation，MAD），剔除了MAD最小的前50%的基因，利用R軟件WGCNA包去除離群的基因和樣本。首先，定義β（軟閾值）的功效以確保標準的無標度網絡，然后使用網絡中每對節點基因與其Pearson相關系數之間的鄰接值構造鄰接矩陣，鄰接矩陣用于計算拓撲重疊矩陣（topological overlap matrix，TOM）和相應的不相似度（1-TOM），基于TOM的差異度量，通過平均連接層次聚類構建樹狀圖，將高度相似的模塊聚集在一起，然后以0.25的高度截止合并。依據模塊劃分結果，計算模塊與性狀Pearson相關系數，選擇與臨床特性相關性最顯著的模塊。

1.4 LASSO-COX

結合患者的生存信息，首先對上述得到的關鍵模塊基因進行LASSO-COX，減少基因之間共線性的影響，防止后續構建的風險模型變量過度擬合。LASSO回歸分析通過引入懲罰系數（λ）將冗余變量的系數壓縮為0，最后剩余的系數非零的變量為最終變量。本研究中使用5折交叉驗證以確定最優懲罰系數，獲取有效基因構建最佳預后模型。再將LASSO回歸分析得到的基因進行多因素COX回歸分析，計算每個多因素回歸系數，構建風險評分方程。

1.5 免疫組織化學檢測

2022年10月至12月在紹興市中醫院收集接受手術治療的THCA患者的腫瘤組織和癌旁組織3例（女性2例，男性1例），進行免疫組化驗證基因的表達。免疫組化采用常規兩步染色法。組織通過石蠟包埋，然后切片脫蠟至水。微波加熱修復抗原。滴加PBS洗凈，滴加一抗。4℃孵育過夜，PBS沖洗干凈，滴加二抗37℃孵育30min，DAB顯色。蘇木精再染色，脫水，透明，封片。所有藥品均購自生工生物工程（上海）股份有限公司。

1.6 風險預后模型的建立與分析

根據上述LASSO-COX的結果，構建基于基因表達的風險評分方程。選擇風險評分的最優截斷值將患者分為高、低風險兩組，進行Kaplan-Meier生存分析，利用R3.6.0軟件繪制模型的ROC曲線，計算曲線下面積（area under curve，AUC），若AUC≥0.8，且P<0.05，則認為該模型具有較好的預測性能。利用R3.6.0軟件繪制模型預測預后的列線圖，以預測模型的C指數（C-index）評價模型的區分度。

1.7 根據模型評分分組進行GSEA分析

利用GSEA研究風險評分分組與京都基因和基因組百科全書（Kyotoencyclopdia of Genes and Genomes，KEGG）通路基因集的相關性。根據風險評分值的最優截斷值將患者分成高、低風險兩組，采用GSEA軟件（version 3.0），設定參數：最小基因集為5，最大基因集為5000，1000次重抽樣，以P< 0.05，FDR<0.25為差異有統計學意義條件，評估相關途徑和分子機制。

2 結果

2.1 DEGs篩選

利用P<0.01和|log2FC|>2為篩選條件，在THCA組織-正常組間篩選到DEGs共計2130個，包含865個上調基因和1265個下調基因。

2.2 WGCNA分析

通過平均連通性圖比較發現基因間聯系軟閾值為4后，設置模塊合并閾值為0.25，繪制聚類樹狀圖得到了9個關鍵模塊。然后根據模塊性狀關系熱圖，發現青綠色模塊與臨床性狀的關聯性最高，確定其為關鍵模塊，該模塊包含403個基因。

2.3 LASSO-COX

整合生存時間、生存狀態和基因表達數據，使用LASSO-COX對青綠色模塊包含的403個基因進一步篩選，設置λ值為0.023，最終獲得了5個關鍵基因LINC02550、STEAP2、ATP2C2、PLEKHG4B和SALL4，并構建模型公式。STEAP2的免疫組化結果顯示該蛋白在腫瘤組織中的表達強于在正常組織中的表達（圖1）。

圖1 STEAP2蛋白在甲狀腺正常組織及腫瘤組織中的表達（DAB染色）

2.4 生存分析

生存分析表明，該風險評分顯著影響THCA患者的預后，（P<0.001），并且高風險評分組的患者生存率較低（圖2）。采用ROC分析獲得1、3、5年曲線下面積分別為0.93、0.88、0.86，均大于0.8，說明該模型具有良好的預測性能（圖3）。

圖2 預后基因風險評分模型的生存分析

圖3 預后基因風險評分模型的ROC曲線

2.5 諾莫圖

基于以上結果利用R3.6.0軟件繪制基于5個基因的綜合風險評分和主要臨床特征的組合預測THCA患者生存的列線圖，結果顯示列線圖對THCA患者具有良好的預后預測能力（C-index=0.96，圖4）。

圖4 風險評分聯合臨床特性建立的預測患者生存諾莫圖

2.6 GSEA分析

5基因構建的預后預測模型具有良好的預測性能，因此認為該模型可能與影響預后的信號通路有關。GSEA分析獲得關鍵通路mTOR信號通路、Hedgehog信號通路、細胞自噬調節及轉化生長因子-β信號通路，且均富集在高風險評分組（圖5）。

圖5 GSEA富集信號通路

3 討論

THCA在內分泌惡性腫瘤患者中占比超過90%，且女性患病率顯著高于男性，因為雌性激素有可能是甲狀腺良惡性細胞的一種有效生長因子[11-12]。基因突變是THCA發生、發展的重要因素之一，并與其疾病診斷、預后評估及分子靶向藥物治療密切相關。早期對THCA進行風險等級評定，不僅避免了對輕度患者過度治療，而且有利于精準預測患者的預后[13]。

本研究結果表明，LINC02550、STEAP2、ATP2C2、PLEKHG4B和SALL4這5個基因與THCA預后緊密相關。關于LINC02550的研究目前還沒有報道。STEAP2作為STEAP家族的成員，是重要的體內金屬還原酶，在維持鐵穩態中發揮重要作用，對不同的癌種有不同的響應[14]。本研究中該基因在甲狀腺腫瘤組織中是下調的，免疫組化結果顯示該蛋白在甲狀腺腫瘤組織中高表達，與已報道研究發現其在乳頭狀THCA中的低表達相反，可能由于該蛋白在甲狀腺不同類型細胞的表達不同[15]。

ATP2C2基因與許多重要的生物學功能有關，如調節鈣穩態、調節語言障礙中的語音短期記憶[16-17]。既往研究發現ATP2C2可能是乳腺癌患者腫瘤微環境狀態的指標[18-19]。PLEKHG4B在細胞之間連接成熟期間參與肌動蛋白細胞骨架重塑，在甲狀腺腫瘤中PLEKHG4B表達是下調的，可能由于其下調，增加了甲狀腺腫瘤細胞間的黏附力形成時間，有益于其轉移[20]。SALL4是公認的致癌基因，該基因不僅通過Wnt/β-catenin和PTEN/PI3K/AKT通路調節細胞生長，還通過細胞凋亡相關的基因，如CYC3和CUL3的表達影響細胞凋亡[21-22]。SALL4在不同的癌種中起不同的作用，在乳腺癌和肝癌中，通過激活Wnt/β-catenin信號，促進癌細胞的增殖、遷移和侵襲；在胃癌中激活轉化生長因子-β/SMAD信號通路，促進癌細胞轉移[23-25]。研究表明，SALL4基因在THCA組織中表達上調，可能參與腫瘤細胞的發展和轉移[26]。

GSEA分析表明mTOR信號通路、Hedgehog信號通路、細胞自噬調節、轉化生長因子-β信號通路富集在高風險評分組。Hedgehog信號通路作用于哺乳類動物的胚胎發育和組織分化，控制細胞的生長與增殖。該信號通路異常激活與多種惡性腫瘤的發生有密切關系，參與腫瘤增殖、遷移、血管生成、上皮-間質轉化、腫瘤干細胞調控等，有研究表明在THCA患者中通過抑制Hedgehog信號通路可以激活TGF-β信號通路，從而抑制腫瘤細胞的凋亡[21]。此外，TGF-β是上皮-間質轉化過程的主要開關，可通過很多信號通路來誘導上皮-間質轉化的發生。

綜上所述，本研究通過生物信息學方法分析THCA基因表達譜數據，篩選了5個關鍵基因可為潛在預后標志物，并通過免疫組化實驗進行表達驗證。此外，建立了生存預后模型，有助于THCA的預后生存預測與基因靶向治療。