WGCNA與機(jī)器學(xué)習(xí)算法識別結(jié)直腸癌發(fā)病機(jī)制中潛在PANoptosis關(guān)鍵基因

李文靜,楊杰,陳冬梅,2

(1.貴州醫(yī)科大學(xué),貴州 貴陽 550000; 2.六盤水市人民醫(yī)院 消化內(nèi)科,貴州 六盤水 553000)

結(jié)直腸癌(CRC)作為一種常見的消化道惡性腫瘤,其發(fā)病率和死亡率位居腫瘤第3位和第2位[1]。盡管CRC的5年生存率已有明顯提高,但目前中國CRC的發(fā)病率和死亡率仍呈上升趨勢且有年輕化傾向[2]。由于CRC的早期癥狀無特異性,易被誤診為痔出血、息肉出血等,部分患者在確診時已處于中晚期。因此,早期診斷是改善CRC患者預(yù)后的關(guān)鍵。結(jié)腸鏡檢查是診斷CRC的最典型方法,具有較高特異性和敏感性,但為侵入性檢查,具有價格高、條件苛刻等缺點(diǎn),患者接受性差[3]。以癌胚抗原(CEA)為代表的非侵入性生物標(biāo)志物診斷CRC的特異性和敏感性不理想[4]。篩選高敏感性和特異性的生物標(biāo)志物對于CRC的篩查、診療和預(yù)后預(yù)測具有重要意義。

PANoptosis是近來定義的一個新概念,強(qiáng)調(diào)細(xì)胞凋亡、焦亡和壞死性凋亡之間的串?dāng)_和協(xié)調(diào)[5]。既往針對CRC與各種形式細(xì)胞死亡相關(guān)的基因已經(jīng)在許多研究中被報道,包括自噬、焦亡、壞死性凋亡相關(guān)基因。然而,與不同類型細(xì)胞死亡之間串?dāng)_和協(xié)調(diào)相關(guān)基因的研究仍不充分。生物信息學(xué)是結(jié)合了生物學(xué)、數(shù)學(xué)和信息技術(shù)的新興學(xué)科。利用生物信息學(xué)方法篩選腫瘤相關(guān)分子生物標(biāo)志物在潛在的高?；颊咴\斷、鑒定和預(yù)后評估中顯示出卓越的性能[6]。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis,WGCNA)作為生物信息學(xué)分析中被高頻應(yīng)用的一種方法,在篩選疾病特征生物標(biāo)志物和潛在防治靶點(diǎn)中表現(xiàn)出明顯優(yōu)勢,在包括腫瘤在內(nèi)的多種疾病中應(yīng)用廣泛[7]。機(jī)器學(xué)習(xí)算法作為人工智能的子集之一,被廣泛應(yīng)用于疾病生物標(biāo)志物、發(fā)病機(jī)制、治療靶點(diǎn)等方面的研究[8]。本研究基于高通量測序數(shù)據(jù),利用WGCNA與機(jī)器學(xué)習(xí)算法相結(jié)合的方法,篩選與CRC發(fā)病相關(guān)的潛在PANoptosis關(guān)鍵基因并構(gòu)建疾病預(yù)測模型,以期為CRC的診斷和治療提供一定指導(dǎo)依據(jù)。

1 資料與方法

1.1 數(shù)據(jù)集的下載及數(shù)據(jù)處理

從美國國立生物技術(shù)信息中心(NCBI)的GEO芯片數(shù)據(jù)庫中選擇了GSE39582和GSE74602 2個芯片數(shù)據(jù)集,從中獲取CRC和正常黏膜組織的基因表達(dá)譜。GSE39582采用美國昂菲公司人基因組U133 Plus 2.0芯片GPL570平臺,包括443個腫瘤樣本和19個對照樣本。GSE74602采用美國昂菲公司人基因組Ref8 2.0芯片GPL6104平臺,包括30個腫瘤樣本和30個對照樣本。下載芯片數(shù)據(jù)壓縮包和探針文件,通過R4.2.3軟件的RMA算法對芯片數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理。采用|logFC|>1和adj.p<0.05作為篩選標(biāo)準(zhǔn),使用R軟件的limma包進(jìn)行差異表達(dá)基因(DEGs)的篩選。

1.2 功能富集分析

分別采用R軟件enrichDO函數(shù)、clusterProfiler包對DEGs進(jìn)行疾病本體(DO)、基因本體(GO)和基因組京都百科全書(KEGG)富集分析。其中GO分析包括了差異基因的細(xì)胞組成(CC)、分子功能(MF)和生物學(xué)過程(BP)。結(jié)果根據(jù)-log10(P)值大小進(jìn)行排序,同時呈現(xiàn)富集基因占總體的比例。

