基于銅死亡相關基因的膠質瘤預后模型構建

張博弘，張祺

膠質瘤是最常見的原發性中樞神經系統腫瘤，約占顱內惡性腫瘤的80%，根據組織和分子特征，可以分為WHO Ⅰ~Ⅳ級。膠質母細胞瘤是惡性程度最高的膠質瘤，盡管采取了包括最大限度的手術切除、放療和化療的標準治療措施，其中位生存期仍只有大約16 個月。WHO Ⅱ~Ⅲ級的膠質瘤也具有侵襲性生長和進展到更高級別的惡性行為，大部分病人最終死于這種疾?。?，2］。探索新的預后模型，識別高?；颊?，從而精準治療，是提高膠質瘤臨床預后的有效措施，除了臨床病理特征外，個體基因特征可望構建有效的預后模型。

銅死亡（cuproptosis），是2022 年3 月Tsvetkov等在《Science》上首次提出的一種銅依賴性的、全新的細胞死亡方式，不同于已知的細胞凋亡、壞死性凋亡、焦亡和鐵死亡，它是通過銅離子與線粒體呼吸中的三羧酸循環中的脂?；煞纸Y合，使得脂?；鞍踪|異常聚集，鐵硫簇蛋白下調，導致蛋白質毒性應激并最終引起細胞死亡［3］。此外，研究者還初步鑒定出12 個銅死亡相關基因（CRGs），包括：FDX1，LIAS，LIPT1，DLD，DLAT，PDHA1，PDHB，ATP7B，MTF1，GLS，CDKN2A 和SLC31A1［3］。近半年來，在多種腫瘤中報道了基于銅死亡相關基因構建的預后模型，準確預測患者預后、腫瘤免疫微環境以及對藥物的治療反應，其中包括：腎透明細胞癌［4］、肝細胞癌［5］、黑色素瘤［6］、頭頸部鱗狀細胞癌［7］、肺癌［8］和乳腺癌［9］。然而，銅死亡相關基因在預測膠質瘤患者預后中的作用仍不清楚。此外，最近，基于新型程序性細胞死亡（焦亡、鐵死亡）相關基因構建的膠質瘤預后模型被相繼報道，但，預測效能仍不能令人滿意。銅死亡是最新報道的程序性細胞死亡方式，探索基于銅死亡相關基因的預后模型有望獲得預測效能更高的預后模型。為探討CRGs 對膠質瘤患者預后預測的價值，本研究擬通過分析公共數據庫中膠質瘤患者的轉錄組表達譜和相應的臨床數據，基于CRGs 構建并驗證預后風險評分模型，進一步構建諾莫圖提高模型臨床適用性，為改善膠質瘤患者預后評估以及精準治療提供新的思路。

1 資料與方法

1.1 數據收集

從GTEx（https：//www.gtexportal.org/）數據庫中下載正常腦組織的轉錄組表達譜數據，共納入1157 個正常腦組織樣本。從TCGA（The Cancer Genome Atlas）（https：//protal.gdc.cancer.gov/）和GEO（Gene Expression Omnibus）（https：//www.ncbi.nl m.nih.gov/geo/）數據庫中下載膠質瘤患者的轉錄組表達譜和相應的臨床數據。獲得了一個TCGA 膠質瘤隊列和一個GEO 膠質瘤隊列（GSE43378）用于后續研究。每個數據集采用“limma”包的normalize BetweenArrays 方法進行背景調整和分位數歸一化處理，通過“Combat”算法消除批次效應后合并兩個數據集。通過刪除缺乏生存數據的樣本，共納入720 個膠質瘤樣本進行研究。臨床變量包括年齡、性別、分級、總生存時間、生存狀態。

1.2 銅死亡相關基因分型

