基于生物信息學對食管鱗狀細胞癌新輔助放化療敏感性生物標志物的發掘*

杜芝霖，伍健，徐鵬，郎錦義

646000 四川 瀘州，西南醫科大學附屬醫院 腫瘤科（杜芝霖、伍健、郎錦義）；610041 成都，四川省腫瘤臨床醫學研究中心，四川省腫瘤醫院·研究所，四川省癌癥防治中心，電子科技大學附屬腫瘤醫院放療中心/放射腫瘤學四川省重點實驗室（徐鵬、郎錦義）

根據全球癌癥統計，2020 年有超過604 000 例食管癌新發病例和544 000 例死亡病例，全球發病率第七，死亡率第六［1］。中國是食管癌的高發地區，2015 年中國癌癥統計，食管癌的發病率排名第三，死亡率排名第四［2］。其中，食管鱗狀細胞癌（esophageal squamous cell carcinoma，ESCC）是食管癌中最常見的組織學類型。

ESCC 的早期癥狀并不明顯，許多患者直到癌癥晚期才被明確診斷。盡管近幾十年出現了新的治療策略和方案，但食管癌的預后仍然很差，5 年生存率約為20%［3-4］。而如今，新輔助放化療（neoadjuvant chemoradiotherapy，NCRT）聯合根治性食管切除術已逐漸成為食管癌的主要治療方法［5］。然而，不同患者接受NCRT 的結果不一樣。獲得病理完全緩解（pathologic complete response，PCR）的患者比未獲得PCR 的患者有更好的生存率［6-7］。但一些臨床參數例如TNM 分期、腫瘤位置和腫瘤大小等并不能預測食管鱗癌患者對NCRT 是否敏感［8］。如何有效地預測食管鱗癌患者對NCRT 治療的敏感性是本研究的目的。

隨著高通量測序技術和生物信息學的發展，對腫瘤發生機制的認識有了很大的提升［9］。有相關研究提出，一些基因水平的變化可以預測NCRT 敏感性，這有助于ESCC 患者的個體化治療［10-11］。如SOX17基因在食管鱗癌患者中的過表達對NCRT更為敏感［12］。研究人員也通過建立35 個基因突變譜的預測模型，來預測食管癌患者接受NCRT 治療的敏感性［13］。還有研究發現了能預測食管癌患者NCRT 治療敏感性的5 個基因［14］。因此基因的異常表達或突變會影響到患者接受NCRT 的敏感性。在本研究中，通過使用生物信息學方法，篩選出與食管鱗癌NCRT 治療敏感性相關的具有預測價值的生物標志物，并采用免疫組化方法進行初步驗證，以期獲得新的治療靶點，為食管癌患者的個體化治療帶來希望。

1 材料和方法

1.1 數據源和數據處理

從GEO 數據庫中篩選與食管鱗癌NCRT 敏感性相關的數據集，納入標準：①提供食管癌患者RNA-Seq 的數據集；②包含食管癌NCRT 后不敏感患者的組織與敏感患者的組織。最終篩選出數據集GSE45670［15］。

回顧性收集了四川省腫瘤醫院2017～2021 年食管鱗癌患者的相關信息。納入標準：① 患者在疾病確診前未進行過任何抗腫瘤治療；② 患者年齡18～75 歲，KPS 評分 ＞ 80；③ 我院初治的中晚期胸段食管鱗癌并且無遠處轉移的患者 （cT1b-cT2N+或cT3-cT4任何N），并于我院行NCRT 及手術治療；④ 所有患者治療前均需完善相關基線檢查，明確病灶情況。所有患者新輔助同步放化療治療后第28天復查增強CT 進行初次篩選NCRT 敏感性。影像學提示腫瘤病灶沒有消退，稱之為疾病穩定（stable disease，SD）或是腫瘤病灶較前增大，稱之為疾病進展（progression disease，PD）的患者歸類于放化療不敏感；影像學提示腫瘤病灶完全消失，稱之為完全緩解（complete response，CR）或者腫瘤病灶較治療前緩解 ＞ 50%的患者歸類于放化療敏感。放化療后約6 周，患者接受食管鱗癌根治性手術切除原發腫瘤以及區域淋巴結。術后病理是否獲得PCR 作為最終劃分食管鱗癌患者接受NCRT 敏感性的標準。