1.3 WGCNA

采用R軟件WGCNA包進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)識別hub基因。首先,計算所有基因的表達(dá)相關(guān)系數(shù),構(gòu)建相似矩陣。將相似矩陣轉(zhuǎn)化為拓?fù)渚仃?采用拓?fù)渲丿B度量(TOM)描述基因之間的關(guān)聯(lián)度。根據(jù)基因間相關(guān)性確定軟閾值(power),構(gòu)建基因聚類樹。根據(jù)基因表達(dá)相關(guān)系數(shù)確定表達(dá)模塊,每個模塊內(nèi)的最少基因數(shù)設(shè)定為100個,模塊間相關(guān)系數(shù)設(shè)定為0.85。特征基因由模塊內(nèi)基因連通性決定。根據(jù)基因表達(dá)趨勢計算皮爾遜相關(guān)系數(shù),取相關(guān)系數(shù)最大的模塊內(nèi)基因作為特征基因。

1.4 機(jī)器學(xué)習(xí)算法

利用R軟件的venn包取DEGs和WGCNA特征基因的交集基因。利用R軟件glmnet包的最小絕對值收斂和選擇算子(LASSO)算法、e1071包的支持向量機(jī)-遞歸特征消除(SVM-RFE)、RandomForest包的隨機(jī)森林算法從交集基因中獲得樞紐基因。

1.5 Panoptosis關(guān)鍵基因的篩選

基于geneCards(https:∥www.genecards.org/),分別以“Pyroptosis”“Necroptosis”和“Apoptosis”為關(guān)鍵詞獲得焦亡、壞死和凋亡相關(guān)基因。利用R軟件的venn包取焦亡、壞死和凋亡相關(guān)基因的交集,獲得PANoptosis相關(guān)基因(PANRGs)。利用R軟件的venn包,取PANoptosis相關(guān)基因和樞紐基因的交集,獲得PANoptosis關(guān)鍵基因。

1.6 診斷列線圖模型的構(gòu)建和驗證

以GSE39582作為訓(xùn)練集,GSE74602為驗證集,基于PANoptosis關(guān)鍵基因在訓(xùn)練集和驗證集中的表達(dá)量,使用R軟件rms包構(gòu)建預(yù)測CRC的列線圖模型。使用校準(zhǔn)曲線評估列線圖模型的預(yù)測能力。采用決策曲線評估列線圖模型的臨床效用。使用R包pROC函數(shù)創(chuàng)建受試者工作特征(ROC)曲線,以確定PANoptosis關(guān)鍵基因和列線圖模型在訓(xùn)練集和驗證集中對CRC的診斷價值。

2 結(jié) 果

2.1 DEGs篩選結(jié)果與功能富集分析

根據(jù)adj.p<0.05和|log2FC|>1篩選標(biāo)準(zhǔn),從GSE39582中篩選出1 503個DEGs,其中上調(diào)744個,下調(diào)759個,見圖1A。DO分析(圖1B)顯示,DEGs主要與良性腫瘤、腎細(xì)胞癌、乳腺癌、腸道疾病等多種腫瘤有關(guān)。GO分析結(jié)果(圖1C)表明,DEGs主要參與核分裂、染色體分離、有絲分裂的核分裂等生物學(xué)過程;主要定位于含膠原蛋白的細(xì)胞外基質(zhì)、細(xì)胞的頂端部分、頂端質(zhì)膜等細(xì)胞成分;主要參與細(xì)胞外基質(zhì)結(jié)構(gòu)成分、糖胺聚糖結(jié)合、細(xì)胞因子活性、CXCR趨化因子受體結(jié)合等分子功能。KEGG分析(圖1D)顯示,DEGs主要富集于細(xì)胞周期、嘌呤代謝、TNF信號通路、PPAR信號通路、P53信號通路等相關(guān)信號通路。

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

使用R軟件“WGCNA”包對GSE39582中的樣本進(jìn)行聚類(圖2A)。隨后進(jìn)行無標(biāo)尺度指數(shù)和平均連通性分析,當(dāng)軟閾值β=4時,網(wǎng)絡(luò)達(dá)到無標(biāo)尺度拓?fù)溟撝?.9,且網(wǎng)絡(luò)結(jié)點(diǎn)平均連接性高(圖2B)。通過動態(tài)樹切割和計算,獲得8個基因模塊。將8個模塊與數(shù)據(jù)集對照組和腫瘤組之間進(jìn)行相關(guān)性分析,得到相關(guān)性熱圖(圖2C),選擇與腫瘤相關(guān)性最強(qiáng)的棕色模塊(r=-0.65,P<0.001)為CRC的關(guān)鍵模塊。根據(jù)基因顯著性的過濾條件為0.5,以及基因與模塊相關(guān)性的過濾條件為0.8,確定棕色模塊中核心基因820個(圖2D)。

