加權基因共表達網絡與機器學習算法確定宮頸癌中免疫浸潤相關的基因簽名與預后分析①

方萌 趙晗池 晴佳（江漢大學醫學院，武漢430056）

宮頸癌（cervical cancer，CC）是世界上第四大最常見的女性癌癥。近年來，其發病率呈現年輕化趨勢［1］。盡管已經證實CC是由HR-HPV（HPV16和18）持續感染引起的［2］。然而，CC的高復發率和異質性導致其預后較差。因此，開發新的診斷和預后生物標記非常迫切。MicroRNA（miRNA）的異常表達被認為是各種癌癥發展的關鍵因素［3］。作為miR?NA家族的成員，miR-141-3p影響多種癌癥進展，并且在CC的發生和發展過程中發揮著重要作用［4-5］。已有研究報道miR-141-3p通過靶向單個基因影響CC細胞的生長和侵襲［6］。但在CC的預后預測方面，基于預后模型的多基因簽名比單基因的準確度更高。考慮到免疫微環境對腫瘤發生發展的關鍵作用，深入簽名基因的免疫浸潤有可能幫助改善CC的預后預測。先前研究尚未探索miR-141-3p的靶基因與免疫微環境之間的相關性。因此，全面分析miR-141-3p的靶基因與總體生存率之間的相關性，建立可靠的免疫相關的多基因簽名以準確預測CC病人的生存至關重要。機器學習在快速、準確的確定基因標記物中起到了重要作用［7-9］。并且，諾謨圖比普通的風險評分模型擁有更好的預后預測能力。在以往的研究中，已結合實驗和生物信息學方法來研究肝癌與其他疾病的預后標記物與免疫調節分子機制［10-12］。本文通過TCGA數據庫分析miR-141-3p在CC中的預后性能。使用WGCNA、COX回歸、LASSO COX回歸建立并驗證了與免疫浸潤相關的多基因簽名與諾謨圖，以預測患者的生存。

1 材料與方法

1.1 數據預處理CC（TCGA數據庫中簡寫為CESC）的mRNA，miRNA表達數據和臨床數據可從TCGA數據庫（https：//tcga-data.nci.nih.gov/tcga/）下載。表達數據進行log2轉換和歸一化處理，排除無隨訪消息或<1 d的時間臨床數據。

1.2 miR-141-3p靶基因的篩選使用miRWalk2.0（http：//zmf.umm.uni-heidelberg.de/apps/zmf/mir?walk2/）中的12個數據庫（Microt4，miRWalk，mirbridge，miRanda，miRDB，miRMap，Pictar2，PITA，MiRNAMap，RNAhybrid，RNA22和Targetscan）來確定miR-141-3p的靶基因。然后，使用R語言中的“limma”程序包分析CESC數據。以FDR<0.01，|log2FC|≥2為閾值篩選差異表達基因（DEGs）。最后，通過Venn Plot分析DEGs與預測基因之間的重疊基因（候選基因）。

1.3 加權基因共表達網絡分析通過R軟件中的“WGCNA”軟件包分析候選基因的表達和臨床文件。通過計算基因對之間的Pearson相關性，模塊特征基因與性狀（包括疾病狀態）的相關性，確定與預后相關的基因模塊。

1.4 基于LASSO COX回歸的風險預測模型通過單變量Cox回歸分析探討了每個模塊基因對總體生存的影響。將P<0.05和與生存相關的模塊基因整合到機器學習算法（LASSO回歸）中，以鑒定預后風險特征。使用R包“timeROC”繪制tROC曲線，預測1年、3年和5年總體生存。

1.5 基于支持向量機免疫細胞浸潤基于支持向量機的CIBERSORT算法用于計算每個患者中22個免疫細胞的比例。通過wilcoxon檢驗評估了高危和低危患者之間免疫細胞浸潤的差異（P<0.05）。TIMER用于分析簽名基因與免疫細胞浸潤之間的相關性。

1.6 基因集富集分析通過“clusterprofiler”軟件包（Perm=1 000，minGSSize=100和P.value=0.01）分析了高危和低危患者之間的信號通路多樣性。“limma”程序包用于計算基因的logFC值。

2 結果

2.1 miR-141-3p的臨床價值miR-141-3p在腫瘤組織中的表達顯著上調（P=0.003 2，圖1A）。ROC曲線表明miR-141-3p的表達可以很好地將腫瘤與正常樣品區分開（AUC=0.871 19，圖1B）。以303個腫瘤組織中miR-141-3p的平均值作為臨界值，將所有患者分為低表達和高表達患者。排除生存時間缺或生存時間為0的情況得到的患者樣本為280個。生存分析結果表明，miR-141-3p的低表達可顯著改善患者的預后，而高表達對應不良的預后（P=0.022 0，圖1C）。tROC曲線分析表明，miR-141-3p對患者的預后具有很強的預測能力。1年、3年和5年的AUC分別為0.762、0.715和0.741（圖1D）。