1.2 差異表達基因（differentially expressed genes，DEGs）的獲得及蛋白互作（protein-protein interaction，PPI）分析

使用R 軟件的“limma”包對GSE45670 數據集的NCRT 敏感組和NCRT 不敏感組的 DEGs 進行篩選，閾值設置為P＜ 0.05 和LogFC ＞ 1，并使用“ggplot2”包進行了熱圖和火山圖的繪制［16-17］。

將篩選出的DEGs 導入String 數據庫（https://string-db.org）進行PPI 網絡分析。通過設置PPI 最低互作得分（0.4），去除該網絡中不相關的節點，得到初次篩選的PPI 網絡［18］。將其導入Cytoscape 軟件中將PPI 網絡可視化。通過Cytoscape 軟件的MCODE 插件進一步過濾，設置MCODE 分數 ＞ 4，節點數 ＞ 5［19］。使用該軟件中Cytohubba 和CytoNCA 兩個插件，運用Degree 算法，篩選出核心差異基因［20-21］。

1.3 構建加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）

利用WGCNA 分析方法對數據集GSE45670構建共表達網絡［22］。首先，篩選出方差前25%的5 837 個基因。其次，計算每對基因之間的皮爾遜相關性，得到一個相似性矩陣。第三，使用“pick Soft Threshold”函數來選擇理想的軟閾值（β = 5）。第四，利用冪函數將相似度矩陣轉換為鄰接矩陣。然后，將鄰接矩陣轉換為拓撲重疊矩陣（topoligical overlap matrix，TOM）。最后，建立了基于TOM 的不相似度（1-TOM）的層次聚類。每個模塊樹狀圖至少包含50 個基因，將相似表達模塊聚類為相同模塊的閾值設為0.25，最終得到模塊。

1.4 篩選候選基因進行基因富集分析

使用R 軟件的“VennDiagram”包，將DEGs 和WGCNA 中相關度最高的模塊內的基因進行交叉重疊后，結合PPI 網絡中獲得核心基因和候選基因。

對候選基因進行基因通路富集分析，使用R 包“clusterProfiler”和“enrichplot”包對其進行注釋，對可能參與的生物過程和信號通路進行相應的分析［23-24］。

1.5 對候選基因進行生存分析

從TCGA 數據庫中篩選食管鱗癌患者行Kaplan-Meier 生存分析，獲得食管鱗癌患者基因表達水平與生存時間之間的關系。在Log-rank 檢驗中，P＜ 0.05 為差異有統計學意義。

1.6 基質金屬蛋白酶（matrix metalloproteinase，MMP）9 和醛固酮還原酶（aldo-keto reductase，AKR）1C3 的免疫組化分析

對已明確診斷為ESCC 患者的活檢樣本進行免疫組化分析。采用抗體MMP9（1∶100）和抗體AKR1C3（1∶200）進行蘇木精染色。目標蛋白的表達量通過使用Image-Pro Plus 圖像軟件進行分析。MMP9 陽性蛋白主要在細胞質和細胞膜中表達為褐色。AKR1C3 陽性蛋白表現出中度到較強的細胞質反應，顏色為褐色。通過測量出每張圖片的累積光密度值 （integrated option density， IOD）以及區域面積值（area），再計算出平均光密度值（mean density），即 mean density = IOD/area，此值反映了目標蛋白的單位面積濃度。最后取每個樣本的 3 個隨機區域mean density 的平均值即為此樣本的值。

1.7 統計學方法

利用R 軟件（4.2.1）、Cytoscape 軟件（3.8.0）、ImageProPlus 軟件、String 在線網站（https: //string-db.org/）進行數據分析。采用Kaplan-Meier 法進行生存分析，并采用Log-rank 檢驗。利用TCGA 數據集，來獲取相關候選基因的表達水平的高低與食管鱗癌患者總生存期之間的關系。并且對其進行了相應曲線的繪制。免疫組化的結果采用GraphPad Prism 7.0 進行統計分析。采用Mann-WhitneyU檢驗或雙尾非配對檢驗來確定顯著性水平。本研究中使用的所有數據均滿足統計檢驗的假設。所有定量數據均以均值±標準差呈現。P＜ 0.05 為差異有統計學意義。