從文獻報道中提取出12 個銅死亡相關基因，根據預后相關CRGs 的表達譜，通過一致性聚類分析把膠質瘤樣本分為不同的銅死亡亞型：CRGcluster A 和CRGcluster B。一致性聚類是一種無監督聚類方法，利用“ConsensusClusterPlus”R 包完成分析，本研究使用的聚類算法是K-means 聚類算法，距離是基于歐氏距離（euclidean）。聚類的標準如下：第一，累積分布函數（CDF）曲線逐漸平緩增加；第二，沒有一個組的樣本量很??；最后，聚類后，組內相關性增加，組間相關性降低。接著，利用prcomp 函數進行主成分分析（PCA）。然后，采用R 包“survival”和“survminer”，通過Kaplan-Meier法對不同亞型進行生存分析。利用R 包“limma”，篩選CRGcluster A 和CRGcluster B 的差異表達基因（DEGs），篩選條件為FC>1.5 和adjustedP<0.05，進一步通過單因素Cox 回歸分析鑒定出1831 個預后相關的DEGs。

1.3 構建預后風險評分模型

首先，按照1∶1 比例，把膠質瘤樣本隨機分為訓練集（n=360）和測試集（n=360）。在訓練集中，基于預后相關的DEGs，使用R 包“glmnet”，通過Lasso Cox 回歸分析建立預后模型，并根據最小偏似然偏差確定預后模型。采用多因素Cox 分析選擇候選基因，得到風險評分計算公式如下：

風險評分=Σ（Expi*Coefi）

其中，Expi 和Coefi 分別表示各基因的表達量和風險系數。根據風險評分中位數將患者分為低風險組和高風險組，并對兩組患者進行Kaplan-Meier 生存分析。利用R 包“timeROC”建立時間依賴ROC 曲線以評估預后模型的預測性能。

1.4 構建諾莫圖評分系統

采用R 包“rms”，根據臨床特征和風險評分構建預測預后的諾莫圖。在諾莫圖評分系統中，每個變量都對應一個分數，通過將所有變量的分數相加得到總分。校準圖用于描述1、3 和5 年生存率的預測值與實際觀察到的結果之間的差異。

1.5 統計學分析

所有統計分析均采用R（Version 4.1.3）進行。如果未特別注明，P<0.05 被認為具有統計學意義。

2 結 果

2.1 膠質瘤中銅死亡相關基因的差異表達和預后分析

從TCGA、GEO 和GTEx 數據庫中下載膠質瘤和正常對照樣本的轉錄組數據和相應的臨床數據，共獲得1877 個樣本，其中膠質瘤樣本720 個，正常對照樣本1157個。本研究共納入12個CRGs，包括：FDX1，LIAS，LIPT1，DLD，DLAT，PDHA1，PDHB，ATP7B，MTF1，GLS，CDKN2A 和SLC31A1。通過比較膠質瘤和正常組織的CRGs 的mRNA表達水平差異，發現10 個CRGs 在膠質瘤顯著上調，2 個CRGs 顯著下調，差異都具有統計學意義（圖1）。進一步通過單因素Cox回歸和Kaplan-Meier分析發現，包括ATP7B、SLC31A1、FDX1、LIAS、LIPT1、DLD、DLAT 和PDHA1 在內的8 個CRGs 具有預后意義（表1）。提示CRGs 在膠質瘤發生、發展中具有潛在的作用。

表1 膠質瘤中8 個具有預后意義的CRGs

圖1 膠質瘤銅死亡相關基因差異表達分析 膠質瘤和正常組織的CRGs 的mRNA 表達水平，比較采用Wilcoxon rank sum test，****P<0.0001 vs.正常組

2.2 膠質瘤銅死亡亞型的鑒定

根據8 個預后相關CRGs 的表達譜，采用一致性聚類算法，對膠質瘤樣本進行分型。結果表明，K=2 似乎是最佳選擇，可以把整個隊列分成CRGcluster A 和CRGcluster B（圖2A-B），且主成分分析（principal component analysis，PCA）表明兩個亞型之間具有較好的區分度（圖2C）。進一步生存分析，Kaplan-Meier 曲線顯示CRGcluster A 的膠質瘤患者具有較長的總體生存時間（OS）（圖2D）。采用R 包“limma”比較不同銅死亡亞型的基因表達譜，得到差異表達基因（DEGs），進一步通過單因素Cox 回歸分析獲得1831 個預后相關的DEGs。

