方 萌，池晴佳，趙 晗，郭 敬，劉宇煒*

1江漢大學醫(yī)學院基礎醫(yī)學部，湖北 武漢 430056；2武漢理工大學理學院力學系，湖北 武漢 430070

肝細胞癌（hepatocellular carcinoma，HCC）是原發(fā)性肝癌的主要病理類型，在世界范圍內(nèi)，原發(fā)性肝癌是癌癥相關死亡的第四大原因，其預后通常較差［1］。在現(xiàn)有的診療體系下，目前大多數(shù)HCC 患者被診斷時已為晚期階段。新的診斷和預后生物學標志物的開發(fā)和運用，對于提高HCC患者的生存率具有重要意義。

SOX 轉(zhuǎn)錄因子家族是人體非常重要的一類轉(zhuǎn)錄因子，具有高度保守的高遷移率（high mobility group protein，HMG）序列［2-4］，在胚胎發(fā)育和細胞分化中發(fā)揮重要作用，近年研究顯示它還與腫瘤的發(fā)生、侵襲轉(zhuǎn)移有關［3，5-6］。據(jù)報道，SOX轉(zhuǎn)錄因子家族已在各種腫瘤組織中作為腫瘤抑制因子或啟動子發(fā)揮作用［7］。在臨床HCC 患者中，SOX 基因家族成員異常表達，其簽名基因與腫瘤等級和腫瘤階段密切相關。新建立的SOX簽名可以在訓練、測試和獨立驗證隊列中強有力地預測患者的總體生存［8］。SOX1 降低與HCC 的不良預后和腫瘤進展相關［9］。Guo 等［10］發(fā)現(xiàn)HCC 組織中SOX9 蛋白的過度表達對腫瘤的進展和不良預后具有預測價值。SOX12 作為SOX 基因家族的一員，其功能同樣復雜，影響各種類型惡性腫瘤疾病的發(fā)生發(fā)展［11-12］。有報道SOX12 在HCC 中過表達，并通過上調(diào)Twist1 促進腫瘤細胞的轉(zhuǎn)移［13］。但尚不清楚HCC 中SOX12 表達的臨床病理和預后意義。

HCC的轉(zhuǎn)移與免疫微環(huán)境的改變密切相關，深入研究SOX12 與免疫浸潤的關系可能幫助改善HCC 預后的預測。然而，目前的研究對SOX12 在HCC 患者的預后價值及其與免疫細胞浸潤的相關性尚不清楚。加權基因共表達網(wǎng)絡（weighted gene co?expression network，WGCNA）和最小絕對收縮和選擇算法（least absolute shrinkage and selection oper?ator，LASSO）在快速、準確識別這些基因標志物，構建風險預后模型中有著顯著的效果［14-15］。本課題組前期研究了肝癌及其他疾病的預后標志物與免疫調(diào)節(jié)分子機制［16］。本研究擬借助WGCNA與LASSO算法建立HCC風險預后模型，并探究SOX12表達與HCC的免疫浸潤和基因突變間的相關性，為判斷其預后價值提供借鑒。

1 資料和方法

1.1 資料

腫瘤基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫（https：//tcga?data.nci.nih.gov/tcga/）下載患者基因表達譜的3 級數(shù)據(jù)及臨床信息，研究基因表達與生存之間的關系?；虮磉_（Gene Expression Omnibus，GEO）數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo/）查找肝細胞癌芯片，獲取數(shù)據(jù)集（GSE76427和GSE14520）用于確認和驗證TCGA 數(shù)據(jù)集結(jié)果。GEPIA 2（http：//gepia2.cancer?pku.cn/）和ENCORI（http：//starbase.sysu.edu.cn）在線數(shù)據(jù)庫進行單基因的篩選。

1.2 方法

1.2.1 評估免疫細胞浸潤水平

使用ESTIMATE 算法確定所有樣本的免疫評分，基 于CIBERSORT（http：//cibersort.stanford.edu/）的反卷積算法計算訓練集和驗證集中22 種免疫細胞在每個患者中的浸潤比例。使用Wilcoxon 檢驗評估上述數(shù)據(jù)集中高低SOX12表達的免疫細胞浸潤差異，P<0.05為具有顯著的浸潤差異。對上述細胞進行后續(xù)分析，評估其浸潤水平對患者預后的影響。

1.2.2 獲取免疫相關基因

利用TCGA數(shù)據(jù)庫中HCC的相關數(shù)據(jù)獲得全基因組。從ImmPort 數(shù)據(jù)庫（https：//immport.niaid.nih.gov）中檢索免疫相關基因（immune?related gene，IRG）數(shù)據(jù)。通過韋恩圖分析全基因組和IRG 之間的重疊免疫相關基因，再用單因素Cox 回歸篩選出P<0.05的基因進行進一步分析。