2 結 果

2.1 數據收集

從GEO 數據庫中篩選得到的數據集GSE45670，包含了28 個食管鱗癌患者的樣本，其中有17 個食管鱗癌NCRT 治療不敏感的患者和11 個食管鱗癌NCRT 治療敏感的患者。而臨床樣本共收集了26 例活檢組織標本，包括11 例NCRT 治療敏感者，也是影像學上提示腫瘤病灶CR 的患者，和15 例NCRT治療不敏感者，其中11 例患者影像學提示腫瘤病灶SD，4 例患者影像學提示腫瘤病灶PD。

2.2 DEGs 的篩選及PPI 網絡分析

以LogFC ＞ 1，P＜ 0.05 為閾值，共篩選得到350個差異基因，其中上調基因166 個，下調基因184 個（圖1）。

圖1 差異基因篩選結果Figure 1. Volcano Plot of DEGs

運用String 在線數據庫構建了350 個DEGs 的PPI 網絡，再導入Cytoscape 軟件中，運用MCODE 插件進行篩選，再使用Degree 算法，運用該軟件中的CytoHubba 和CytoNCA 插件最終確定MMP9是該網絡的核心基因（圖2）。

圖2 由差異基因構成的PPI 網絡Figure 2. Protein-Protein Interaction Networks Composed of Differentially Expressed Genes

2.3 加權基因共表達網絡的構建

WGCNA 是近年來廣泛使用的一種技術，旨在避免僅僅分析DGEs 而造成忽視核心基因的可能性［25-26］。通過選取相關系數為0.85 時對應的軟閾值為 5（圖3A）構建無尺度網絡進行 WGCNA，再通過動態雜化剪切合并相似的模塊，獲得13 個模塊（圖3B）。其中沒有被分配到任何集群中的灰色模塊，在本研究中沒有進行進一步的分析。結果顯示，紫色模塊與食管鱗癌患者接受NCRT 敏感性的相關性最高，共159個基因，差異具有統計學意義（P= 0.03；圖3C）。

圖3 構建WGCNA 模塊Figure 3. Construction of WGCNA Modules

2.4 獲得候選基因

將紫色模塊中的159 個基因與350 個DEGs 進行交集，提取到5 個重疊基因（圖4）。將這5 個重疊基因與PPI 網絡中的核心基因結合，得到6 個候選基因，包括AKR1C3、GPX2、RRBP1、GSTM3、AKR1C1和MMP9。

圖4 DEGs 和紫色模塊基因的韋恩圖Figure 4．Screening of Candidate Genes: Venn Diagram of Intersections between DEGs and the Purple Module

2.5 候選基因的功能富集分析

這些候選基因共得到335 個相關的GO 條目。選出其中10 個GO 條目（圖5A）。同時，分別從生物過程、細胞組分和分子功能這三個方面選取了相關的條目（圖5B）。基因和基因組的KEGG 信號通路分析結果發現，這些基因主要富集于“氧化應激反應”中，也可能與“細胞解毒”“細胞間橋”等途徑有關（圖5C）。

圖5 候選基因的功能富集分析Figure 5. Functional Enrichment Analysis of the Candidate Genes

2.6 候選基因的生存分析

從TCGA 的食管癌數據庫中篩選出81 個食管鱗癌患者的臨床數據來對這6 個候選基因進行生存分析，通過生成最佳截止點，形成基因高表達和低表達兩組。從圖6 可知，AKR1C3和AKR1C1基因高表達組患者與低表達組患者相比，總生存率（overall survival，OS）明顯較差，差異具有統計學意義（P= 0.011，P= 0.035），而食管鱗癌患者的OS 并不受其它4 個基因的表達影響，差異不具有統計學意義（均P＞ 0.05）。

圖6 TCGA 數據庫中ESCC 中候選基因的Kaplan-Meier 生存曲線Figure 6. Kaplan-Meier Survival Curves of Candidate Genes in ESCC from TCGA Database

2.7 通過免疫組化檢測食管鱗癌標本中AKR1C3和MMP9 的表達

