內質網應激相關基因能預測骨肉瘤預后并與腫瘤免疫微環境相關

張 丹,周逸馳

(1.武漢大學中南醫院,武漢 430000;2.武漢市第四醫院 脊柱二骨腫瘤科,武漢 430000)

骨肉瘤起源于間葉組織,是一種侵襲性惡性腫瘤,高發于兒童與青少年,具有病情進展迅速、易轉移、易復發、預后不佳等特點[1-2],約有30%患者在確診5年后出現進展性轉移[3]。盡管立體定向放射治療、碳離子放射治療、標準化學治療、手術切除等綜合性治療方案近年來有了較大發展,但骨肉瘤患者的5年生存率仍不理想[4]。內質網( Endoplasmic reticulum, ER )參與蛋白質合成、加工與轉運,鈣離子的儲存以及脂質和類固醇激素的合成與轉運等多種生理活動。當細胞受到各種理化因素刺激時,內質網功能的內穩態體系被打破,由正常狀態轉變為應激狀態,稱為內質網應激(Endoplasmic reticulum stress, ERS ) 。ERS主要表現為ER腔內錯誤折疊與未折疊蛋白大量沉積、鈣離子平衡失調等;細胞通過激活未折疊蛋白反應(Unfold protein response, UPR)等信號通路引導和參與ER中異常蛋白質的再折疊和降解,防止細胞功能進一步受損。然而持續或過強的ERS則會誘導細胞凋亡信號通路的啟動。多種研究表明內質網應激與腫瘤的發生及進展有關,本研究基于內質網應激相關基因以生物信息學方法構建了骨肉瘤患者的風險模型,并探索了其與腫瘤免疫微環境的關系。

1 資料與方法

1.1 數據下載

訓練集的轉錄組數據及相應的臨床信息下載于UCSC Xena數據庫;驗證集的轉錄組數據與臨床信息下載于GEO數據庫(GSE21257,GSE39058),兩組數據使用R軟件包“inSilicoMerging”[5]進行了合并,隨后使用R軟件包“limma”的removeBatchEffect函數對合并后的數據去除批次效應。在GeneCards網站(www.genecards.org)搜索條目“endoplasmic reticulum stress”,并提取relevance score≥7.0的基因作為內質網應激相關基因。

1.2 風險模型的構建與驗證

使用R軟件包“survival”評估了每個基因的預后顯著性,提取logrankP<0.05的基因;隨后使用維恩圖(Jvenn)[6]對預后相關基因與ERS相關基因取交集。使用R軟件包“glmnet”的LASSO-cox方法提取了風險特征基因,此外還設置了5折交叉驗證以獲得最優模型。對上述獲得的基因使用R軟件包“survival”進行多因素COX回歸分析,最終篩選P<0.05的風險特征基因構建風險模型。風險評分=β1×expression of gene 1 +β2×expression of gene 2 +β3×expression of gene 3 + … +βn×expression of gene n 。使用五種臨床特征(種族、性別、年齡、轉移狀態及原發腫瘤部位)對風險評分進行分組并進行比較以驗證風險模型的獨立性。

1.3 風險模型與臨床特征的相關性

使用R軟件包“rms”建立了列線圖,并評估了生存時間、生存狀態、風險評分、臨床特征在樣本中的預后顯著性。隨后使用了R軟件包“survival”和“stdca.R”進行3、5年的決策曲線分析(Decision curve analysis, DCA)以評估該風險模型的臨床效用;最后評估了風險評分與5種臨床特征(種族、性別、年齡、轉移狀態及腫瘤原發部位)之間的關系。

1.4 生存分析

使用R軟件包“survival”與“survminer”以Kaplan-Meier生存曲線評估了高、低風險組之間預后差異的顯著性;隨后用R軟件包“pROC”的roc函數進行了1、3、5年的受試者工作特征曲線(Receiver operating characteristic curve, ROC)分析,并使用ci函數評估曲線下面積(Area under curve, AUC)和置信區間;隨后使用R軟件包“pheatmap”繪制了風險預后熱圖,顯示不同風險組樣本的生存狀態分布及風險特征在不同分組之間的表達情況。