1.2.3 WGCNA的構建

通過測量一定功率和網(wǎng)絡連接強度的無標度拓撲擬合指數(shù)，保證相應的軟閾值功率和無標度網(wǎng)絡。對1～30 的閾值功率進行網(wǎng)絡拓撲分析，確定WGCNA 的規(guī)模獨立性和平均連通性，并采用層次聚類功能將表達譜相似的基因聚類成模塊。為了進一步確定共表達網(wǎng)絡中的功能模塊，進行了拓撲重疊測度（TOM）：

TOMij表示基因i 和j 的共享相鄰基因的重疊。基于TOM 的不相似度法，進一步平均連鎖層次聚類，具有相似表達的被歸類到相同的基因模塊中。

1.2.4 基于LASSO算法的風險預后模型

單變量Cox回歸分析用于探索基因?qū)傮w生存的影響。如果基因P<0.05，則整合到LASSO 算法中，生成嚴格等于0 的回歸系數(shù)，為了數(shù)據(jù)降維，去除一些權值較低的變量，防止協(xié)變量共線性導致的過擬合。計算每個樣品的風險評分（RS）：

其中N 為基因的個數(shù)，Expi為基因表達譜，Ci為LASSO分析確定的回歸系數(shù)。以Kaplan?Meier曲線評估風險評分與總生存率之間的關系。計算受試者工作特征（receiver operating characteristic，ROC）曲線的曲線下面積（area under curve，AUC）。

1.2.5 突變景觀分析

在R軟件中使用maftools軟件包，分析高低風險組和SOX12高低表達對基因突變頻率的影響。

1.3 統(tǒng)計學方法

利用R 軟件（https：//www.r?project.org/）進行統(tǒng)計學分析。使用“survival”軟件包進行生存分析，評估免疫評分與總生存期之間的關系。取中值作為分界點。使用“pROC”軟件包生成與時間相關ROC（the time?dependent ROC，tROC）曲線。使用“forest?plot”包進行單變量和多變量Cox 回歸分析。使用“glmnet”軟件包進行LASSO分析。在整個研究過程中，兩組獨立樣本間采用t檢驗，P<0.05 為差異具有統(tǒng)計學意義。

2 結(jié)果

2.1 SOX12表達在HCC中的臨床價值

利用TCGA?HCC 免疫相關基因和突變基因的交集產(chǎn)生了297 個基因。GEPIA 2（http：//gepia.can?cer?pku.cn/）和ENCORI（http：//starbase.sysu.edu.cn/index.php）數(shù)據(jù)庫分析了這些基因的生存和差異表達，篩選出生存和差異同時滿足P<0.05 的8 個基因。其中CCR7、VAV1、OX40、PSMD4 已經(jīng)在文獻中進行了大量報道［17-21］，OXTR、PMCH 僅用于動物實驗研究［22-23］，而有關PNOC的信息量較少。因此，最終確定了滿足條件的單基因SOX12（圖1）。

圖1 單基因SOX12的確定Figure 1 Determination of single gene SOX12

通過TCGA數(shù)據(jù)庫基因數(shù)據(jù)（表1）分析，與相鄰正常組織相比，HCC組織中SOX12的表達顯著上調(diào)（圖2A）。根據(jù)Kaplan?Meier 分析，高SOX12表達的HCC 患者其總生存期較差（P=0.024，圖2B）。由tROC 可知，SOX12的表達對患者預后具有較強的預測能力，其1年、3年和5年的AUC分別為0.679、0.636和0.617（圖2C）。根據(jù)疾病分期和腫瘤分級情況顯示，SOX12 的表達在不同分期中具有明顯差異（圖2D、E、F）。為了確保結(jié)果的可靠性，使用GEO數(shù)據(jù)庫（表2）中的GSE76427（n=116）進行了驗證（圖3）。

圖3 GEO?LIHC隊列中SOX12臨床價值分析Figure 3 Analysis of the clinical value of SOX12 in GEO?LIHC cohort

表2 GEO數(shù)據(jù)庫臨床信息表Table 2 The clinical data in the GEO database（n）

圖2 TCGA?LIHC隊列中SOX12的臨床價值分析Figure 2 Analysis of the clinical value of SOX12 in the TCGA?LIHC cohort

表1 TCGA數(shù)據(jù)庫臨床信息表Table 1 The clinical data in the TCGA database（n）

通過對TCGA和GEO數(shù)據(jù)進行單因素和多因素Cox回歸分析，結(jié)果顯示，SOX12表達是一個獨立的預后生物標志物（圖4A），這一顯著性差異在GEO數(shù)據(jù)中也得已驗證（圖4B）。以上結(jié)果表明，SOX12表達水平可能影響HCC的發(fā)展，可作為HCC預后不良的獨立生物標志物。