圖2 膠質瘤銅死亡亞型鑒定 A、B：一致性聚類分析K=2 時的一致性矩陣熱圖，同時一致性累積分布函數圖提示K=2 時CDF 下降坡度小，最佳的K 值為2；C：主成分分析顯示兩個銅死亡亞型之間轉錄組差異顯著；D：KM 曲線，CRGcluster A 類患者的OS 時間更長，比較采用chi-square test：χ2=38.04，P<0.001

2.3 構建預后風險評分模型

采用R 包“caret package”把樣本按1∶1 隨機分成訓練集（n=360）和測試集（n=360）。在訓練組集中，對1831 個預后相關的DEGs 進行LASSO Cox回歸分析，根據最小偏似然偏差確定了8 個基因的預后模型（圖3A~B）：風險評分=MGME1×0.4304+GPR156×0.2167+ARNTL×0.3229+RAB38×0.2424+IGFBP5 × 0.1675 + SERPINA5 × 0.1579 + FABP5 ×0.1809+COL22A1×0.1415。計算訓練集中所有患者的風險評分，根據中位風險評分，把患者分為高風險組（180 例）和低風險組（180 例），風險評分的分布圖顯示，隨著風險評分的增加，生存時間減少，死亡率增加（圖3C），Kaplan-Meier 生存曲線顯示，低風險組患者總體生存期較長（圖3D）。此外，接受者操作特征曲線（ROC）顯示，曲線下面積（AUC）分別為：0.872（1 年）、0.933（3 年）、0.889（5 年），表明該模型具有良好的預測性能（圖3E）。

圖3 在訓練集中構建預后風險評分模型 A~B：LASSO 回歸系數分布圖和10 倍交叉驗證選取最優λ 值；C：訓練集中患者風險評分、總體生存時間、生存狀態的分布；D：Kaplan-Meier 生存曲線；E：模型的ROC 曲線及AUC 值

2.4 驗證風險評分模型

計算測試集中所有樣本的風險評分，根據中位風險評分，把測試集樣本分為高風險組（183 例）和低風險組（177 例），與訓練集中結果一致，風險評分的分布圖顯示，隨著風險評分的增加，生存時間減少，死亡率增加（圖4A），Kaplan-Meier 生存曲線顯示，低風險組患者總體生存期較長（圖4B），此外，ROC 曲線顯示AUC 分別為：0.881（1 年）、0.894（3 年）、0.874（5 年），表明該模型對測試集膠質瘤患者的同樣具有良好的預測能力（圖4C）。

圖4 在測試集中驗證風險評分模型 A：測試集中患者風險評分、總體生存時間、生存狀態的分布；B：Kaplan-Meier 生存曲線；C：模型的ROC 曲線及AUC 值

2.5 構建諾莫圖預測膠質瘤患者生存率

考慮到風險評分在臨床應用上缺乏便利性，我們建立了一個包含風險評分和臨床參數的諾莫圖，用于預測膠質瘤患者1 年、3 年和5 年的生存率（圖5A），其中預測因素包括患者的性別、年齡和風險評分，同時校準圖表明，本諾莫圖具有理想模型相似的預測性能（圖5B）。

3 討 論

越來越多的證據顯示，同樣WHO 級別的膠質瘤，其生物學行為和臨床預后卻相差甚遠，有的患者能長期存活或具有較長的生存期，而有的患者卻表現為高度惡性的轉歸，因此準確的預后評估并選擇性地給予不同的治療策略是十分必要的。目前臨床應用的膠質瘤預后風險分層標準主要有美國的RTOG9802 標準［10］和歐洲的EORTC22844標準［11］。然而，這兩個標準仍存在不足的地方，有時候甚至出現較大的偏差。隨著膠質瘤基因組測序的普及，海量的生物數據提供了更多的預后信息，補充甚至模糊了傳統的WHO 分級標準。近幾年來，許多研究者借助公共數據庫，基于不同的基因表達譜，通過多種分析方法，在不同腫瘤中構建了預后模型。Chen 等［12］基于11 個鐵死亡相關基因表達譜構建了一個膠質瘤預后風險評分模型，其1 年的AUC 值是0.879。Zhang 等［13］基于10 個焦亡相關基因表達譜構建了一個膠質瘤預后風險評分模型，其1年的AUC值是0.872。Feng等［14］基于6個自噬相關基因表達譜構建了一個低級別膠質瘤預后風險評分模型，其1 年的AUC 值是0.840。

