齊延旭, 馬 川, 朱 勇, 黃晟棟, 趙華強, 張風(fēng)河
(山東大學(xué)齊魯醫(yī)學(xué)院口腔醫(yī)學(xué)院,口腔醫(yī)院口腔頜面外科,山東省口腔組織再生重點實驗室,山東省口腔生物材料與組織再生工程實驗室,山東 濟南 250012)
自噬是一個高度保守的真核細(xì)胞循環(huán)過程,在細(xì)胞生存和維持中發(fā)揮著重要作用[1-2]。已有的研究表明,自噬在腫瘤發(fā)生、發(fā)展中是一把雙刃劍,在腫瘤的形成和進展中起著相反的作用[3-6]。在口腔鱗狀細(xì)胞癌(OSCC)中,自噬通過Beclin-1/Atg7/Atg12-Atg5通路抑制OSCC進展[7];另一方面,也有研究顯示,自噬促進了OSCC的干性和耐藥性,線粒體自噬可能是導(dǎo)致OSCC化療耐藥性的原因[8]。目前已開發(fā)出許多促進OSCC自噬的藥物,并取得了良好的治療效果。然而,這些研究主要集中在自噬對OSCC進展和治療的影響,很少有研究者研究自噬在OSCC預(yù)后中的作用[9]。
本研究通過分析OSCC中差異表達的自噬相關(guān)基因(ARGs),建立OSCC患者的風(fēng)險評分模型,驗證其獨立預(yù)后價值,并構(gòu)建列線圖,以期對OSCC患者預(yù)后評估和臨床治療提供新思路。
我們從人類自噬數(shù)據(jù)庫 (human autophagy database,HADb)(http://www.autophagy.lu/)中 獲 得232個自噬相關(guān)基因,然后從癌癥基因組圖譜(the cancer genome atlas,TCGA)數(shù)據(jù)庫(https://portal.gdc.cancer.gov/)下載372個(340例OSCC患者和32例正常人)FPKM標(biāo)準(zhǔn)化的轉(zhuǎn)錄組測序(RNA-seq)數(shù)據(jù)和OSCC隊列的臨床病理和生存信息。從基因表達匯編(GEO)數(shù)據(jù)庫(http://www.ncbi.nlm.nih.gov/geo/)下載GSE41613數(shù)據(jù)集(196例OSCC)作為驗證集。
引用“l(fā)imma”包對TCGA下載數(shù)據(jù)進行數(shù)據(jù)處理[10],利用WilcoxTest統(tǒng)計分析OSCC樣本與正??谇唤M織樣本中差異表達的自噬相關(guān)基因。之后利用錯誤發(fā)現(xiàn)率(false discovery rate,F(xiàn)DR)方法對每個基因通過WilcoxTest統(tǒng)計分析得到的P值進行校正,得到校正后的P值(adj.P);并計算自噬基因在腫瘤組織中表達量的均值,將該均值除以正常組織中表達量的均值得到foldchange。設(shè)置篩選標(biāo)準(zhǔn):adj.P<0.05,|log2(foldchange)|>1,即篩選出腫瘤中表達量是正常組織中表達量2倍以上的自噬基因。
使用“org.HS.eg.db”包對自噬相關(guān)基因進行Symbol-ID轉(zhuǎn)換,使用“DOSE”“clusterProfiler”和“enrichplot”包進行基因本體論(gene ontology,GO)和京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes,KEGG)富集分析[11],挑選出每個基因的ID所對應(yīng)的GO功能和KEGG信號通路,將這些差異基因與上述挑選出的GO功能和KEGG信號通路在背景中所對應(yīng)的基因進行比較,得到P值,檢驗這些差異自噬基因在這些GO功能和KEGG信號通路上是否有影響。最后使用“barplot”包進行可視化,設(shè)置以P<0.05作為納入標(biāo)準(zhǔn),定義在差異表達基因中顯著富集的GO功能和KEGG信號通路。
在TCGA訓(xùn)練集中,利用log2(x+1)轉(zhuǎn)換將OSCC隊列中RNA-Seq數(shù)據(jù)進行標(biāo)準(zhǔn)化,將標(biāo)準(zhǔn)化后差異表達的ARGs表達量與生存信息合并,使用“survival”包進行單因素Cox比例風(fēng)險回歸分析,以確定與OS顯著相關(guān)的ARGs。然后使用“glmnet”包進行套索(LASSO)回歸分析[12],LASSO根據(jù)每個樣本的表達量,約束參數(shù)防止模型過擬合,選擇最佳“λ”值,從而篩選重要預(yù)后基因,為每個預(yù)后基因分配相應(yīng)的回歸系數(shù),構(gòu)建風(fēng)險評分模型。風(fēng)險評分計算公式如下:
風(fēng)險評分(risk score)=Gene1表達水平×β1+Gene2表達水平×β2+……+Gene n的表達水平×βn(β為回歸系數(shù),n為預(yù)后相關(guān)ARGs的總數(shù)目)。
將所有TCGA訓(xùn)練集中的OSCC患者按照風(fēng)險評分中位值分為高風(fēng)險組和低風(fēng)險組。利用“survival”包繪制高、低風(fēng)險組的風(fēng)險評分分布曲線、生存情況分布圖及預(yù)后相關(guān)基因表達量熱圖。繪制Kaplan-Meier(K-M)生存曲線評價2組患者的生存差異,采用對數(shù)秩(log-rank)檢驗評估其顯著性。為了評估預(yù)后模型的預(yù)測準(zhǔn)確性,用“survival ROC”包構(gòu)建了1、3、5年生存率的時間依賴受試者工作特征(receiver operating characteristic,ROC)曲線,并計算曲線下面積(area under the curve,AUC)。對GEO驗證集(GSE41613)采用同樣的方法檢驗風(fēng)險評分模型的普遍性和可靠性。
為確定風(fēng)險評分是否為OSCC的獨立預(yù)后因素,利用“survival”包進行單因素及多因素Cox回歸分析,比較TCGA訓(xùn)練集中風(fēng)險評分與其他臨床病理因素(年齡、性別、病理分級、腫瘤分期、T分期、N分期)的預(yù)后相關(guān)性,若單因素及多因素Cox回歸分析均顯示與預(yù)后顯著相關(guān)(P<0.05),則表示該因素為OSCC獨立預(yù)后因素。
為了預(yù)測OSCC患者1、3、5年生存的概率,我們利用“rms”包,結(jié)合年齡、性別、病理分級、腫瘤分期、T分期、N分期及風(fēng)險評分,在訓(xùn)練集中構(gòu)建列線圖,然后繪制校準(zhǔn)曲線,以圖形方式評估實際生存率與預(yù)測生存率之間的一致性[13]。
利用CIBERSORT算法繪制訓(xùn)練集中OSCC患者免疫細(xì)胞的比例[14];隨后通過“corrplot”包進行Pearson相關(guān)性分析,確定免疫細(xì)胞之間的共表達模式;最后通過非參數(shù)WilcoxTest進行檢驗比較高、低風(fēng)險組免疫細(xì)胞浸潤的差異,探究高、低風(fēng)險組患者預(yù)后差異的潛在免疫機制。
所有統(tǒng)計分析均采用R4.0.5軟件(https://www.r-project.org/)進行。采用Wilcoxon秩和檢驗評估風(fēng)險評分與臨床病理特征的相關(guān)性。在所有的統(tǒng)計檢驗中,以P<0.05為差異有統(tǒng)計學(xué)意義。
我們首先從HADb數(shù)據(jù)庫中獲得232個ARGs,然后從TCGA數(shù)據(jù)庫中下載OSCC隊列中372例患者的RNA-seq數(shù)據(jù)和臨床及預(yù)后數(shù)據(jù)。以FDR<0.05,|log2(foldchange)|>1為篩選標(biāo)準(zhǔn),得到37個差異表達的ARGs,其中11個為下調(diào)ARGs,26個為上調(diào)ARGs(圖1A、B)。腫瘤組織與正常組織間差異表達的ARGs的表達情況見圖1C。
對37個差異表達的ARGs進行功能富集分析。GO富集分析表明,在生物學(xué)過程中,ARGs主要富集于神經(jīng)元死亡、凋亡信號通路的調(diào)節(jié),細(xì)胞對未折疊蛋白的反應(yīng)等;在細(xì)胞成分中,ARGs主要富集于整合素復(fù)合物、參與細(xì)胞黏附的蛋白質(zhì)復(fù)合物和自噬體膜;在分子功能方面,ARGs主要富集在細(xì)胞因子活性、受體配體活性、信號受體激活劑活性等方面(圖1D)。在KEGG通路中,ARGs主要富集于細(xì)胞凋亡、鉑類藥物耐藥、EGFR酪氨酸激酶抑制劑耐藥、ErbB信號通路及癌癥中的PD-L1表達和PD-1檢查點通路等(圖1E)。
圖1 OSCC和正常對照樣本中ARGs的表達差異和基因功能富集分析Figure 1 The expression difference and gene functional enrichment analysis of ARGs in OSCC and normal samples
經(jīng)單因素Cox回歸分析,在TCGA數(shù)據(jù)集中共篩 選 出6個ARGs(BID、NKX2-3、DDIT3、VEGFA、FADD、BIRC5)作為預(yù)后相關(guān)基因(P<0.05),森林圖顯示了每個與預(yù)后相關(guān)的ARG的危險比(harzard ratio,HR)和相應(yīng)的置信區(qū)間(圖2A),這些預(yù)后基因被納入進一步的LASSO回歸分析(圖2B)以篩選最佳建?;?,圖2C為6個基因的系數(shù)分布圖,以此構(gòu)建風(fēng)險評分模型。在這些ARGs中,BID、DDIT3、VEGFA、FADD、BIRC5的HR>1,提示高表達與患者預(yù)后不良有關(guān),而NKX2-3的HR<1,其低表達預(yù)示著患者的不良預(yù)后。風(fēng)險評分計算如下:
圖2 差異表達的ARGs的單因素Cox回歸分析及LASSO回歸分析Figure 2 Univariate Cox regression and LASSO regression analysis of ARGs with differential expression
風(fēng) 險 評 分=(0.197×BID表 達 量)+(-0.409×NKX2-3表達量)+(0.091×DDIT3表達量)+(0.226×VEGFA表達量)+(0.052×FADD表達量)+(0.072 7×BIRC5表達量)。
我們在TCGA數(shù)據(jù)集內(nèi)部和GEO驗證集中驗證了ARGs風(fēng)險評分模型的性能。將TCGA數(shù)據(jù)集中的OSCC患者按照風(fēng)險評分中位數(shù)分為高風(fēng)險組和低風(fēng)險組。低風(fēng)險組NKX2-3的表達水平高于高風(fēng)險組,另一方面,高風(fēng)險組中其他自噬相關(guān)預(yù)后基因表達水平較高,將2組的風(fēng)險評分與生存數(shù)據(jù)進行可視化(圖3A)。時間依賴性ROC曲線顯示,ARGs預(yù)后模型對OS的預(yù)測效果較好,1、3、5年的AUC值分別為0.61、0.66、0.67(圖3B)。生存分析結(jié)果顯示,高風(fēng)險組的總體生存率明顯低于低風(fēng)險組(圖3C)。在GEO驗證集中用同樣的方法進行驗證,得到的結(jié)果與TCGA測試集相似。在GEO數(shù)據(jù)集中,風(fēng)險評分分布、生存狀態(tài)和基因表達模式如圖4A所示,除保護性基因NKX2-3外,其余5個ARGs在高風(fēng)險組的表達水平均高于低風(fēng)險組。GEO數(shù)據(jù)集中風(fēng)險評分1、3、5年的AUC值分別為0.65、0.67、0.68(圖4B)。高風(fēng)險組的總體生存率顯著低于低風(fēng)險組(圖4C)。
圖3 基于ARGs的風(fēng)險評分模型在TCGA數(shù)據(jù)集顯示出良好的預(yù)測能力Figure 3 ARGs-based risk score model displays favorable prediction ability in the TCGA datasets
圖4 基于ARGs的風(fēng)險評分模型在GEO驗證集顯示出良好的預(yù)測能力Figure 4 ARGs-based risk score model displays favorable prediction ability in the GEO verification datasets
風(fēng)險評分在單因素(HR=2.522,95%CI=1.496~4.249,P<0.001)及多因素(HR=3.131,95%CI=1.933~5.073,P<0.001)Cox回歸分析中均顯示,其與OSCC患者的預(yù)后顯著相關(guān)。另外,在其他臨床病理因素中,年齡在單因素(HR=1.028,95%CI=1.010~1.046,P<0.01)及多因素(HR=1.023,95%CI=1.006~1.039,P<0.01)Cox回歸分析中與患者OS顯著相關(guān)(圖5A、B)。以上結(jié)果說明風(fēng)險評分和年齡這兩個因素具有獨立預(yù)后價值,可作為OSCC的獨立預(yù)后因素。ROC曲線驗證了ARGs風(fēng)險評分模型的預(yù)測準(zhǔn)確性,風(fēng)險評分AUC值為0.627,高于其他臨床特征(圖5C)。
圖5 獨立預(yù)后因素評估Figure 5 Evaluation of independent prognostic factors
為了確定自噬相關(guān)風(fēng)險評分模型是否與OSCC進展有關(guān),我們分析了TCGA數(shù)據(jù)集中風(fēng)險評分與臨床病理變量之間的相關(guān)性。≤65歲與>65歲2組間(P>0.05,圖6A),男性與女性性別間(P>0.05,圖6B)以及病理分級G1~2與病理分級G3~4 2組間(P>0.05,圖6C)風(fēng)險評分差異無統(tǒng)計學(xué)意義。Ⅲ、Ⅳ期的風(fēng)險評分高于Ⅰ、Ⅱ期,且差異具有統(tǒng)計學(xué)意義(P<0.001,圖6D),T3~4期的風(fēng)險評分顯著高于T1~2期(P<0.01,圖6E),N1~3期的風(fēng)險評分顯著高于N0期(P<0.05,圖6F)。上述結(jié)果表明,風(fēng)險評分納入的ARGs與OSCC患者腫瘤分期、T分期及N分期顯著相關(guān)。
圖6 OSCC患者組間風(fēng)險評分分布Figure 6 Risk score distribution among OSCC patients
為了定量估計臨床環(huán)境下OSCC患者的生存概率,我們在TCGA訓(xùn)練集中構(gòu)建了一個列線圖,整合ARGs風(fēng)險評分模型和臨床病理特征(年齡、性別、病理分級、腫瘤分期、T分期、N分期)來預(yù)測患者1、3、5年的生存概率。通過圖7A,我們能夠定量得出列線圖的預(yù)測效果,樣本分?jǐn)?shù)越高,預(yù)后越差。此外,我們構(gòu)建了校準(zhǔn)曲線,以圖形化評估列線圖預(yù)測與實際生存率之間的一致性。在數(shù)據(jù)集的1年(圖7B)、3年(圖7C)、5年(圖7D)的生存概率中,實際生存和預(yù)測生存是吻合的。
圖7 構(gòu)建列線圖預(yù)測OSCC患者1、3、5年生存率Figure 7 Nomogram is constructed to predict 1-,3-,and 5-year survival rate of OSCC patients
如圖8A所示,利用CIBERSORT算法檢測到OSCC患者中22種免疫細(xì)胞的比例。通過Wilcoxon秩和檢驗比較高、低風(fēng)險組免疫細(xì)胞浸潤的差異,發(fā)現(xiàn)高風(fēng)險組M0期巨噬細(xì)胞顯著增加(P<0.001),CD8+T細(xì)胞、M1期巨噬細(xì)胞顯著低于低風(fēng)險組(P<0.001)。其中CD8+T細(xì)胞在高、低風(fēng)險組中浸潤差異如圖8B所示,高風(fēng)險組CD8+T細(xì)胞浸潤比例中位值為4.26%,低風(fēng)險組中位值為7.94%(P<0.001)。
圖8 免疫細(xì)胞浸潤模式Figure 8 Immune cell infiltration pattern
生物標(biāo)志物一直被認(rèn)為是癌癥診斷和靶向治療的重要工具。在過去的幾十年里,大量的生物標(biāo)志物被用于OSCC的早期檢測或治療過程中[15]。目前已經(jīng)應(yīng)用的OSCC患者靶向治療方案包括表皮生長 因 子受 體 (epidermal growth factor receptor,EGFR)抑制劑和西妥昔單抗(cetuximab,Cet),但對EGFR抑制劑治療的耐藥性是臨床使用該藥物的一大阻礙。此外,低成本和高產(chǎn)量的標(biāo)準(zhǔn)生物標(biāo)志物診斷模型仍然缺乏[16]。越來越多的研究表明,自噬在OSCC的腫瘤發(fā)生、進展和耐藥中起著重要作用[17-21]。一些ARGs被認(rèn)為是OSCC的潛在生物標(biāo)志物,如Liang等[22]指出,高表達的細(xì)胞質(zhì)自噬蛋白p62可能是OSCC患者預(yù)后不良的標(biāo)志物。Hou等[23]構(gòu)建了一個基于13個ARGs的OSCC預(yù)后模型,包括4個風(fēng)險基因和9個保護性基因,在不同數(shù)據(jù)集的驗證中均展現(xiàn)了良好的預(yù)測性能。本研究在保證預(yù)測準(zhǔn)確性的基礎(chǔ)上,納入了更少的ARGs構(gòu)建預(yù)后模型,可能在未來的前瞻性研究和臨床實際應(yīng)用中有更好的可操作性;此外,本研究還引入了免疫細(xì)胞浸潤分析,深層次地挖掘該模型影響預(yù)后的潛在機制。
本研究分析了OSCC和正常口腔組織中差異表達的ARGs,GO和KEGG富集分析表明,差異表達的ARGs主要富集在凋亡中,提示這些ARGs與凋亡顯著相關(guān),與之前的研究一致,自噬介導(dǎo)的細(xì)胞凋亡通過AKT/mTOR通路促進OSCC進展[24-25]。這些ARGs主要參與鉑耐藥途徑,其可能在OSCC治療策略中發(fā)揮重要作用[26]。另外,這些ARGs富集于ErbB信號通路和EGFR酪氨酸激酶抑制劑耐藥,這些都與OSCC進展相關(guān)[27],提示這些ARGs可能為OSCC提供潛在的治療靶點。通過單因素Cox回歸分析及LASSO回歸分析,將篩選出的6個ARGs(BID、NKX2-3、DDIT3、VEGFA、FADD、BIRC5)納 入風(fēng)險評分模型,并在TCGA數(shù)據(jù)集內(nèi)部和GEO驗證集中進行了評價。這6個ARGs中大多與OSCC或其他惡性腫瘤的發(fā)生和預(yù)后密切相關(guān),有研究表明,在二乙基亞硝胺(diethylnitrosamine,DEN)誘導(dǎo)的原發(fā)性肝癌模型中,通過BID基因缺失抑制肝細(xì)胞死亡途徑可以保護動物免受腫瘤發(fā)生[28];BID蛋白是接受輔助治療的結(jié)腸癌患者的獨立預(yù)后變量[29]。惡性轉(zhuǎn)化風(fēng)險高的口腔白斑與NKX2-3基因的啟動子異常甲基化有關(guān)[30];在人黑色素瘤細(xì)胞系中能夠觀察到NKX2-3基因的高甲基化[31]。DDIT3通過上調(diào)胃癌中的CEBPβ促進腫瘤干細(xì)胞(cancer stem cells,CSCs)的干性,并與預(yù)后不良有關(guān)[32];DDIT3基因在內(nèi)質(zhì)網(wǎng)應(yīng)激誘導(dǎo)的細(xì)胞凋亡和自噬中起重要作用[33-34]。VEGFA基因在OSCC組織中過表達,并且與患者較差的總生存率相關(guān)[35];同時,血管生成增加往往會導(dǎo)致遠(yuǎn)處轉(zhuǎn)移和預(yù)后不良[36]。FADD基因拷貝數(shù)和蛋白表達與中國臺灣地區(qū)OSCC患者的淋巴結(jié)轉(zhuǎn)移密切相關(guān)[37]。BIRC5/Survivin蛋白的細(xì)胞質(zhì)表達與OSCC患者的總體低生存率有關(guān),而核表達與較高的增殖率相關(guān)[38]。綜上所述,風(fēng)險評分模型中納入的預(yù)后ARGs主要通過調(diào)節(jié)自噬和凋亡、內(nèi)質(zhì)網(wǎng)應(yīng)激、血管生成等生物學(xué)過程參與腫瘤的發(fā)生、發(fā)展。同時,單因素、多因素Cox回歸分析提示風(fēng)險評分可作為獨立的預(yù)后因素,ROC分析也表明,與其他臨床危險因素相比,ARGs風(fēng)險評分模型具有更好的預(yù)測性能。本研究還分析了風(fēng)險評分與臨床病理變量之間的相關(guān)性,風(fēng)險評分的高低與腫瘤分期、T分期、N分期密切相關(guān),但它們之間的因果關(guān)系仍需進一步探索。為了準(zhǔn)確預(yù)測OSCC患者的生存率,本研究將風(fēng)險評分與多種臨床風(fēng)險因素(年齡、性別、病理分級、腫瘤分期、T分期、N分期)結(jié)合,建立列線圖用于生存概率的定量評估,校準(zhǔn)曲線也證明了預(yù)測生存率與實際生存率之間的一致性。
最后,本研究分析了高、低風(fēng)險組間免疫浸潤細(xì)胞的差異,揭示高風(fēng)險組預(yù)后差的潛在免疫機制。如,高風(fēng)險組CD8+T細(xì)胞顯著低于低風(fēng)險組(P<0.001),CD8+T細(xì)胞對于細(xì)胞內(nèi)病原體和腫瘤的保護性免疫很重要,并與患者免疫治療效果相關(guān),CD8+T細(xì)胞減少往往預(yù)示著預(yù)后不佳[39]。結(jié)合圖1E中KEGG富集的信息,OSCC中差異表達的ARGs與PD-L1表達和PD-1檢查點通路相關(guān),T細(xì)胞表面表達的PD-1與腫瘤細(xì)胞表面表達的PD-L1結(jié)合會向T細(xì)胞傳遞負(fù)向調(diào)控信號,抑制CD8+T細(xì)胞的產(chǎn)生和其破壞腫瘤細(xì)胞的能力[40]。因此,在本研究中,我們推測高風(fēng)險組差異表達的ARGs可能參與調(diào)控PD-1和PD-L1的免疫檢查點通路進而影響了CD8+T細(xì)胞免疫浸潤,最終導(dǎo)致了不良的預(yù)后。另外,Carleton等[41]發(fā)現(xiàn)自噬失活A(yù)tg5-/-小鼠的腫瘤浸潤淋巴細(xì)胞 (tumor infiltrating lymphocytes,TIL)中超過80%的CD8+TIL獲得效應(yīng)記憶表型,產(chǎn)生更多的干擾素γ(interferon-gamma,IFN-γ)和腫瘤壞死因子(tumor necrosis factor,TNF),而CD4+TIL效應(yīng)記憶無變化,提示自噬缺失增強了CD8+T細(xì)胞效應(yīng)記憶及功能;并發(fā)現(xiàn)自噬失活可阻止代謝下調(diào)H3K4me3水平,從而使免疫應(yīng)答基因啟動子的H3K4me3表達量得以維持。據(jù)此我們也可大膽假設(shè)本研究納入的預(yù)后相關(guān)基因可能通過同樣或類似的表觀遺傳修飾的生物學(xué)過程影響OSCC患者的免疫應(yīng)答及CD8+T的產(chǎn)生。這些猜測需要進一步的實驗加以論證,同時也為接下來研究自噬調(diào)控腫瘤發(fā)生、發(fā)展的潛在機制提供指引。
綜上所述,本研究綜合分析了ARGs,并且成功構(gòu)建了基于6個預(yù)后相關(guān)ARGs的風(fēng)險評分模型和列線圖,其可以為OSCC患者的預(yù)后提供準(zhǔn)確可靠的預(yù)測方法,幫助臨床醫(yī)生優(yōu)化個性化治療策略,同時也為OSCC腫瘤進展的潛在分子機制提供了一些思路。然而,本研究也存在一些局限性,首先,此次研究是回顧性的,因此風(fēng)險評分模型和列線圖需要在前瞻性研究和多中心臨床試驗中進一步驗證;另外,本研究所鑒定的ARGs特異性生物學(xué)特性及其在OSCC腫瘤發(fā)生、發(fā)展中的潛在機制有待進一步的實驗研究。