A.樣本-性狀聚類樹;B.軟閾值參數(shù)與尺度獨(dú)立性、平均連接性分析(最佳軟閾值為4);C.臨床性狀與基因模塊的相關(guān)性熱圖,單元格中的數(shù)字分別表示相關(guān)性系數(shù)r和對應(yīng)的P值;D.棕色模塊基因散點(diǎn)圖

2.3 特征差異基因的確定及功能富集分析

根據(jù)基因差異表達(dá)分析和WGCNA的結(jié)合,共有1 503個DEGs和820個關(guān)鍵模塊基因,二者取交集,得到525個交集基因(圖3)。對這525個交集基因進(jìn)行富集分析,DO富集分析共富集到422條DO條目,涉及腸道疾病、胃癌、結(jié)腸炎、炎癥性腸病。GO分析共富集到1 637條GO條目,這些交集基因在醇類代謝過程、四氯化碳運(yùn)輸、細(xì)胞頂端部分、頂端質(zhì)膜、陰離子跨膜轉(zhuǎn)運(yùn)體活性、次級活性跨膜轉(zhuǎn)運(yùn)體活性等方面顯著富集。KEGG分析表明,這些交集基因涉及緊密連接、視黃醇代謝、藥物代謝等途徑。

圖3 特征差異基因的確定及功能富集分析

2.4 基于3種不同機(jī)器學(xué)習(xí)算法篩選樞紐基因

在GSE39582中,利用SVM-RFE對上述525個特征差異基因進(jìn)行篩選,“最佳點(diǎn)”在0.003 39時5折交叉驗證誤差最小,共篩選出8個特征基因(圖4A)。隨機(jī)森林算法獲得30個特征基因(圖4B)。LASSO回歸分析顯示,在選擇最小均方誤差對數(shù)λ時,回歸模型中有9個系數(shù)不為零的特征基因(圖4C)。將SVM-RFE、RandomForest和LASSO分析所篩選的基因取交集,共得到6個樞紐基因(HubGenes,圖4D),即細(xì)胞周期蛋白依賴性激酶1(CDK1)、二肽酶1(DPEP1)、水通道蛋白8(AQP8)、溶質(zhì)載體家族30成員10(SLC30A10)、半胱氨酸天冬氨酸蛋白水解酶 7(CASP7)和半胱氨酸天冬氨酸蛋白水解酶 8(CASP8)。

A.SVM-RFE分析;B.Random Forest分析;C.LASSO回歸分析;D.SVM-RFE、Random Forest和LASSO交叉分析的韋恩圖

2.5 PANoptosis關(guān)鍵基因的篩選與功能富集分析

以“Pyroptosis”為關(guān)鍵詞,在geneCards中獲得焦亡相關(guān)基因436個。以“Necroptosis”為關(guān)鍵詞,在geneCards中獲得壞死相關(guān)基因667個。以“Apoptosis”為關(guān)鍵詞,在geneCards中獲得凋亡相關(guān)基因14 989個。取交集,得到PANoptosis相關(guān)基因(PANRGs)96個(圖5A)。進(jìn)一步取PANRGs和HubGenes的交集,獲得4個PANoptosis關(guān)鍵基因,即CDK1、DPEP1、CASP7和CASP8(圖5B)。

A.焦亡、壞死和凋亡相關(guān)基因的交集,獲得PANoptosis相關(guān)基因(PANRGs);B.取PANRGs和樞紐基因(HubGenes)的交集,獲得PANoptosis關(guān)鍵基因

2.6 CRC預(yù)測模型構(gòu)建與驗證

基于4個PANoptosis關(guān)鍵基因,在訓(xùn)練集GSE39582中采用Logistic回歸分析構(gòu)建nomogram圖(圖6A)。校準(zhǔn)圖顯示,校準(zhǔn)預(yù)測曲線與標(biāo)準(zhǔn)曲線貼合較好。決策(DCA)曲線顯示,列線圖模型在預(yù)測CRC發(fā)生的臨床效能上表現(xiàn)良好。在驗證集GSE74602中構(gòu)建的列線圖、校準(zhǔn)圖和DCA曲線也證實(shí)了上述結(jié)果(圖6B)。