對GSE45670 數據集中6 個基因的表達進一步分析，P＜ 0.01 被認為是基因表達的顯著性差異，發現AKR1C3、AKR1C1、MMP9和RRBP1在NCRT 不敏感的組織和NCRT 敏感組織中的表達差異更為顯著（圖7）。通過對本研究收集的活檢標本進行免疫組化分析，與NCRT 敏感的組織相比，NCRT 不敏感組織中AKR1C3 蛋白表達上調（圖8）；MMP9 蛋白則在NCRT 敏感組織中表達上調（圖9）。

圖7 GSE45670 數據集中6 個候選基因的接受NCRT 食管鱗癌不敏感組織和接受NCRT 食管鱗癌敏感組織之間的基因表達差異Figure 7．Expressions of 6 Candidate Genes in Non-responder and Responder to NCRT in GSE45670 Dataset

圖8 AKR1C3 在NCRT 敏感組織和 NCRT 不敏感組織表達水平的比較（DAB 法，×20）Figure 8．Expression Levels of AKR1C3 in Responder and Non-responder to NCRT (DAB, ×20)***P ＜ 0.001.

圖9 MMP9 在 NCRT 敏感組織和 NCRT 不敏感組織表達水平的比較（DAB 法，×20）Figure 9．Expression Levels of MMP9 in Responder and Non-responder to NCRT（DAB，×20）

3 討 論

本研究利用生物信息學對基因集GSE45670 進行分析。首先通過DEGs 分析篩選出 350 個差異基因。其次，基于WGCNA 生物信息學方法，得出與食管鱗癌放化療治療中異質性顯著相關的159 個基因，與350 個DEGs 進行交叉，獲得了5 個重疊基因（AKR1C1、AKR1C3、GPX2、GSTM3、RRBP1）。 同時，基于350 個DEGs，生成了相應的PPIs 網絡。利用MCODE、CytoHubba 和CytoNCA 這三個插件，將其導入Cytoscape 軟件進行可視化，得到了1 個核心基因（MMP9）。最后，得到了6 個候選基因，包括AKR1C1、AKR1C3、GPX2、GSTM3、RRBP1和MMP9。通過進一步對這6 個候選基因進行功能富集分析發現，“氧化應激反應”途徑在其中起著非常重要的作用。放射治療通過多種方式誘發腫瘤細胞的死亡，例如產生活性氧（reactive oxygen species，ROS）、DNA損傷和亞細胞細胞器中的應激反應。ROS 是由電離輻射所誘發的，不僅與蛋白質、DNA 損傷、細胞膜表面酶和基因組不穩定相關，還引發損傷修復，可以破壞氧化還原穩態，從而導致氧化應激反應，最終導致細胞死亡［27-29］。并且化療會誘導產生大量的ROS，從而導致癌細胞死亡。該通路對于癌癥的發展、治療的耐藥性和基因組的不穩定性等方面都有著不可忽視的作用。

同時，通過設置P＜ 0.01 為閾值，得出AKR1C3/1、MMP9和RRBP1在GSE45670 數據集中NCRT 敏感組和NCRT 不敏感組的表達差異更為顯著。

RRBP1 是一種內質網膜蛋白，在核糖體結合、新生蛋白易位以及內質網應激和未折疊蛋白反應中起重要作用。對于新生蛋白核糖體的結合和轉運，RRBP1 是必不可少［30］。許多研究證實，RRBP1 在許多惡性腫瘤中過表達，例如宮頸癌、肺癌、膀胱癌和結直腸癌，并且影響患者的預后［31-34］。同時，也有研究證實，RRBP1 可能通過調節YAP1 的表達而成為順鉑耐藥的主要驅動因素［35］。這些證據都說明，RRBP1基因水平的變化與癌癥的發生發展具有一定的相關性。