2022 年3 月，Tsvetkov 等首次在《Science》上提出了一種具有銅依賴性的、全新的細胞死亡方式-“銅死亡（cuproptosis）”，研究指出Cu2+直接與三羧酸循環中的脂?；煞纸Y合，導致脂?；鞍桩惓＞奂丸F硫簇蛋白丟失，進而引發蛋白質毒性應激反應，最終致使細胞死亡［3］，并初步列出了12 個銅死亡相關基因（CRGs）。受此啟發，基于CRGs 構建的預后模型在多個腫瘤類型相繼被報道［4-9］。為構建基于CRGs 的膠質瘤預后模型。本研究首先探討了12 個CRGs 在腫瘤和正常樣本的差異表達情況和預后意義，基于8 個預后相關CRGs 的表達譜，對膠質瘤樣本進行銅死亡分型（CRGcluster A 和CRGcluster B），然后進行銅死亡分型的差異基因分析和預后分析，得到1831 個預后相關差異基因，接著通過LASSO 回歸分析，在訓練集中構建了一個包含8 個基因的預后風險評分模型，其1 年的AUC 值是0.872，并在測試集中進行了驗證，其1 年的AUC 值是0.881，優于已報道模型的預測性能。

本研究所構建的模型包含以下8 個基因：MGME1、GPR156、ARNTL、RAB38、IGFBP5、SERPINA5、FABP5、COL22A1。MGME1（線粒體基因組維持外切酶1）對線粒體DNA（mtDNA）的有效合成至關重要，其功能喪失性突變會影響mtDNA 的正常維持，并導致線粒體疾?。?5］。ARNTL 是控制晝夜節律的時鐘基因之一，其高度甲基化可以促進鼻咽癌的發生并抑制其對順鉑的敏感性［16］。RAB38 屬于小GTP 結合蛋白家族，通過調控線粒體代謝以及線粒體介導的細胞死亡促進膠質母細胞瘤惡性進展［17］。IGFBP5 是胰島素樣生長因子結合蛋白家族的一員，參與線粒體功能的維持，沉默IGFBP5 的表達會導致神經母細胞瘤細胞線粒體凋亡的發生［18，19］。SERPINA5 屬于絲氨酸蛋白酶抑制劑超家族，可以作為低級別膠質瘤的預后因子，補充當前基于組織學的腫瘤分類系統［20］。FABP5 是脂肪酸結合蛋白家族的重要成員，在脂肪酸代謝中起著重要作用，一方面，FABP5 的缺失將損害心臟細胞的線粒體功能［21］；另一方面，FABP5 參與氧化應激時α 突觸核蛋白誘導的線粒體損害［22］。相比之下，關于GPR156 和COL22A1的研究相對較少，與腫瘤的關系還不是很清楚。模型的大部分基因都涉及線粒體功能的調控，一致的是，銅死亡的發生與線粒體密切相關，這也反映了本研究結果的可靠性。

值得注意的是，本研究存在一些局限性。首先，本研究在訓練集構建模型，通過測試集進行驗證，但這一模型并沒有在臨床隊列中得到進一步的驗證；其次，銅死亡與膠質瘤的關系需要進一步通過基礎實驗驗證；最后，本諾莫圖僅包含了風險評分和3 個臨床參數（性別、年齡和膠質瘤WHO分級），如果能納入更多的參數，將進一步提高預測性能，例如IDH1 突變情況，1p19q 缺失情況和MGMT 甲基化情況，這些都是與膠質瘤預后和治療反應密切相關的參數，但是由于公共數據庫并沒有詳細提供這些數據，導致沒辦法進行相關分析，下一步我們將收集自己醫院的膠質瘤樣本，詳細收集這方面的數據，再進行相關的分析。