A.訓(xùn)練集GSE39582中基于PANoptosis關(guān)鍵基因構(gòu)建的列線圖模型、校準(zhǔn)曲線和DCA曲線;B.驗證集GSE74602中基于PANoptosis關(guān)鍵基因構(gòu)建的列線圖模型、校準(zhǔn)曲線和DCA曲線

2.7 PANoptosis關(guān)鍵基因和預(yù)測模型的診斷效能評價

在訓(xùn)練集GSE39582中,CDK1、DPEP1、CASP7和CASP8的AUC分別為0.926、0.953、0.991和0.873,而預(yù)測模型的AUC為0.998。在驗證集GSE74602中,CDK1、DPEP1、CASP7和CASP8的AUC分別為0.999、0.931、0.941和0.651,而預(yù)測模型的AUC為1.000。說明預(yù)測模型能夠準(zhǔn)確地區(qū)分正常和腫瘤組織樣本。見圖7。

A.在訓(xùn)練集GSE39582中CDK1、DPEP1、CASP7、CASP8和列線圖模型區(qū)分腫瘤組織和正常組織的ROC曲線;B.在驗證集GSE74602中CDK1、DPEP1、CASP7、CASP8和列線圖模型區(qū)分腫瘤組織和正常組織的ROC曲線

3 討 論

CRC早期診斷困難,治療效果及預(yù)后均不理想,且詳細(xì)的分子生物學(xué)機(jī)制尚未完全闡明。探究新的分子標(biāo)志物和分子治療靶點(diǎn)對研究CRC的發(fā)病機(jī)制、預(yù)防和治療具有重要意義。近年來,微陣列和生物信息學(xué)分析已被廣泛用于鑒定不同癌癥診斷和治療的潛在靶標(biāo)[9]。本研究基于差異基因表達(dá)篩選與WGCNA共鑒定出525個疾病特征DEGs,主要與涉及腸道疾病、胃癌、結(jié)腸炎、炎癥性腸病等疾病。此外,它們可以通過調(diào)控醇類代謝過程、四氯化碳運(yùn)輸、藥物代謝等生物學(xué)功能和途徑介導(dǎo)CRC的病理過程。這些結(jié)果證實(shí)了這些重要DEGs介導(dǎo)CRC的分子機(jī)制與既往文獻(xiàn)報道結(jié)果存在一些吻合[10],在一定程度上證實(shí)了本研究結(jié)果的準(zhǔn)確性。

機(jī)器學(xué)習(xí)算法能夠從新視角來分析大群體基因測序或微陣列數(shù)據(jù)。隨機(jī)森林、LASSO和SVM等機(jī)器學(xué)習(xí)算法解決了既往高維數(shù)據(jù)構(gòu)建的預(yù)測模型可能具有的過度擬合風(fēng)險等問題[11]。與使用組學(xué)生物標(biāo)志物組合的數(shù)學(xué)方法相比,機(jī)器學(xué)習(xí)算法定義風(fēng)險等級在疾病預(yù)測中表現(xiàn)出更優(yōu)的性能[12]。因此,集合多種機(jī)器學(xué)習(xí)算法開發(fā)與CRC病理機(jī)制相關(guān)的生物標(biāo)志物來預(yù)測CRC有重要臨床意義。本研究利用LASSO、SVM-RFE和RandomForest算法從525個特征差異基因中鑒定出了6個CRC的樞紐基因,即CDK1、DPEP1、AQP8、SLC30A10、CASP7和CASP8。