MMP9 是一種參與多種生物活性的基質蛋白酶，是金屬蛋白酶類中分子量最大的一種。MMP9在作用于細胞外基質（extracellular matrix，ECM）的構成和膜蛋白裂解中扮演著重要的角色，它可以切割多種ECM 蛋白來使得ECM 重構。MMP9 對基底膜的降解過程也起著非常重要的作用，基底膜中含有膠原，包括IV 型膠原，而這些膠原可被MMP9 降解，而基底膜破裂通常是腫瘤形成過程中腫瘤發生侵襲和轉移的關鍵階段［36-37］。目前，有部分研究已經證實，MMP9 在一些惡性腫瘤中過表達，例如非小細胞肺癌、宮頸癌和胰腺癌［38-40］。同時，一些研究提出，MMP9 與放射治療有關，部分原因是ROS 的形成可以破壞鋅離子和Cys99 之間的相互作用，從而激活MMP9，增加MMP9 的表達，影響放射治療效果［41-42］。血漿MMP9 水平在肺癌、乳腺癌和肝癌放療期間會升高［43-44］。在晚期非小細胞肺癌患者中，血清MMP9 水平在放射治療后顯著降低，但在SD和PD 的患者中，放射治療后的血清MMP9 水平沒有變化［45］。ESCC 患者NCRT 后，健康食管組織近端甚至遠端MMP9 水平升高［46］。在最近的一項研究中，研究發現MMP9 的表達與小細胞肺癌患者的順鉑耐藥有關［47］。因此，MMP9 與癌癥病理廣泛相關，不僅在腫瘤侵襲、轉移和血管生成中發揮作用，同時還介導了腫瘤的微環境。

AKR1C3 是細胞內的一種多功能酶，是AKRs的一員。AKR1C1/3基因具有較高的序列同源性。AKR1C3 主要在內分泌器官中表達，參與腎上腺和腫瘤中類固醇的新生物合成，通過調節各種激素的合成以及其與相應受體結合的量，影響了腫瘤的發生發展［48］。AKR1C3 表達也會影響到與上 皮 間 充 質轉 化 （epithelial-mesenchymal transitions， EMT）相關的許多轉錄因子，而EMT 也是促進腫瘤轉移、腫瘤侵襲的主要作用機制之一［49］。有相關研究提出，基因AKR1C3 在乳腺癌和前列腺癌中表達上調，是其預后不良的標志［50-51］。在惡性腫瘤細胞中，AKR1C3基因表達上調可消除ROS，積累PGF2，進一步刺激MAPK 信號，防止G2/M 期阻滯或減少凋亡，促進癌細胞增殖，而導致其抗輻射［52-53］。AKR1C3 在非小細胞肺癌中表達較高，而在小細胞肺癌中表達則相反。而小細胞肺癌在臨床中通常對放射治療非常敏感［54］。同時，AKR1C3 表達上調可維持腫瘤細胞內的還原環境，不僅導致ROS 濃度降低，還可能與多西他賽耐藥和蒽環類藥物耐藥有關［55-56］。因此，AKR1C3 不僅與多種惡性腫瘤的發生和發展密切相關，并且與多種抗腫瘤治療也有一定的聯系。

同時，本研究通過在食管鱗癌組織標本進行免疫組化分析，發現MMP9 和AKR1C3 蛋白表達水平差異與ESCC 患者對NCRT 敏感性具有一定的相關性。AKR1C3 陽性蛋白在食管鱗癌接受NCRT治療不敏感的組織中有著中度到較強的細胞質反應，如圖8 所示。MMP9 陽性蛋白在食管鱗癌接受NCRT 敏感的組織中的細胞質和細胞膜中有著更高的表達，如圖9 所示。生存分析顯示，AKR1C1 和AKR1C3 高表達的食管癌患者的OS 低于AKR1C1和AKR1C3 低表達的食管癌患者。但是，我們的研究也有相應的缺陷。首先，與ESCC 患者NCRT 敏感性相關的樣本數量是有限的，需要在更大的隊列中評估這些基因的疾病預測效用。其次，本研究提出的生物標志物應進行進一步實驗驗證，以明確其確切的生物學行為。

綜上所述，本研究結果表明，基因MMP9和AKR1C3的異常表達與ESCC 患者接受NCRT 敏感性的生物學特性有著密切的關系，這為實現食管癌患者個體化治療提供了線索。

作者聲明：本文全部作者對于研究和撰寫的論文出現的不端行為承擔相應責任；并承諾論文中涉及的原始圖片、數據資料等已按照有關規定保存，可接受核查。

學術不端：本文在初審、返修及出版前均通過中國知網（CNKI）科技期刊學術不端文獻檢測系統的學術不端檢測。

同行評議：經同行專家雙盲外審，達到刊發要求。

利益沖突：所有作者均聲明不存在利益沖突。

文章版權：本文出版前已與全體作者簽署了論文授權書等協議。