2.2 SOX12高低表達的腫瘤免疫微環(huán)境

利用CIBERSORT算法觀察22種腫瘤免疫細胞在HCC組織中的聚類情況（圖5A、6A），在SOX12高表達和低表達組中上述腫瘤免疫細胞比例和亞群分布有顯著性差異（圖5B～C、6B～C）。此外，在TCGA數(shù)據(jù)中，靜息記憶性CD4+T 細胞、濾泡輔助性T 細胞、調(diào)節(jié)性T細胞、巨噬細胞M0、巨噬細胞M1、激活的肥大細胞在SOX12 高低表達中具有顯著性差異（圖5D）。在GEO 數(shù)據(jù)中，初始B 細胞、靜息記憶性CD4+T 細胞、濾泡輔助性T 細胞、巨噬細胞M2 在SOX12 高低表達中具有明顯差異（圖6D）。通過對TCGA 和GEO 數(shù)據(jù)集的免疫細胞交集，獲取交集的兩個免疫細胞：靜息記憶性CD4+T 細胞和濾泡輔助性T細胞進行后續(xù)分析。

圖5 TCGA?LIHC隊列分析HCC中的免疫浸潤情況Figure 5 Immune infiltration of HCC in the TCGA?LIHC cohort

圖6 GEO?LIHC隊列分析HCC中的免疫浸潤情況Figure 6 Immune infiltration of HCC in the GEO?LIHC cohort

2.3 WGCNA確定與免疫標記物有關的模塊

通過TCGA和IRG交集獲取的1 333個基因，進行單因素分析篩選出P<0.05的基因共1 090個（圖7A）。利用WGCNA 分析SOX12 表達與免疫標志物的相關基因，選取靜息記憶性CD4+T 細胞的標志物CCR5、CD69、FOXP3 及濾泡輔助性T 細胞的標志物CXCR5、IL21、CD44作為臨床性狀分析，經(jīng)過聚類篩選，共鑒定出7個不同的免疫相關基因表達模塊（圖7C）。在黃色、藍色、青綠色、棕色、綠色、紅色、灰色模塊中分別包含66、283、404、106、41、20 和170 個免疫相關基因。將HCC 患者的模塊特征基因與免疫標志物進行相關性分析（圖7F），發(fā)現(xiàn)SOX12表達與青綠色模塊顯著相關，選取青綠色模塊內(nèi)基因用于預后模型分析。

圖7 WGCNA分析Figure 7 WGCNA analysis

2.4 預后模型分析

青綠色模塊的404個基因確定為與SOX12生存相關的IRG。經(jīng)過LASSO 回歸分析選定出18 個顯著影響患者生存預后的IRG（ACKR1、CALCR、CCR3、CD8A、HMOX1、IL15RA、IL12A、IL18RAP、KLRK1、NOD2、PGF、PLXNC1、RNASE2、S100A9、SLC11A1、SSTR2、SPP1、TNFRSF11B）（圖8A、B）。根據(jù)風險評分系數(shù)，將HCC 患者分為高低風險組（圖8C）。由Kaplan?Meier 分析顯示出高風險組比低風險組的總生存率較差（圖8D）。tROC曲線分析表明，風險預后模型具有顯著的預后預測效果，其1年、3 年和5 年的AUC 分別為0.823、0.811 和0.824（圖8E）。

圖8 LASSO分析Figure 8 LASSO analysis

2.5 高低表達和風險組之間不同的突變頻率

利用TCGA數(shù)據(jù)確定SOX12基因表達和IRG的高低風險是否與基因的突變頻率相關。結(jié)果顯示，高SOX12 表達組中TP53 基因突變比例（40%）明顯高于低SOX12 表達組（25%）（圖9A、B）。這一結(jié)果和IRG 高低風險組結(jié)果類似（圖9C、D）。上述數(shù)據(jù)表明高低SOX12 表達和基因簽名的高低風險分組影響基因突變的頻率，說明SOX12可能是HCC的標志基因。

圖9 HCC中基因突變的瀑布圖Figure 9 Waterfall chart of genetic mutations in HCC

3 討論

