黎蘊(yùn)琪,蔣利和,陳闖*(. 廣西醫(yī)科大學(xué)附屬腫瘤醫(yī)院,南寧 53002;2. 右江民族醫(yī)學(xué)院,廣西 百色533000;3. 廣西大學(xué),南寧 530004)
結(jié)直腸癌(colorectal cancer,CRC)是世界上最常見的消化系統(tǒng)惡性腫瘤之一。據(jù)統(tǒng)計(jì),在全球范圍內(nèi),CRC發(fā)病率排第三,病死率排第二[1]。結(jié)腸腺癌(colorectal adenocarcinoma,COAD)是CRC中最常見的病理分型。目前,根治性手術(shù)治療是早期COAD最主要的治療方法,而對(duì)中晚期患者來(lái)說(shuō),化療是最常用且有效的選擇,可限制腫瘤細(xì)胞的擴(kuò)散和轉(zhuǎn)移,進(jìn)而控制病情。但化療藥物耐藥性的產(chǎn)生和腫瘤的高復(fù)發(fā)率仍是COAD臨床實(shí)踐中的主要問(wèn)題[2-3]。
近年來(lái),生物信息學(xué)在挖掘癌癥生物標(biāo)志物[4]和抗腫瘤藥物分子靶點(diǎn)[5]等研究方面發(fā)揮重要作用。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)主要運(yùn)用于分析大樣本基因表達(dá)數(shù)據(jù),將基因根據(jù)表達(dá)模式相似性分為不同模塊[6],這種方法跟普通聚類方法的不同在于其將基因表達(dá)值的相關(guān)系數(shù)取了N次冪,從而使得相關(guān)系數(shù)分布更加符合無(wú)尺度網(wǎng)絡(luò)分析,這更符合生物學(xué)規(guī)律[7]。相比于只關(guān)注差異表達(dá)的基因,WGCNA充分利用信息,把基因與表型的關(guān)聯(lián)轉(zhuǎn)換為基因集與表型的關(guān)聯(lián),免去了多重假設(shè)檢驗(yàn)校正的問(wèn)題。本研究基于WGCNA結(jié)合差異表達(dá)基因(DEGs)的方法,構(gòu)建高性能的預(yù)后預(yù)測(cè)模型,同時(shí)對(duì)這些靶點(diǎn)進(jìn)行化療藥物敏感性分析,旨在篩選出高風(fēng)險(xiǎn)腫瘤患者并對(duì)其進(jìn)行精準(zhǔn)的靶向治療,以達(dá)到個(gè)體COAD治療療效最大化。
微陣列數(shù)據(jù)集GSE39582和GSE41258以及這些數(shù)據(jù)集的相應(yīng)臨床數(shù)據(jù)從GEO數(shù)據(jù)庫(kù)(https://www.ncbi.nlm.nih.gov/geo)中下載,通過(guò)R軟件進(jìn)行數(shù)據(jù)歸一化處理,利用limma包去除批次效應(yīng),通過(guò)ggord包繪制PCA圖。GEO數(shù)據(jù)用于篩選DEGs,構(gòu)建預(yù)后模型。從癌癥基因組圖譜(TCGA)數(shù)據(jù)庫(kù)(http://cancergenome.nih.gov/)下載COAD的RNA-seq數(shù)據(jù)和相應(yīng)臨床信息。該數(shù)據(jù)集包括424個(gè)COAD樣本,用于驗(yàn)證預(yù)后模型和生存分析。
采用limma軟件包尋找mRNA的差異表達(dá),并通過(guò)pheatmap包繪制熱圖,“AdjustedP<0.05且|log2FC|>1”被定義為閾值,得到的DEGs被用于后續(xù)分析。
在評(píng)估GSE39582和GSE41258的基因表達(dá)譜后,使用WGCNA包構(gòu)建無(wú)標(biāo)度的共表達(dá)網(wǎng)絡(luò),計(jì)算基因的皮爾遜相關(guān)系數(shù)以獲得相關(guān)矩陣,選擇適當(dāng)?shù)能涢撝祦?lái)測(cè)量基因之間的連接性,將鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣(TOM),并計(jì)算相應(yīng)的差異(1-TOM)。具有相似表達(dá)模式的基因被分為同一個(gè)共表達(dá)模塊。此外,我們計(jì)算了模塊特征基因的差異性,對(duì)模塊進(jìn)行了層次聚類,并合并了相似的模塊(abline=0.25)。
將DEGs和關(guān)鍵模塊基因取交集,得到共同基因用于進(jìn)一步分析。將這些基因的表達(dá)譜和臨床信息合并,再隨機(jī)分為訓(xùn)練集(70%)和測(cè)試集(30%),使用訓(xùn)練集建立模型并在測(cè)試集進(jìn)行驗(yàn)證。在訓(xùn)練集中,通過(guò)單變量Cox比例風(fēng)險(xiǎn)回歸分析篩選與預(yù)后相關(guān)的基因,納入Lasso回歸分析,防止模型過(guò)度擬合,然后納入多因素Cox回歸分析,得出COAD預(yù)后風(fēng)險(xiǎn)模型,計(jì)算每個(gè)患者的風(fēng)險(xiǎn)評(píng)分,依據(jù)訓(xùn)練集風(fēng)險(xiǎn)評(píng)分的中位值,分別將訓(xùn)練集和測(cè)試集分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。運(yùn)用Kaplan-Meier生存分析和ROC曲線評(píng)估該預(yù)后風(fēng)險(xiǎn)模型的準(zhǔn)確性。利用生存時(shí)間和基因風(fēng)險(xiǎn)模型分別繪制散點(diǎn)圖和高低風(fēng)險(xiǎn)熱圖,利用風(fēng)險(xiǎn)評(píng)分結(jié)合臨床信息繪制列線圖。
為了進(jìn)一步確認(rèn)預(yù)后關(guān)鍵基因的潛在功能和識(shí)別高、低風(fēng)險(xiǎn)組的相關(guān)通路,使用GSEA軟件(4.0.1)進(jìn)行基因組富集分析,以P<0.05為閾值,展示KEGG相關(guān)通路的富集結(jié)果。
CellMiner是以美國(guó)國(guó)家癌癥研究所癌癥研究中心(NCI)所列出的60種癌細(xì)胞為基礎(chǔ)建立的數(shù)據(jù)庫(kù)[8]。從CellMiner數(shù)據(jù)庫(kù)下載NCI-60中基因mRNA表達(dá)數(shù)據(jù)和藥物活性相關(guān)數(shù)據(jù),探究5個(gè)關(guān)鍵基因的表達(dá)與藥物敏感性之間的關(guān)系。此外,通過(guò)數(shù)據(jù)庫(kù)GSE81005分析5-氟尿嘧啶(5-FU)耐藥的結(jié)直腸癌細(xì)胞系和野生型結(jié)直腸癌細(xì)胞系中HSD17B2mRNA的表達(dá)高低。
使用Wilcoxon檢驗(yàn)分析連續(xù)變量。使用Kaplan-Meier 生存分析和 Cox 單因素分析進(jìn)行預(yù)后分析。Pearson相關(guān)性分析用于評(píng)估5個(gè)關(guān)鍵基因的表達(dá)與藥物敏感性之間的相關(guān)性。使用R軟件(3.6.3版本)進(jìn)行統(tǒng)計(jì)分析,P<0.05為差異有統(tǒng)計(jì)學(xué)意義。
GSE39582和GSE41258數(shù)據(jù)集共包含73個(gè)正常樣本和752個(gè)COAD組織樣本。從這些樣本中鑒定出490個(gè)DEGs,其中上調(diào)基因313個(gè),下調(diào)基因177個(gè)(見圖1A、B)。數(shù)據(jù)批次效應(yīng)情況通過(guò)批次去除前后的PCA圖進(jìn)行評(píng)估(見圖1C、D)。
圖1 基因表達(dá)差異分析結(jié)果Fig 1 Gene expression differential analysis
為了找出與COAD相關(guān)的關(guān)鍵基因,筆者使用正常樣本和COAD樣本的基因表達(dá)矩陣進(jìn)行了WGCNA分析。將軟閾值β設(shè)置為5,平均連通性接近0(見圖2C、D)。將具有相似表達(dá)模式的DEGs聚集到相同的模塊中,最終得到11個(gè)共表達(dá)模塊(見圖2B),灰色模塊被認(rèn)為是無(wú)法被分配給任何模塊的基因集合。綠松石色模塊與正常樣本表型相關(guān)系數(shù)為0.64,綠色模塊與腫瘤樣本表型相關(guān)系數(shù)為0.53,通過(guò)分析基因與模塊的相關(guān)性(MM)和基因與臨床特征(正常和腫瘤)的相關(guān)性(GS)之間的關(guān)聯(lián),發(fā)現(xiàn)綠色和綠松石色模塊中MM和GS成正相關(guān)(見圖2E、F),說(shuō)明這些與性狀高度相關(guān)的基因,在關(guān)鍵模塊中也扮演著舉足輕重的角色。因此選擇綠色和綠松色為目的模塊,并將DEGs和目的模塊基因取交集篩選出247個(gè)基因用于進(jìn)一步分析。
圖2 WGCNA分析Fig 2 WGCNA analysis
將表達(dá)和生存數(shù)據(jù)合并后的740個(gè)GEO數(shù)據(jù)集樣本分為訓(xùn)練集(70%)和測(cè)試集(30%)。首先使用訓(xùn)練集的生存數(shù)據(jù)對(duì)247個(gè)DEGs進(jìn)行單因素Cox比例風(fēng)險(xiǎn)回歸分析,鑒定出7個(gè)對(duì)預(yù)后有顯著影響的基因(P<0.01),分別為GZMB、HSD17B2、PDK4、PIGR、PXMP2、GINS1、EMP1。將上述基因納入Lasso回歸分析降維(見圖3A、B),進(jìn)一步使用多變量Cox比例風(fēng)險(xiǎn)回歸分析,最終得到5個(gè)與COAD預(yù)后相關(guān)的風(fēng)險(xiǎn)基因(見圖3C)。基于這5個(gè)關(guān)鍵基因進(jìn)行風(fēng)險(xiǎn)預(yù)后模型的構(gòu)建:Risk score=-0.162GZMB+0.210HSD17B2+0.184PDK4-0.141PIGR-0.516PXMP2。生存分析結(jié)果均顯示低風(fēng)險(xiǎn)組的患者生存狀況明顯優(yōu)于高風(fēng)險(xiǎn)組(P<0.001,見圖4A、B)。訓(xùn)練集1年、2年和3年的隨訪預(yù)測(cè)AUC分別為0.690、0.709、0.699(見圖4C),測(cè)試集的隨訪預(yù)測(cè)AUC,見圖4D。從生存時(shí)間和風(fēng)險(xiǎn)評(píng)分繪制的熱圖和散點(diǎn)圖中看出,隨著風(fēng)險(xiǎn)得分的增加,死亡的患者也增加,存活時(shí)間相對(duì)減少,由此可見模型有相對(duì)較好的預(yù)測(cè)能力(見圖5)。
圖3 LASSO及多變量Cox回歸分析進(jìn)一步篩選預(yù)后基因Fig 3 Further screening of prognostic genes by the LASSO and multivariate Cox regression analysis
圖4 訓(xùn)練集(A、C)、測(cè)試集(B、D)Kaplan-Meier曲線和時(shí)間依賴性ROC曲線Fig 4 Distribution of time-dependent ROC curves and Kaplan-Meier curves of the training set(A and C)and the test set(B and D)
圖5 訓(xùn)練集(A)、測(cè)試集(B)風(fēng)險(xiǎn)評(píng)分的分布、患者的生存狀態(tài)和5個(gè)基因的表達(dá)熱圖Fig 5 Distribution of risk scores,patient survival status,and the five-gene expression heat map of the training set(A)and the test set(B)
進(jìn)一步通過(guò)單因素和多因素Cox回歸分析評(píng)估COAD患者風(fēng)險(xiǎn)評(píng)分和臨床病理信息與總生存期的關(guān)系。訓(xùn)練集和測(cè)試集的結(jié)果均顯示風(fēng)險(xiǎn)評(píng)分可以作為獨(dú)立危險(xiǎn)因素來(lái)預(yù)測(cè)COAD患者的生存預(yù)后(P<0.01),多因素Cox回歸分析結(jié)果顯示,風(fēng)險(xiǎn)評(píng)分在訓(xùn)練集(見圖6B)和測(cè)試集(見圖6D)中的HR分別為1.581和1.616。
圖6 訓(xùn)練集(A、B)、測(cè)試集(C、D)關(guān)于臨床特征因素和風(fēng)險(xiǎn)評(píng)分的獨(dú)立預(yù)后分析Fig 6 Independent prognostic analysis of clinical characteristics and risk score of both the training set(A and B)and the test set(C and D)
接下來(lái),在TCGA驗(yàn)證集(424例樣本)中使用上述模型公式計(jì)算每位患者的風(fēng)險(xiǎn)評(píng)分,將患者分為高風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)組,并通過(guò)KM曲線比較預(yù)后。結(jié)果表明,該模型可以區(qū)分患者的預(yù)后(P=0.013,見圖7)。
圖7 TCGA驗(yàn)證集的生存分析Fig 7 Survival on TCGA validation set
筆者使用臨床病理特征(年齡、性別、腫瘤分級(jí)、T分期、N分期)和5個(gè)基因的表達(dá)在GEO測(cè)試集具有完整臨床資料的503例樣本中建立了建立列線圖來(lái)預(yù)測(cè)患者的生存率。結(jié)果顯示列線圖可以預(yù)測(cè)1年、2年和3年的生存率,接近理想模型(見圖8A)。校準(zhǔn)圖顯示出較佳的預(yù)測(cè)準(zhǔn)確性(見圖8B),并且預(yù)測(cè)的生存率與實(shí)際生存率基本重合。
圖8 預(yù)后風(fēng)險(xiǎn)模型列線圖(A)及校準(zhǔn)曲線圖(B)Fig 8 Nomogram(A)and calibration(B)curve of the prognostic risk model
為了探索5個(gè)關(guān)鍵基因與COAD的潛在生物學(xué)機(jī)制,筆者使用GSEA在訓(xùn)練集和測(cè)試集中對(duì)高、低風(fēng)險(xiǎn)組之間的DEGs進(jìn)行了功能富集分析。訓(xùn)練集和測(cè)試集結(jié)果均顯示低風(fēng)險(xiǎn)組的基因集在代謝以及細(xì)胞周期等相關(guān)途徑有明顯的富集作用(見圖9),提示低風(fēng)險(xiǎn)組可能通過(guò)細(xì)胞周期通路影響腫瘤細(xì)胞凋亡進(jìn)而影響患者生存。
圖9 訓(xùn)練集(A)、測(cè)試集(B)GSEA富集分析Fig 9 GSEA enrichment analysis of the training set(A)and the test set(B)
最后,筆者研究了5個(gè)關(guān)鍵基因在60種人類癌細(xì)胞系(NCI-60)中的表達(dá)水平及其與化療藥物之間的相關(guān)性。這些基因和多種腫瘤的化療藥物的耐藥性增加有關(guān)(P<0.05),其中HSD17B2mRNA的表達(dá)增加與多種治療COAD的化療藥物的耐藥性增加有關(guān)(見圖10A~E),包括伊立替康(cor=-0.262,P=0.043,中晚期CRC的一線用藥[9])、5-氟脫氧尿苷(cor=-0.318,P=0.013)、氟尿苷(cor=-0.334,P=0.009,治療結(jié)腸癌的肝轉(zhuǎn)移[10])、依托泊苷(cor=-0.308,P=0.017)、雷替曲塞(cor=-0.298,P=0.021)。另外,通過(guò)數(shù)據(jù)庫(kù)GSE81005對(duì)5-FU耐藥的CRC細(xì)胞系和野生型CRC細(xì)胞系進(jìn)行分析,發(fā)現(xiàn)5-FU耐藥的CRC細(xì)胞系中HSD17B2mRNA的表達(dá)明顯高于野生型CRC細(xì)胞系(P<0.001,見圖10F)。這些結(jié)果表明HSD17B2在CRC的發(fā)生發(fā)展和耐藥性中發(fā)揮著重要生物學(xué)作用。
圖10 HSD17B2的表達(dá)與伊立替康(A)、5-氟脫氧尿苷(B)、氟尿苷(C)、依托泊苷(D)和雷替曲塞(E)的藥物敏感性的相關(guān)性分析以及HSD17B2在CRC野生型細(xì)胞及其5-FU誘導(dǎo)耐藥細(xì)胞系中的表達(dá)(F)Fig 10 Correlation between HSD17B2 expression and drug sensitivity of irinotecan(A),5-fluorodeoxyuridine(B),fluorouridine(C),etoposide(D)and raltitrexed(E); expression of HSD17B2 in CRC wild type cells and its 5-FU-induced resistant cell lines(F)
COAD是全球范圍內(nèi)最常見的癌癥類型之一,惡性程度和病死率較高[1,11]。目前COAD的治療主要采取手術(shù)切除及化療等傳統(tǒng)治療方法,對(duì)于早期腫瘤效果較好,但對(duì)于中、晚期腫瘤只能延長(zhǎng)患者生命,緩解臨床癥狀,不能從根本上提高患者的生活質(zhì)量,患者預(yù)后仍較差。因此,尋找COAD中用于生存評(píng)估和靶向治療的生物標(biāo)志物仍至關(guān)重要。
本研究通過(guò)分析兩個(gè)GEO數(shù)據(jù)集來(lái)鑒定DEGs,構(gòu)建了WGCNA共表達(dá)網(wǎng)絡(luò),從中找出5個(gè)預(yù)后關(guān)鍵基因(GZMB、HSD17B2、PDK4、PIGR和PXMP2),構(gòu)建COAD預(yù)后預(yù)測(cè)模型,并對(duì)其進(jìn)行了驗(yàn)證,列線圖和ROC曲線顯示該模型具有較好的臨床實(shí)用潛力。藥物敏感性的分析顯示,這些關(guān)鍵基因與多種化療藥物的耐藥性相關(guān)。先前的研究表明,顆粒酶(GZMs)是細(xì)胞毒性T淋巴細(xì)胞和自然殺傷細(xì)胞的關(guān)鍵效應(yīng)分子,在控制細(xì)胞內(nèi)病原體和癌癥方面發(fā)揮著關(guān)鍵作用[12-13]。顆粒酶B(GZMB)作為抗腫瘤和抗感染因子,在多種腫瘤中顯著表達(dá)并與良好預(yù)后相關(guān),如非小細(xì)胞肺癌、CRC、骨肉瘤等[14-16]。17-β羥基類固醇脫氫酶2型(HSD17B2)是一種參與雌激素和雄激素調(diào)節(jié)的酶[17],其可催化睪酮和雌二醇(E2)分別轉(zhuǎn)化為雄烯二酮和雌酮(E1)[18]。而雌二醇被報(bào)道可減少CRC細(xì)胞的增殖和促進(jìn)凋亡[19]。本研究構(gòu)建的風(fēng)險(xiǎn)模型也提示,隨著HSD17B2的表達(dá)增加,COAD患者的預(yù)后更差,這可能與雌二醇的失活相關(guān)。Lee等[20]發(fā)現(xiàn)HSD17B2可以作為獨(dú)立預(yù)測(cè)因子,其過(guò)表達(dá)與直腸癌患者的不良預(yù)后相關(guān)。丙酮酸脫氫酶激酶4(PDK4)是細(xì)胞能量代謝的關(guān)鍵調(diào)節(jié)劑,可以參與m6A調(diào)節(jié)的糖酵解和ATP生成[21]。研究表明PDK4的高表達(dá)可以促進(jìn)腫瘤細(xì)胞的增殖、遷移和侵襲,例如CRC、宮頸癌、肝癌、膀胱癌細(xì)胞等[21-23]。Woolbright等[23]發(fā)現(xiàn),抑制PDK4表達(dá)可能會(huì)加劇順鉑誘導(dǎo)的細(xì)胞死亡。聚合免疫球蛋白受體(PIGR)是黏膜免疫系統(tǒng)的關(guān)鍵組成部分,在非小細(xì)胞肺癌、鼻咽癌、上皮性卵巢癌等的表達(dá)降低[24-27]。PIGR被證實(shí)了和消化道腫瘤密切相關(guān)。研究發(fā)現(xiàn),PIGR可以啟動(dòng)Yes-Dap12-Syk-Rac1/CDC42-MEK/ERK信號(hào)級(jí)聯(lián)以促進(jìn)腫瘤細(xì)胞的增殖和轉(zhuǎn)化[28],富含PIGR的細(xì)胞外囊泡可以促進(jìn)肝癌細(xì)胞的侵襲性[29]。Ohkuma等[30]發(fā)現(xiàn),PIGR的mRNA和蛋白水平是胰腺癌的獨(dú)立預(yù)后因素,其高表達(dá)與患者預(yù)后不良有關(guān)。PXMP2是高等真核生物中最豐富的一種同源三聚體過(guò)氧化物酶體膜蛋白(PMP),具有通道形成活性[31-32],參與廣泛的代謝途徑,包括脂質(zhì)、氨基酸和羥基酸、嘌呤和活性氧物質(zhì)的轉(zhuǎn)化[33]。但PXMP2基因與人類疾病的關(guān)系尚未被充分了解。
對(duì)高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組進(jìn)行GSEA功能富集分析,發(fā)現(xiàn)低風(fēng)險(xiǎn)組在細(xì)胞周期、DNA復(fù)制、RNA聚合酶和代謝相關(guān)通路顯著富集。細(xì)胞周期失調(diào)和細(xì)胞周期蛋白依賴性激酶(CDK)激活導(dǎo)致的持續(xù)細(xì)胞增殖是癌癥的標(biāo)志[34]。目前,已經(jīng)有3種選擇性CDK4/6抑制劑獲得美國(guó)食品藥品監(jiān)督管理局的批準(zhǔn)[35],多種CDK4/6抑制劑正處于治療癌癥的臨床試驗(yàn)中[36]。此外,通過(guò)藥物敏感性分析發(fā)現(xiàn)這5個(gè)關(guān)鍵基因與多種化療藥物相關(guān),特別是HSD17B2與結(jié)腸癌化療藥物(伊立替康、5-氟脫氧尿苷、依托泊苷、雷替曲塞等)存在顯著負(fù)相關(guān)。此外,5-FU耐藥的細(xì)胞系中的HSD17B2表達(dá)顯著高于野生組。5-FU是CRC治療中使用最廣泛的化療藥物之一。盡管近年越來(lái)越多的化療藥物出現(xiàn),但以5-FU為基礎(chǔ)的化療(如 XELOX 和 FOLFOX)仍然是目前中晚期 CRC 的一線化療藥物。但由于化療耐藥性,其總體反應(yīng)率僅為 10%~15%[37],化療耐藥和嚴(yán)重的不良反應(yīng)已成為影響化療治療效果的主要因素[38]。為了克服化學(xué)耐藥性和減少不良反應(yīng),許多藥物在臨床試驗(yàn)中進(jìn)行了測(cè)試,但迄今為止尚未取得重大進(jìn)展,需要更多努力尋找能夠克服5-FU耐藥性且毒性較小的新藥以改善CRC患者。本研究提示HSD17B2可能是5-FU化療耐藥的潛在生物標(biāo)志物,針對(duì)HSD17B2的靶向治療或許能在提高化療敏感性的同時(shí)改善患者的預(yù)后??傊?,預(yù)后模型中的5個(gè)靶基因可能通過(guò)影響細(xì)胞周期和細(xì)胞對(duì)化療藥物的敏感性,進(jìn)而影響COAD的進(jìn)展和預(yù)后。
綜上,本研究篩選的5個(gè)靶基因是潛在COAD的預(yù)后生物標(biāo)志物,HSD17B2更是有望作為降低COAD化療藥物耐藥性的靶點(diǎn)。然而,本研究是基于多個(gè)公共數(shù)據(jù)庫(kù)的數(shù)據(jù)分析,具有局限性,這些關(guān)鍵基因在結(jié)腸腺癌預(yù)后以及化療耐藥中的分子機(jī)制仍需要進(jìn)一步的基礎(chǔ)實(shí)驗(yàn)來(lái)證實(shí)。
本研究篩選了5個(gè)樞紐基因可能為COAD的潛在預(yù)后標(biāo)志物,并發(fā)現(xiàn)HSD17B2可能是化療耐藥的潛在生物標(biāo)志物,為進(jìn)一步研究COAD的發(fā)病機(jī)制和治療提供潛在靶點(diǎn)。