1.5 免疫分析

使用ESTIMATE算法[7]評估了高低風險組之間免疫評分、基質評分和ESTIMATE評分之間的差別;使用ssGSEA方法[8]比較了28種不同的免疫細胞在高低風險分組樣本之間的浸潤豐度。隨后使用了MCP-counter算法[9]評估了8種不同的免疫細胞及兩種基質細胞(血管內皮細胞與纖維母細胞)在高低風險分組之間的浸潤豐度;并使用相關性散點圖評估了風險評分與免疫評分、基質評分及ESTIMATE評分的相關性。

1.6 功能富集分析

使用R軟件包“limma”[10]對高、低風險組之間進行了基因差異表達分析,隨后使用Metascape[11]對差異表達基因(Differentially expressed genes, DEGs)進行了進一步的分析。使用了R軟件包“clusterProfiler”[12]對DEGs進行了基因本體論(Gene ontology, GO)與京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes, KEGG)聯合功能富集分析,“GOplot”包[13]計算每個條目的zscore并對結果進行可視化。基于分子特征數據庫 (http：//www.gsea-msigdb.org/gsea/downloads.jsp) 下載了GO生物過程基因集以評估相關途徑及分子機制,使用R軟件包“GSVA”[14]計算了每個樣本在每個基因集中的富集得分。基于MSigDB數據庫,使用R軟件包“clusterProfiler”包進行了基因集富集分析(Gene set enrich analysis, GSEA)。

1.7 統計學分析

2 結果

2.1 構建ERS相關風險模型

從UCSC Xena數據庫(https：//xenabrowser.net/datapages/)下載了88例患者的臨床信息與相應基因表達數據作為訓練集,其中3例無隨訪信息予以剔除;GEO數據庫(https：//www.ncbi.nlm.nih.gov/geo/,GSE21257,GSE39058)共下載95例患者信息與相應基因表達數據作為驗證集,訓練集與驗證集臨床數據詳見表1。訓練集數據使用單因素COX回歸分析顯示共有7 714個基因與骨肉瘤預后相關,將該預后相關基因集與內質網應激相關基因集(n=770)取交集,提取出79個內質網應激-預后相關基因(見圖1a),這79個基因在人類染色體的坐標、表達情況及相互聯系如圖1b所示。隨后使用LASSO回歸分析(LASSO regression algorithm)進一步提取出26個候選基因,最佳λ值選擇為0.041 4(見圖1c-1d)。基于這26個候選基因,采用多因素COX回歸分析最終確定了6個風險特征基因：TMED10、VEGFA、MAPK10、TOR1B、PTGIS、SERPINH1,其中VEGFA、PTGIS及SERPINH1與骨肉瘤患者的不良預后相關,而TMED10、MAPK10及TOR1B與骨肉瘤患者的良好預后相關(見圖1e)。風險模型構建公式：風險評分= 0.239×VEGFA表達量-2.45362×TMED10表達量-0.94029×MAPK10表達量-1.45051×TOR1B表達量+ 0.79041×PTGIS表達量+ 1.03296×SERPINH1表達量。ROC分析顯示(見圖1f),該風險模型1年ROC曲線下面積(Area under curve, AUC)為0.95(95%CI：1.00-0.88),3年AUC為0.97(95%CI：1.00-.094),5年AUC為0.98(95%CI： 1.00-0.95)。訓練集決策曲線分析(Decision curve analysis, DCA)顯示,風險評分的預測的準確性高于其余臨床預測指標(年齡、性別、種族、原發腫瘤部位及轉移狀態,見圖1g)。以風險評分的平均值將訓練集患者分為高風險組(n=43)與低風險組(n=42),K-M生存曲線顯示,低風險組預后顯著優于高風險組(Log-rankP<0.001,見圖1h)。

表1 訓練集與驗證集的臨床數據特征Table 1 Clinical data characteristics of the training and validation sets

圖1 訓練集數據建立ERS基因相關風險模型Fig.1 Establishment of the ERS genes related risk model in the training cohort

2.2 基因差異表達分析與富集分析

訓練集中,對高、低風險組進行基因表達差異分析(設定差異倍數為1.5,顯著閾值為0.05),共有341個差異表達基因(Differential expressed genes, DEGs),其中有138個上調基因,203個下調基因(見圖2a)。Metascape分析顯示,DGEs共分為6個模塊,其主要富集于淋巴細胞介導的免疫力、細胞外基質組成、骨化、凋亡程序的調控及應對缺氧(見圖2b,2c)。GO-KEGG聯合富集分析結果顯示,DEGs主要富集于循環免疫球蛋白介導的體液免疫反應、補體激活的經典途徑、免疫球蛋白復合體、金色葡萄球菌感染(見圖2d)。GSVA結果顯示,肥大細胞活化參與的免疫反應、成熟B細胞分化參與的免疫反應、淋巴細胞活化參與的免疫反應、T細胞活化參與的免疫反應、B細胞活化參與的免疫反應、巨噬細胞活化參與的免疫反應、免疫反應激活富集于低風險組,而缺氧誘導因子1 α信號通路、對缺氧反應的內在凋亡信號通路、細胞對缺氧反應的調節、凋亡染色體凝縮、翻譯后蛋白質靶向內質網膜、內質網未折疊蛋白質反應富集于高風險組(見圖2e)。GSEA分析顯示,DEGs主要富集于獲得性免疫反應、先天性免疫反應、淋巴細胞與非淋巴細胞間的免疫調節作用及補體系統(見圖2f)。上述分析表明,與免疫細胞激活相關的信號通路主要富集于低風險組,而低風險組預后優于高風險組,表明低風險組的免疫景觀及免疫狀態可能優于高風險組;同時,與內質網過度應激的相關通路及缺氧反應的相關通路富集于高風險組,提示內質網應激及缺氧反應可能為導致骨肉瘤患者預后不佳的原因。

2.3 風險模型的獨立性

本研究基于五種臨床特征驗證了該風險模型的獨立性,種族(P=0.55)、性別(P=0.4)、年齡(P=0.5)、轉移狀態(P=0.16)及原發腫瘤部位(P=0.99、P=0.88、P=0.74)之間的風險評分均無統計學差異(見圖3a)。隨后,以性別(見圖3b,女性P<0.001,男性P<0.001)、轉移狀態(見圖3c,轉移P=0.002,非轉移P<0.001)、年齡(見圖3d,18歲以下P<0.001,18歲以上P=0.013)及種族(見圖3e,白種人P<0.001,非白種人P<0.001)分組進行生存分析(見圖3b),該風險模型仍舊顯示出了強大的預后預測能力,低風險組預后均顯著優于高風險組。

圖2 高低風險組DEGs相關分析Fig.2 DEGs related analysis of the high risk and low risk group

圖3 風險模型獨立性評估Fig.3 Independence assessment of the risk model

2.4 列線圖的構建

整合風險評分與臨床特征(性別、種族、原發腫瘤部位、轉移狀態)構建了列線圖以更加準確地預測骨肉瘤患者的預后(見圖4a),根據風險評分與上述臨床特征對骨肉瘤預后的影響,給予相應的評分,隨后分別在訓練集與驗證集中以校準曲線驗證了該列線圖結果,校準曲線顯示了該列線圖良好的預后能力(見圖4b,4c)。模型總體C指數為0.967 3(95%CI：0.944 8-0.989 8),P<0.001。

2.5 風險特征基因的特點

圖5a表明,該風險模型成功將訓練集的骨肉瘤患者分為高風險組與低風險組,在高風險組中VEGFA、PTGIS及SERPINHI高表達,而在低風險組中TMED10、MAPK10及TOR1B高表達。圖5b表明,VEGFA、PTGIS及SERPINHI在高風險組中的表達量高于低風險組(P<0.05),而TMED10、MAPK10及TOR1B在低風險組中的表達量高于高風險組(P<0.05)。生存分析(見圖5c)顯示,TMED10、MAPK10及TOR1B的高表達與良好預后相關(Log-rankP<0.05),而VEGFA、PTGIS及SERPINHI的高表達與不良預后相關(Log-rankP<0.05)。

2.6 免疫分析

ESTIMATE分析顯示,訓練集的低風險組的基質評分、免疫評分及ESTIMATE評分均高于高風險組(P<0.05),進一步表明低風險組與免疫反應相關(見圖6a)。相關性散點圖分析表明,風險評分與基質評分(P<0.001,R=-0.39),免疫評分(P<0.05,R=-0.26)及ESTIMATE評分(P<0.001,R=-0.365)呈負相關(見圖6b)。MCP-counter分析(見圖6c)顯示,T細胞、CD8 T細胞、B細胞系、NK細胞、單核細胞系及內皮細胞在低風險組中的浸潤豐度高于高風險組(P<0.05)。ssGSEA分析(見圖6d,6e)顯示,28種不同的免疫細胞主要富集于低風險組,其中活化B細胞、活化CD8 T細胞、記憶B細胞、1型輔助T細胞、2型輔助T細胞、髓源性抑制細胞、NK細胞及NK T細胞的富集評分有統計學差異(P<0.05)。上述結果表明低風險組具有較高的免疫景觀與免疫評分,而高風險組則具有較低的免疫景觀與免疫評分。

圖4 構建列線圖與校準曲線Fig.4 Establishment of nomograph and calibration curve

圖5 風險特征基因的特點Fig.5 Characteristics of risk genes

2.7 驗證風險模型

在驗證集中,風險評分同樣將骨肉瘤患者分為高風險組與低風險組,熱圖可見VEGFA、PTGIS及SERPINH1在高風險組高表達,而TMED10、MAPK10及TOR1B在低風險組高表達(見圖7a)。圖7b表明,VEGFA、PTGIS及SERPINH1在高風險組表達高于低風險組(P<0.05),而TMED10、MAPK10及TOR1B在低風險組表達高于高風險組(P<0.05)。圖7c表明低風險組預后優于高風險組。ROC分析顯示(見圖7d),1年AUC為0.59(95%CI： 0.72-0.45),3年AUC為0.63(95%CI： 0.76-0.49),5年AUC為0.62(95%CI： 0.77-0.48),顯示該預后模型在驗證集中也有較好的預測能力。上述分析結果與訓練集一致,進一步表明該風險模型具有優良的預測能力。

圖6 免疫分析Fig.6 Immune analysis

圖7 驗證風險模型的魯棒性Fig.7 Verification of the robustness of the risk model

3 討 論

骨肉瘤是兒童和青少年最常見的惡性骨腫瘤[15],盡管近年來綜合治療方案快速發展,但骨肉瘤的5年生存率仍不理想,目前迫切需要制定有效的風險分層方法[16-17]。內質網應激(Endoplasmic reticulum stress,ERS)在人類癌癥領域中的作用得到了廣泛關注,為骨肉瘤的個體化療法帶來了新的希望。內質網是蛋白質合成和細胞內鈣儲存最重要的細胞器,并參與多種細胞信號通路的調控[18-19]。在細胞內和細胞外應激誘導下,未折疊或錯誤折疊蛋白的積累可激活未折疊蛋白反應(Unfolded protein response, UPR),從而恢復內質網的內環境平衡[20]。PERK/ATF4/CHOP、IRE-1/XBP1和ATF6信號通路在調控ERS誘導的細胞生理反應中起重要作用[21-22],然而過度、持續的內質網應激激活也可誘導細胞凋亡[23]。同時,ERS參與了多種疾病的進程,如帕金森病、阿爾茲海默癥、糖尿病及腫瘤[24-26],對ERS的深入研究可為骨肉瘤的臨床治療提供重要參考。在惡性腫瘤中,癌細胞不斷暴露于缺氧的腫瘤微環境和細胞內DNA損傷應激中,ERS激活的風險較高,PERK-ATF4分支可激活上皮-間充質轉化(Epithelial-To-mesenchymal transition, EMT)反應,促進乳腺癌細胞轉移[27]。在肝細胞癌和前列腺癌中,來自UPR的轉錄因子可以結合并激活VEGF的啟動子,導致內皮細胞的增殖和遷移[28-29]。此外,ERS激活也能誘導癌細胞休眠, 細胞休眠通過將癌細胞阻滯在細胞周期的G0/G1期,從而阻斷癌細胞的分裂和增殖[30],處于靜止狀態生存的癌細胞對放化療不敏感,最終導致腫瘤復發,患者預后不良 。此外,相關研究證明IRE-1-XBP1信號通路可以抑制抗腫瘤免疫應答,為腫瘤的形成提供機會[31]。上述結論提示ERS可能在腫瘤免疫微環境(Tumor immune microenvironment, TIME)的重構中起重要作用,是免疫治療的潛在靶點。雖然多項研究證實了ERS在骨肉瘤中的治療價值,但ERS在骨肉瘤中的生物學功能尚未完全了解,特別是其在調節TIME的作用機制。

本研究基于UCSC Xena數據庫及GeneCards數據庫以生物信息學方法構建了一個ERS基因相關骨肉瘤風險模型,該風險模型包括6個風險特征基因：P24跨膜轉運蛋白10(Transmembrane P24 trafficking protein 10, TMED10)、血管內皮生成因子A(Vascular endothelial growth factor A, VEGFA)、絲裂原激活的蛋白激酶(Mitogen-activated protein kinase, MAPK10)、Torsin家族1成員B(Torsin family 1 member B, TOR1B)、前列腺素I2合酶(Prostaglandin I2 synthase, PTGIS)及絲氨酸蛋白酶抑制劑1(Serine proteinase inhibitor, SERPINH1),其中TMED10、MAPK10及TOR1B與骨肉瘤患者的良好預后相關,而VEGFA、PTGIS及SERPINH1與骨肉瘤患者的不良預后相關。MAPK10是肝細胞肝癌(Hepatocellular carcinoma, HCC)的重要抑癌基因,其表達丟失提示HCC患者的預后較差[32],另外MAPK10低表達也與宮頸癌患者的不良預后相關[33];VEGFA在多種腫瘤組織中高表達,不僅與腫瘤血管生成相關,還能直接或間接參加腫瘤免疫反應[34];PTGIS高表達可促進腫瘤相關巨噬細胞(Tumor-associated macrophages, TAMs)和調節性T細胞(T-regulatory cells, Tregs)在腫瘤微環境中的浸潤,并與肺癌、卵巢癌及胃癌患者的不良預后有關[35],另外PTGIS過表達也能作為結直腸癌肝轉移的預測指標[36];SERPINH1是多種腫瘤不良預后的危險因素和預測因子,同時其過表達與腫瘤免疫抑制狀態有關,SERPINH1可作為泛癌免疫治療的潛在靶點[37-38]。本研究的結果表明ERS過度激活與惡性臨床病理特征和基因組改變密切相關,同時這些結果揭示了ERS在骨肉瘤中的意義及其影響ERS激活的分子機制。該風險模型將骨肉瘤患者分為了高風險組與低風險組,生存分析表明低風險組預后優于高風險組。同時ROC分析、DCA均表明該風險模型能準確地預測骨肉瘤患者預后,隨后進一步驗證了該風險模型的獨立性。

本研究對高、低風險組患者的基因進行差異分析,并對DEGs進行富集分析,結果揭示了應對缺氧、ERS等相關通路富集于高風險組,而免疫反應相關通路富集于低風險組。TIME在腫瘤患者預后方面具有重要作用,因為腫瘤的發展與周圍基質的修飾有關,而免疫細胞是腫瘤基質的關鍵組成部分[39]。ESTIMATE算法可以根據基因表達值推斷腫瘤純度,以及腫瘤中免疫細胞和基質細胞的比例[7],因此能夠準確地反應TIME。本研究結果表明,具有高免疫評分、基質評分及ESTIMATE評分的骨肉瘤患者具有良好預后。另外,本研究還采用了另外兩種方法—MCP counter算法及ssGSEA來評估高低風險組腫瘤組織中的免疫細胞浸潤豐度。MCP counter算法結果表明高風險組中有5/8的不同免疫細胞豐度顯著降低,與ESTIMATE結果一致,表明免疫景觀(Immune landscape)在高風險組下調。ssGSEA分析概述了28種免疫相關細胞的豐度,結果表明高風險組患者處于相對較低的免疫狀態,進一步證實了ESTIMATE和MCP counter的結果。在腫瘤微環境中,ERS在腫瘤的免疫機制中也發揮了重要作用,其不僅決定了腫瘤細胞的命運,而且還間接觸發了免疫抑制反應,例如ERS的中性粒細胞通過凝集素樣氧化低密度脂蛋白受體-1(Lectin-like oxidized low-density lipoprotein receptor-1,LOX-1)的表達獲得免疫抑制活性[40],口腔鱗狀細胞癌的ERS可傳遞給中性粒細胞,導致中性粒細胞上調其表面受體LOX-1并導致其免疫功能受到抑制[41],綜上所述,缺氧、過度ERS等腫瘤免疫微環境的變化可能引起了高風險組患者的免疫抑制狀態,而免疫反應的激活可能與低風險組患者的良好預后相關。

雖然已有學者從腫瘤微環境、免疫細胞浸潤和能量代謝等方面構建骨肉瘤風險模型[42-44],但本研究仍具有獨特的優點：(1)本研究聚焦于ERS,構建了由6個風險特征基因組成的風險模型,并將骨肉瘤患者分為了具有顯著生存差異的高、低風險組;(2)本研究探索了不同分組之間的DEGs,并進一步探索了DEGs的富集通路,結果表明這些通路與免疫反應有關;(3)本研究驗證了該風險模型的準確性與獨立性,并進一步在GEO隊列中驗證了該風險模型;第四,本研究探索了這六種風險特征基因的特點,其結果在訓練集與驗證集中一致;最后,免疫分析表明高風險組具有較低的免疫景觀及免疫評分,推測過度ERS可能導致免疫抑制狀態,從而引起高風險組的不良預后。

4 結 論

以生物信息學方法挖掘公共數據庫,篩選出6個內質網應激相關風險特征基因,并以此為基礎構建了骨肉瘤患者的風險模型,為骨肉瘤患者的預后提供了新思路、新方法,具有一定的新穎性;本研究所得到的基因可作為骨肉瘤基礎研究和治療的潛在靶點,有一定的應用價值;不足之處在于本研究仍具有一定的局限性：一方面,本研究結果來源于生物信息學分析,并沒有得到進一步的實驗驗證;另一方面,本研究使用的數據從開放數據庫下載獲得的,且回顧性研究的低證據水平的性質仍然存在,還需要更多的前瞻性研究來進一步證實ERS在骨肉瘤中的預后價值;最后,關于該風險模型與骨肉瘤患者免疫狀態的關系需要進一步的實驗驗證。