轉(zhuǎn)錄因子SOX 基因家族是一類新發(fā)現(xiàn)的編碼轉(zhuǎn)錄因子的超基因家族，先前研究發(fā)現(xiàn)SOX12在胚胎發(fā)育和細胞特征維持中起著至關重要的作用［24］。越來越多的證據(jù)表明，SOX12 基因的突變、缺失或過表達與多種類型惡性腫瘤的形成和發(fā)展密切相關［25-26］。研究人員發(fā)現(xiàn)SOX12表達有助于維持HCC的腫瘤細胞特性［27］，并通過激活HCC 中上皮?間質(zhì)轉(zhuǎn)化過程而引起腫瘤的轉(zhuǎn)移［28］。但SOX12 作為HCC中一個重要的生物標志物，其與HCC患者的生存預后關系尚不清楚。相比先前的研究，本文通過差異表達分析、ROC 曲線分析和生存分析，首次揭示了SOX12 在HCC 中的診斷和預后價值。與相鄰的非腫瘤組織相比，HCC組織中的SOX12表達顯著上調(diào)，并且其高表達與患者的不良生存密切相關，可能成為影響HCC預后的關鍵標志物。

通過TCGA 數(shù)據(jù)庫基因數(shù)據(jù)分析，與相鄰的正常組織相比，HCC 組織中SOX12 的表達顯著上調(diào)。根據(jù)Kaplan?Meier 分析，高SOX12 表達的HCC 患者其總生存期較差（P=0.024）。由tROC 可知，SOX12的表達對患者的預后具有較強的預測能力，其1、3、5 年的AUC 分別為0.679、0.636 和0.617。根據(jù)疾病分期和腫瘤分級情況顯示，SOX12 的表達在不同分期中具有顯著差異。為了確保結(jié)果的可靠性，使用GEO 數(shù)據(jù)GSE76427 進行了驗證。對TCGA 和GEO 數(shù)據(jù)進行單因素和多因素Cox 回歸分析顯示，SOX12 表達與患者生存預后顯著相關。這一結(jié)果表明，SOX12 表達可能影響HCC 的發(fā)展，可作為HCC 預后不良的獨立生物標志物。之前研究結(jié)果也證實SOX12 與多種癌癥的不良預后顯著相關［29］，其在腫瘤相關的免疫浸潤中起重要作用［30］，并且與TNM 分期顯著相關［13］，本文結(jié)果與這些研究報道相一致。基于WGCNA與LASSO 回歸分析選定出18 個顯著影響患者生存預后的IRG，并獲得了相應的風險評分。

盡管臨床上使用TNM 分期和分子特征對不同的腫瘤亞型、階段和臨床預后進行分類，但新出現(xiàn)的證據(jù)表明免疫環(huán)境可能對于腫瘤類型分類和分期以及預測患者生存率非常有用，甚至更為優(yōu)越［5］。先前已有報道，T細胞存在于HCC的免疫細胞浸潤中，T細胞的免疫浸潤程度與HCC患者的存活率密切有關［31-32］，本研究觀察到與文獻報道相同的現(xiàn)象。此外本研究還發(fā)現(xiàn)SOX12 高低表達與CD4+T 細胞和濾泡輔助性T 細胞的免疫浸潤水平具有顯著差異，提示SOX12 在免疫微環(huán)境中的作用是通過這兩類免疫細胞浸潤體現(xiàn)的。隨后，我們研究了SOX12基因高低表達組和高低風險組的基因突變頻率，發(fā)現(xiàn)TP53、CTNNB1 等基因突變在高SOX12基因表達組中更常見。據(jù)報道稱TP53突變體是在癌細胞表達中最常見的蛋白質(zhì)突變體［33］。值得一提的是，CTNNB1在散發(fā)性腫瘤不顯示突變，而在可遺傳的患者所占比例較高［34］，提示某些基因突變與SOX12表達和腫瘤進展相關。這與本研究結(jié)果一致，后期將進一步通過實驗驗證以逐步完善SOX12 表達對HCC 影響的證據(jù)。

根據(jù)風險評分系數(shù)將HCC 患者分為高低風險組，Kaplan?Meier 分析顯示出高風險組比低風險組的總生存率較差。tROC 曲線分析表明基于SOX12表達的風險預后模型1、3、5 年的AUC 分別為0.823、0.811、0.824，比之前報道的肝癌預后模型表現(xiàn)更好［35-36］。雖然前期研究人員分別構建了4 基因、6 基因肝癌預后模型［37-38］，但這兩項研究沒有定量評估預后模型的生存預測能力。而本預后模型具有很好的預測能力，有望幫助加強肝癌的預后預測能力。

總之，利用WGCNA 與LASSO 回歸得到的基于SOX12 表達的風險預后模型具有很好的預后預測能力。此外，SOX12高低表達組中靜息記憶性CD4+T 細胞和濾泡輔助性T細胞的免疫浸潤水平表現(xiàn)出顯著差異，SOX12的高表達不僅會影響HCC患病基因的突變頻率，而且與HCC 的不良預后顯著相關。本研究提示SOX12 可能成為新型的HCC 預后生物標志物，為HCC后續(xù)免疫治療和靶向治療的進一步開展提供參考。