圖1 miR-141-3p的臨床價值Fig.1 Clinical value of miR-141-3p

2.2 鑒定預后相關基因模塊miR-141-3p靶基因包含1 485個交集基因（圖2A）。選擇軟閾值=5以符合無標度網絡法則（圖2B）。基于基因和樣本數據的層次聚類建立的網絡熱圖（圖2C）。在9個模塊和多個臨床表型參數的相關分析結果中，發現黃色模塊同時與T.stage（cor=0.28，P=1e-06），N.stage（cor=-0.23，P=5e-05）和AJCC階段（cor=0.22，P=2e-04）具有很強的相關性（圖2D）。因此，確定黃色模塊為與CC預后相關的功能性miR-141-3p下游基因集。

圖2 WGCNA分析Fig.2 WGCNA analysis

2.3 風險預測模型與諾謨圖隨機將280例CC患者的樣本分為訓練集（n=140）和驗證集（n=140）。通過單變量COX回歸確定了與訓練集中的患者OS顯著相關的31個基因（P<0.05）。利用LASSO回歸進一步優化結果，最后有5個靶標顯著影響患者的生存率（表1，圖3A）。基于LASSO COX分析的結果，建立了5個靶標簽名風險預后模型。模型公式為：RS=0.012 1×FOXA1+0.007 7×DMBX1+0.149 3×TMEM98+0.014 2×RHPN1+0.026 4×SRMS基 于 訓練集的生存分析表明，與高風險相比，低風險可以顯著改善患者的預后（P=0.000 54）。tROC曲線分析表明，該風險預后模型具有顯著的預后預測效果，其1年、3年和5年AUC分別為0.810、0.821和0.679（圖3B）。此外，基于驗證集K-M曲線分析（P=0.000 86）和tROC分析也證實了上述結果（圖3C）。

表1 與OS相關的5個預后靶標Tab.1 5 prognostic targets related to OS

圖3 風險預測模型構建與驗證Fig.3 Risk prediction model construction and verification

通過單變量和多變量COX回歸評估了患者的風險評分和多種臨床特征。分析顯示，在訓練和驗證組中，風險評分和AJCC分期與患者的OS顯著相關（表2）。通過組合風險評分和AJCC分期而建立的諾謨圖（C-index=0.83，圖4A）在預測患者生存率方面明顯優于單獨使用風險評分（C-index=0.77）和AJCC分期（C-index=0.70）。基于訓練集的研究發現。1年、3年和5年OS預測表明，由諾謨圖預測的生存率與最佳預測性能非常匹配（圖4B）。與以上結果一致，諾謨圖（C-index=0.68，圖4C）與預后因素風險評分（C-index=0.67）和AJCC分期（C-index=0.56）相比，具有更好的預后和預測能力（驗證集）。同時，校準曲線還驗證了驗證集中諾謨圖的預后預測性能（圖4D）。

表2 預后信息的單因素和多因素分析Tab.2 Univariate and multivariate analyses of prognostic informations

圖4 建立預后模型的諾謨圖Fig.4 Establishment of nomogram of prognostic model

2.4 低危和高危人群免疫細胞浸潤和GSEA通路差異基于CIBERSORT的280個腫瘤樣品中22個免疫細胞的浸潤率（圖5A、B）。CD4+T細胞，巨噬細胞M0和M2是腫瘤微環境中3個最豐富的免疫細胞。有趣的是，巨噬細胞和T細胞的浸潤率最大，約占免疫細胞的65%。在TCGA訓練集中，高危和低危患者之間存在5種類型（P<0.05）的免疫細胞浸潤差異（圖5C）。根據TCGA驗證集，在不同風險組中存在6種類型（P<0.05）的免疫細胞浸潤差異（圖5D）。并且通過TIMER獲得了預后標志介導的免疫浸潤變化（圖6）。結果表明，樹突狀細胞的浸潤與5個基因在預后標志物中的表達顯著相關。基于全基因組GSEA分析，在高危和低危患者中獲得了6條重要的KEGG通路（圖7）。

圖5 基于CIBERSORT算法分析22種免疫細胞浸潤Fig.5 Analysis of 22 immune cells infiltration according to CIBERSORT algorithm

圖6 基于TIMER的5種預后基因表達與腫瘤浸潤性免疫細胞的相關性分析Fig.6 Correlation analysis between 5 prognostic gene expression and tumor infiltrating immune cells based on TIMER

圖7 GSEA獲得了高危組和低危組間存在差異的6條通路Fig.7 GSEA obtained 6 pathways with differences be?tween high-risk group and low-risk group

3 討論