腫瘤的發(fā)生發(fā)展是一個與多基因相關(guān)的多步驟、多階段復(fù)雜過程。目前研究認(rèn)為,腫瘤是細(xì)胞增殖和死亡失衡所致。其中,細(xì)胞死亡在腫瘤發(fā)生發(fā)展過程中起負(fù)調(diào)控作用,可以遏制腫瘤細(xì)胞的快速生長。在程序性細(xì)胞死亡的擬議形式中,細(xì)胞凋亡、焦亡和壞死性凋亡已被廣泛研究,這些過程涉及負(fù)責(zé)細(xì)胞死亡的啟動、轉(zhuǎn)導(dǎo)和執(zhí)行的復(fù)雜分子機(jī)制[13]。然而,關(guān)于細(xì)胞死亡的早期研究通常集中于這些獨(dú)立機(jī)制下的獨(dú)特程序和功能。PANoptosis是最近發(fā)現(xiàn)的一種程序性細(xì)胞死亡類型,可促進(jìn)腫瘤細(xì)胞死亡,并可能降低異常凋亡途徑對腫瘤化學(xué)抗性的影響[14]。作為一種依賴于PANoptosome復(fù)合物的炎癥程序性細(xì)胞死亡,PANoptosis含有焦亡、細(xì)胞凋亡和壞死性凋亡所必需的分子,可激活不同的細(xì)胞死亡途徑[15]。即使其他程序性細(xì)胞死亡途徑被抑制或突變,PANoptosis也可以提供另一種途徑來殺死癌細(xì)胞并降低癌細(xì)胞獲得耐藥性的風(fēng)險[16]。因此,探討關(guān)于PANoptosis對癌癥影響的機(jī)制對于開發(fā)癌癥治療新策略至關(guān)重要。Huang等[17]篩選了與CRC預(yù)后特征相關(guān)的PANoptosis相關(guān)基因,并提出基于PANoptosis相關(guān)基因構(gòu)建的風(fēng)險模型有望成為CRC的潛在預(yù)后標(biāo)志物。本研究進(jìn)一步取樞紐基因與PANoptosis相關(guān)基因的交集,得到4個與CRC相關(guān)的PANoptosis關(guān)鍵基因,即CDK1、DPEP1、CASP7和CASP8。CDK1是絲氨酸/蘇氨酸蛋白家族成員之一,已被證實(shí)參與細(xì)胞周期的調(diào)控。研究證明,CDK1可作為預(yù)測Ⅱ期CRC復(fù)發(fā)的指標(biāo),且其高表達(dá)與較高的遠(yuǎn)處轉(zhuǎn)移率和較低的病因特異性生存率顯著相關(guān)[18]。DPEP1是由肽鍵連接的兩個氨基酸殘基組成的一種二肽化合物,可在胰腺、乳腺、結(jié)直腸等多種組織中表達(dá)[19]。研究[20]顯示,DPEP1通過調(diào)節(jié)基因轉(zhuǎn)錄調(diào)控通路來發(fā)揮致癌作用。黃金華等[21]認(rèn)為,DPEP1在CRC不同病理階段的表達(dá)有助于CRC的早期診斷。CASP7是凋亡通路關(guān)鍵基因,在控制細(xì)胞死亡中充當(dāng)環(huán)境傳感器,其遺傳變異已被證實(shí)與CRC患者放化療敏感性和副反應(yīng)及生存密切相關(guān)[22]。CASP8誘導(dǎo)了廣泛的生物學(xué)過程,已被證實(shí)是一種控制焦亡、細(xì)胞凋亡和壞死性凋亡的分子開關(guān)[23]。

與單獨(dú)使用一種機(jī)器學(xué)習(xí)算法相比,整合多種機(jī)器學(xué)習(xí)算法構(gòu)建的疾病預(yù)測標(biāo)志物具有更優(yōu)異的性能[24]。Liu等[25]發(fā)現(xiàn),整合LASSO和SVM算法構(gòu)建的預(yù)后標(biāo)志物組合表現(xiàn)出對乳腺癌生存率預(yù)測更強(qiáng)的效能。本研究基于上述4個與CRC相關(guān)的PANoptosis關(guān)鍵基因構(gòu)建了列線圖模型,發(fā)現(xiàn)該列線圖模型在預(yù)測CRC發(fā)生的臨床效能上表現(xiàn)良好。進(jìn)一步采用訓(xùn)練集和驗證集評估了列線圖模型診斷CRC的ROC曲線,發(fā)現(xiàn)AUC值均大于0.9,且其預(yù)測效能優(yōu)于各單基因,說明基于這4個PANoptosis關(guān)鍵基因的預(yù)測模型在CRC診斷方面具備潛在的優(yōu)勢。然而,本研究仍存在幾個局限性。首先,本研究僅納入了GEO數(shù)據(jù)庫中2個數(shù)據(jù)集樣本,未與其他數(shù)據(jù)庫進(jìn)行聯(lián)合比較,無法排除轉(zhuǎn)錄本可能存在的選擇偏倚。其次,未對4個PANoptosis關(guān)鍵基因在CRC細(xì)胞中的功能和潛在機(jī)制進(jìn)行評估,仍需要額外的體內(nèi)和體外實(shí)驗來進(jìn)行功能和臨床驗證。此外,在臨床樣本中對所篩選的基因和構(gòu)建的模型進(jìn)行驗證是十分必要的。

綜上所述,通過綜合運(yùn)用WGCNA和LASSO、SVM-RFE和RandomForest算法鑒定了4個PANoptosis關(guān)鍵基因,開發(fā)并驗證了基于這4個基因的列線圖預(yù)測模型。結(jié)果顯示,構(gòu)建的列線圖模型可能成為診斷CRC的有價值工具。本研究為后續(xù)進(jìn)一步深入探討CRC的發(fā)病機(jī)制,開發(fā)新型診斷生物標(biāo)志物提供了新思路。