首先通過TCGA數據集挖掘了miR-141-3p的臨床價值。結果表明，miR-141-3p在CC組織中的表達明顯上調，與LI等［6］的研究一致。miR-141-3p高表達水平表明CC患者生存期較差。tROC曲線表明，miR-141-3p的表達對患者的預后具有強大的預測能力。然后，對miR-141-3p的潛在靶標進行了全面分析。WGCNA指出目標基因模塊與T.Stage，N.stage和AJCC階段顯著相關。

越來越多的研究表明，基因標記在預測腫瘤的預后中起著重要的作用［13-15］。迄今為止，還未有研究通過將WGCNA和LASSO Cox回歸方法相結合以鑒定miR-141-3p的預后靶基因標志物來提高CC預后的預測。在本次研究中，通過對miR-141-3p的候選基因進行了WGCNA和LASSO Cox回歸確定了5個與預后相關的靶基因。在以前的報道中，FOXA1、DMBX1、TMEM98、RHPN1和SRMS與CC細胞遷移發展顯著相關［16-20］。

作為一種機器學習算法，LASSO COX已有大量研究［21-23］。對于生存數據，最常用的統計模型是Cox比例風險回歸模型［24］。而其他傳統機器學習算法如隨機森林、人工神經網絡、貝葉斯算法并沒有專門開始針對生存數據的模塊。本研究的重點是預測CC患者的總體生存。作為機器學習算法的LAS?SO COX算法能夠識別潛在的風險因素，評估預測、擬合優度，以及結果的臨床相關性［25］。故本文采用LASSO COX回歸作為機器學習算法。LASSO COX模型的優點在于可以直接解釋患病風險和生存之間的關系。另一方面，機器學習技術是無假設和數據自適應，這意味著它們可以有效地用于對復雜數據進行建模。當然，LASSO COX算法也存在一定的局限。例如，該算法要求輸入數據進行適當的處理，并且計算過程中會對參數進行適當的調整，以避免算法性能下降［26］。

研究表明，多基因簽名在預測腫瘤的預后中起著至關重要的作用。LIANG等［27］提出了一種3-miRNA簽名，用于CC的獨立預后預測。XIE等［28］還確定了1個8基因簽名來預測CC患者的預后。此外，LI等［29］使用組蛋白家族基因簽名來預測CC患者的生存。同樣，DING等［30］提出了1個3基因簽名，以通過Cox回歸改善CC的生存預測。先前的研究使用單變量Cox回歸分析方法來預測CC患者的生存。但當前對CC預后簽名的研究用到的算法不夠深入，對簽名基因的免疫浸潤角色也極少關注。本文通過結合WGCNA、LASSO Cox回歸以及支持向量機算法得到的5個基因簽名和諾謨圖，有效改善了CC預后的預測。WGCNA使用無監督聚類來識別基因模塊，默認方法是分層聚類。分層聚類是一種用于確定多維空間中相似數據點聚類的常用方法，它可以使用動態樹切割方法來確定模塊［31-32］。動態樹切割方法成功地識別了使用靜態切割方法無法識別的分支，產生預后相關的基因簇。為了結果的準確性，本實驗進一步通過單變量COX回歸分析基因簇以獲得與總體生存相關的基因，利用LASSO COX算法則解決了肝癌從高維數據獲得的預測模型可能具有過擬合的風險，構建了預后的特征。基于風險評分和臨床特征，通過多變量COX回歸建立了諾謨圖。而研究結果也證實，諾謨圖在預測患者生存率方面明顯優于單獨使用風險評分（C-index=0.77）和AJCC分期（C-index=0.70）。因而，基于預后簽名的風險評分顯示出生存預測具備潛在優勢。

確定5個基因的特征后，構建5個基因的預后模型并研究其預后價值。發現訓練組和驗證組的高風險和低風險水平的預后顯著不同。所有高危人群的OS均低于低危人群。tROC分析表明5個基因簽名具有很強的預測能力。腫瘤的進展不僅受腫瘤本身特征的影響，而且還受腫瘤微環境的影響。值得注意的是本實驗的預后標志與腫瘤微環境的免疫浸潤水平顯著相關。根據CIBERSORT分析，發現在不同的風險組中，B細胞，漿細胞和巨噬細胞M1、M2具有較大的豐度差異。此外，GSEA還揭示了高風險組與低風險組相比有6種不同的免疫狀態。因此，研究結果可為CC患者更好的治療提供依據。

總之，綜合運用加權基因共表達網絡、COX回歸、LASSO回歸算法，開發并驗證了基于miR-141-3p的預后靶基因的新型5個基因簽名（FOXA1、DMBX1、TMEM98、RHPN1和SRMS），并進行了預后分析。結果顯示，與免疫浸潤相關的5個基因簽名可以顯著改善CC患者的預后。因此，本次研究提出的5個基因簽名與諾謨圖有潛力為CC患者提供預后評估工具。