線粒體功能相關(guān)基因作為結(jié)直腸癌診斷與預(yù)后標(biāo)志物的應(yīng)用研究

周凱翔，袁 晴，王珍妮，郭姍姍，蘇麗萍，劉 洋，郭 旭，邢金良，黃啟超

(1空軍軍醫(yī)大學(xué)基礎(chǔ)醫(yī)學(xué)院生理與病理生理學(xué)教研室，陜西 西安 710032； 2西北工業(yè)大學(xué)醫(yī)學(xué)研究院，陜西 西安 710072)

結(jié)直腸癌(colorectal cancer，CRC)是世界范圍內(nèi)發(fā)病率和死亡率最高的常見(jiàn)癌癥之一。最新數(shù)據(jù)顯示，CRC的發(fā)病率和死亡率在全球所有癌癥患者中分別位列第三和第二位[1]。與大多數(shù)癌癥一樣，CRC的早期診斷有助于良好的預(yù)后，已有流行病學(xué)調(diào)查結(jié)果表明早期發(fā)現(xiàn)疾病可使五年生存率提高到90%[2-3]。然而由于缺乏敏感有效的診斷方法，CRC通常到晚期才被診斷出來(lái)，極大地影響了CRC患者的臨床治療和生存質(zhì)量。盡管近年來(lái)在CRC診斷和治療方面有了顯著的改善，但被診斷為轉(zhuǎn)移性CRC患者的5年生存率仍然很低，約為12年[4]。因此，迫切需要了解CRC發(fā)生發(fā)展的分子機(jī)制，并尋找新的生物標(biāo)志物用于CRC的早期診斷和預(yù)后評(píng)估。

既往研究證實(shí)CRC是一種依賴于眾多致癌基因和抑癌基因改變的遺傳性疾病[5]。大量研究致力于鑒定CRC相關(guān)基因及其編碼蛋白在包括細(xì)胞增殖、分化、凋亡和轉(zhuǎn)移在內(nèi)的大量生理和病理過(guò)程中的重要作用。然而，CRC的確切分子機(jī)制仍有待深入研究。線粒體作為真核細(xì)胞內(nèi)重要的半自主細(xì)胞器，含有自身的線粒體DNA并編碼呼吸鏈的13種必需蛋白質(zhì)[6-7]。已有研究報(bào)道CRC細(xì)胞以線粒體氧化磷酸化作為其主要能量來(lái)源。此外，在CRC發(fā)生過(guò)程中，線粒體在調(diào)節(jié)細(xì)胞增殖和凋亡中起核心作用。除了由線粒體基因組編碼的37個(gè)基因外，參與調(diào)節(jié)線粒體結(jié)構(gòu)和功能的核基因超過(guò)1 600個(gè)[8]。因此，核基因組和線粒體基因組中許多基因的突變會(huì)導(dǎo)致線粒體功能障礙，從而導(dǎo)致包括腫瘤在內(nèi)的許多嚴(yán)重疾病。然而，線粒體相關(guān)基因在CRC發(fā)生發(fā)展及預(yù)后預(yù)測(cè)中的具體作用尚不明確。因此，闡明線粒體相關(guān)基因的異常表達(dá)機(jī)制，并探索這些基因作為潛在生物標(biāo)志物在CRC患者的有效診斷和預(yù)后評(píng)估中的具體作用，對(duì)揭示CRC新的診斷治療靶點(diǎn)和預(yù)后策略是十分重要的。

本研究基于腫瘤基因圖譜(The Cancer Genome Atlas，TCGA)數(shù)據(jù)庫(kù)中CRC患者的轉(zhuǎn)錄組數(shù)據(jù)集，分析了線粒體功能相關(guān)基因的表達(dá)差異及其潛在的生物學(xué)功能，并將線粒體功能相關(guān)基因納入診斷和預(yù)后評(píng)估的模型，進(jìn)而建立了聯(lián)合指標(biāo)的CRC診斷和預(yù)后預(yù)測(cè)模型，為CRC患者的診斷和預(yù)后評(píng)估提供臨床新思路。

1 材料與方法

1.1 材料

本文數(shù)據(jù)來(lái)源于TCGA數(shù)據(jù)庫(kù)，分別下載了結(jié)腸癌和直腸癌的RNA-seq counts表達(dá)譜數(shù)據(jù)及對(duì)應(yīng)的臨床信息。其中結(jié)腸癌包括478例腫瘤樣本和41例正常樣本，直腸癌包括166例腫瘤樣本和10例正常樣本。經(jīng)數(shù)據(jù)匯總后，本文分析的CRC數(shù)據(jù)包括644例腫瘤樣本和51例正常樣本。

線粒體功能相關(guān)基因信息來(lái)源于MitoMiner數(shù)據(jù)庫(kù)[9]，綜合線粒體蛋白指數(shù)數(shù)據(jù)庫(kù)收錄了1 626個(gè)線粒體功能相關(guān)基因的信息，其中包含了13個(gè)線粒體基因組編碼基因。本文將非線粒體編碼的1 613個(gè)線粒體相關(guān)基因及國(guó)家生物技術(shù)中心信息檢索數(shù)據(jù)庫(kù)線粒體參考基因組中的37個(gè)線粒體基因統(tǒng)稱為線粒體功能相關(guān)基因。對(duì)數(shù)據(jù)進(jìn)行質(zhì)控的結(jié)果表明，本研究所用數(shù)據(jù)集中的基因?qū)€粒體功能相關(guān)基因的覆蓋率為100%。

1.2 方法

1.2.1 表達(dá)譜數(shù)據(jù)的差異基因篩選 對(duì)CRC的RNA-seq counts數(shù)據(jù)進(jìn)行差異表達(dá)基因篩選，方法包括t檢驗(yàn)和倍數(shù)變化(fold change，F(xiàn)C)，之后采用多重檢驗(yàn)校正的方法Benjamini-Hochberg對(duì)P值進(jìn)行校正。本文通過(guò)R軟件的DESeq2程序包對(duì)RNA-seq counts數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理后，再分析差異基因，其中符合P.adj<0.05且|log2FC|≥1的基因在本文被定義為差異表達(dá)基因。

1.2.2 對(duì)差異表達(dá)線粒體功能相關(guān)基因進(jìn)行注釋 為更進(jìn)一步分析差異表達(dá)線粒體功能相關(guān)基因的生物學(xué)通路，本文通過(guò)R軟件的org.Hs.eg.db、clusterProfiler、DOSE、stringr等程序包進(jìn)行基因本體論(gene ontology，GO)功能注釋[10]及京都基因和基因組百科全書(shū)(Kyoto Encyclopedia of Genes and Genomes，KEGG)功能富集分析[11]，其中GO注釋為線粒體功能相關(guān)基因參與的生物學(xué)過(guò)程。采用多重檢驗(yàn)校正的方法Benjamini-Hochberg對(duì)P值進(jìn)行校正，設(shè)置滿足P.adj<0.01的GO術(shù)語(yǔ)和KEGG通路具有顯著性。

1.2.3 隨機(jī)森林構(gòu)建診斷模型 將差異的線粒體功能相關(guān)基因用作構(gòu)建診斷模型的不同特征，隨機(jī)森林的實(shí)現(xiàn)采用R語(yǔ)言caret程序包篩選與診斷相關(guān)的最優(yōu)線粒體相關(guān)特征基因。十折交叉驗(yàn)證的方法，將數(shù)據(jù)按照9∶1的比例隨機(jī)分成訓(xùn)練集和驗(yàn)證集，最后得到最優(yōu)的特征基因。

1.2.4 單因素Cox回歸分析 首先對(duì)臨床數(shù)據(jù)中的年齡、病理分期變量進(jìn)行單因素Cox回歸分析，再將單因素分析中P<0.05的變量進(jìn)行保留。其次，將差異表達(dá)的11個(gè)線粒體基因作為單因素分析中的變量，再將其中P<0.05的變量進(jìn)行保留。選取顯著的獨(dú)立變量納入模型構(gòu)建并使用受試者工作特征(receiver operating characteristic，ROC)曲線下面積(area under curve，AUC)進(jìn)行效能評(píng)價(jià)。

1.2.5 生存分析 將顯著的獨(dú)立變量納入模型，采用R語(yǔ)言的predict函數(shù)對(duì)變量進(jìn)行風(fēng)險(xiǎn)打分，高于中位數(shù)為高風(fēng)險(xiǎn)，低于中位數(shù)為低風(fēng)險(xiǎn)。最后用R語(yǔ)言的Survival和Survminer程序包繪制生存曲線。

2 結(jié)果

2.1 CRC差異表達(dá)基因篩選

在包含CRC組644例(結(jié)腸患者478例，直腸患者166例)和正常組51例的基因表達(dá)譜數(shù)據(jù)中，采用T檢驗(yàn)和FC的差異基因分析方法設(shè)置P.adj<0.05且|log2FC|≥1的閾值，共篩選出15 115個(gè)CRC相關(guān)的差異表達(dá)基因，上調(diào)基因10 762個(gè)，下調(diào)基因4 353個(gè)。其中，差異表達(dá)的線粒體功能相關(guān)基因有257個(gè)，上調(diào)136個(gè)，下調(diào)121個(gè)(圖1)。圖1展示了1 650個(gè)線粒體功能相關(guān)基因在CRC組和正常組中基因表達(dá)值的情況。差異表達(dá)的11個(gè)線粒體基因的詳細(xì)信息見(jiàn)表1，其中MT-ND6基因?yàn)橄抡{(diào)基因，剩余10個(gè)基因?yàn)樯险{(diào)基因。

圖1 線粒體功能相關(guān)基因火山圖

表1 11個(gè)差異表達(dá)線粒體基因

2.2 CRC差異表達(dá)線粒體功能相關(guān)基因的功能分析

對(duì)257個(gè)差異表達(dá)線粒體功能相關(guān)基因進(jìn)行GO功能注釋和KEGG功能富集分析，結(jié)果分別富集到2 910條GO生物學(xué)通路和205條KEGG通路。符合閾值P.adj<0.01的生物學(xué)術(shù)語(yǔ)154條(圖2展示了最顯著的15條)，它們主要與小分子分解代謝過(guò)程、含嘌呤化合物代謝過(guò)程和有機(jī)酸分解代謝過(guò)程等有關(guān)；符合閾值P.adj<0.01的KEGG通路13條(圖3)，它們主要與脂肪酸代謝、精氨酸和脯氨酸代謝及脂肪酸降解等有關(guān)。

縱坐標(biāo)表示不同的GO術(shù)語(yǔ)，橫坐標(biāo)表示注釋到某條GO術(shù)語(yǔ)的基因比例。圖中不同P.adj值表示為不同色階，紅色表示P.adj值較小，藍(lán)色表示P.adj值較大。圖2 差異表達(dá)線粒體功能相關(guān)基因的GO功能注釋

縱坐標(biāo)表示不同的KEGG通路，橫坐標(biāo)表示注釋到某條KEGG通路的基因比例。圖中不同P.adj值表示為不同色階，紅色表示P.adj值較小，藍(lán)色表示P.adj值較大。圖3 差異表達(dá)線粒體功能相關(guān)基因的KEGG功能富集

2.3 診斷特征篩選

十折交叉驗(yàn)證將數(shù)據(jù)按照9∶1的比例隨機(jī)分成訓(xùn)練集和驗(yàn)證集，得到8個(gè)最優(yōu)特征基因(圖4A)。隨后對(duì)其進(jìn)行評(píng)估，繪制出ROC曲線，其AUC值為0.919(圖4B)。

A：隨機(jī)森林篩選后的特征基因；B：特征基因的ROC曲線。圖4 診斷特征篩選結(jié)果

2.4 生存率的獨(dú)立預(yù)測(cè)因素

644例腫瘤樣本中納入生存分析的CRC患者600例，選取P<0.05的因素作為模型構(gòu)建，通過(guò)單因素Cox回歸分析，發(fā)現(xiàn)年齡 ≥ 67歲、Stage Ⅲ & Ⅳ和MT-TL1基因是生存曲線的獨(dú)立危險(xiǎn)因素(表2)。

表2 CRC患者預(yù)后單因素分析

2.5 預(yù)后評(píng)估模型的構(gòu)建及評(píng)估

通過(guò)對(duì)單因素分析顯著的變量，年齡、Stage分期、MT-TL1基因表達(dá)值模型的建立，采用R語(yǔ)言的predict函數(shù)對(duì)多因素進(jìn)行風(fēng)險(xiǎn)打分，高于中位數(shù)為高風(fēng)險(xiǎn)組，低于中位數(shù)為低風(fēng)險(xiǎn)組。通過(guò)風(fēng)險(xiǎn)得分來(lái)繪制3年和5年的ROC曲線(圖5)，結(jié)果表明3年生存的AUC值為0.749，5年生存的AUC值為0.721。

圖5 3年生存和5年生存的ROC曲線

2.6 生存曲線繪制

經(jīng)過(guò)模型構(gòu)建的評(píng)估，并結(jié)合臨床信息將600例CRC劃分高風(fēng)險(xiǎn)組(221例)和低風(fēng)險(xiǎn)組(379例)，繪制總生存曲線(P<0.000 1，圖6)。此結(jié)果表明高風(fēng)險(xiǎn)組生存較差，低風(fēng)險(xiǎn)組生存較好。

圖6 總生存曲線

3 討論

在本研究中，我們綜合分析了來(lái)自TCGA的全部腸癌相關(guān)轉(zhuǎn)錄組測(cè)序數(shù)據(jù)。CRC組織與正常組織共檢測(cè)到257個(gè)線粒體相關(guān)基因表達(dá)，其中上調(diào)基因表達(dá)136個(gè)，下調(diào)基因表達(dá)121個(gè)。功能富集分析表明，線粒體相關(guān)差異基因在小分子分解代謝、含嘌呤化合物代謝和脂肪酸代謝等生物學(xué)過(guò)程中富集。之前研究報(bào)道，將羧酸酯酶1(carboxylesterase 1，CES1)鑒定為一種必需的NF-κB調(diào)節(jié)的脂肪酶，它將肥胖相關(guān)炎癥與脂肪代謝和侵襲性CRC中的能量壓力適應(yīng)聯(lián)系起來(lái)。CES1通過(guò)促進(jìn)脂肪酸氧化的細(xì)胞自主機(jī)制促進(jìn)CRC細(xì)胞存活并防止三?；视偷亩拘苑e聚。研究人員發(fā)現(xiàn)CES1表達(dá)升高與超重CRC患者的較差結(jié)果相關(guān)[12-13]。另外，膽汁酸是由肝臟中的膽固醇合成的生理產(chǎn)物，具有產(chǎn)生膽汁流動(dòng)以及促進(jìn)腸道對(duì)脂質(zhì)、營(yíng)養(yǎng)素和維生素的吸收和運(yùn)輸?shù)淖饔?。在CRC的癌變和發(fā)展過(guò)程中，膽汁酸似乎起著至關(guān)重要的作用。研究人員揭示了石膽酸和脫氧膽酸在N-甲基-N′-硝基-N-亞硝基胍治療后促進(jìn)大鼠CRC腫廇生長(zhǎng)。BERNSTEIN等[14]進(jìn)一步證實(shí)了脫氧膽酸在小鼠實(shí)驗(yàn)中的作用；他們用含2 mL/L脫氧膽酸的飲食喂養(yǎng)18只小鼠，8～10個(gè)月后，其中10只患上CRC。近年來(lái)，越來(lái)越多的研究發(fā)現(xiàn)次級(jí)膽汁酸石膽酸和脫氧膽酸是參與CRC發(fā)生和發(fā)展的主要因素。我們的研究結(jié)果表明，這些線粒體相關(guān)差異基因可能參與了CRC的不同代謝過(guò)程，從而導(dǎo)致腫瘤的發(fā)生和進(jìn)展。

近年來(lái)，許多研究表明機(jī)器學(xué)習(xí)技術(shù)在癌癥診斷領(lǐng)域具有廣泛前景。例如，PU等[15-16]利用支持向量機(jī)方法確定了一個(gè)基于多個(gè)超甲基化CpG位點(diǎn)的診斷模型，準(zhǔn)確率為0.82%。本研究以695個(gè)樣本的257個(gè)差異線粒體相關(guān)基因的表達(dá)譜為原始數(shù)據(jù)，通過(guò)隨機(jī)森林構(gòu)建診斷模型。經(jīng)過(guò)十折交叉驗(yàn)證的方法，結(jié)果表明8個(gè)線粒體相關(guān)基因集合(ABCA8、SLC25A34、MT-ND6、P2RY1、ABCG2、MT-TV、MT-TL1和ADHFE1)對(duì)CRC腫廇生長(zhǎng)有較好的診斷性能，其AUC值為0.919。有趣的是，研究人員最近從一小部分患者中報(bào)道了結(jié)腸腺瘤(帶蒂息肉)中ABCG2相對(duì)于相鄰的正常結(jié)腸黏膜的下調(diào)。這被認(rèn)為有利于致癌物質(zhì)在腸道中的積累，從而促進(jìn)CRC的形成。

許多研究已經(jīng)基于腫瘤相關(guān)差異基因構(gòu)建了一些預(yù)后模型來(lái)預(yù)測(cè)腫瘤患者的預(yù)后。例如，QIU等[17-18]基于生物信息學(xué)分析發(fā)現(xiàn)了口腔舌鱗狀細(xì)胞癌患者的16個(gè)基因預(yù)后特征，其預(yù)測(cè)5年生存率的AUC為0.872。LIU等[19-20]通過(guò)單因素和多因素Cox比例風(fēng)險(xiǎn)回歸模型研究胃癌預(yù)后相關(guān)基因，他們發(fā)現(xiàn)一組由9個(gè)基因組成的預(yù)后基因標(biāo)記，其與胃癌患者的預(yù)后有顯著相關(guān)性。在本研究中，為了評(píng)估CRC患者線粒體相關(guān)差異基因與臨床結(jié)局之間的關(guān)系，我們首先進(jìn)行了單因素Cox比例風(fēng)險(xiǎn)回歸，結(jié)果發(fā)現(xiàn)了年齡、腫瘤病理分期和MT-TL1基因與生存相關(guān)。最后，我們通過(guò)風(fēng)險(xiǎn)預(yù)測(cè)模型建立了基于上述三個(gè)因素的生存預(yù)測(cè)模型。經(jīng)風(fēng)險(xiǎn)評(píng)分分層后，預(yù)后特征的Kaplan-Meier曲線顯示，CRC患者的臨床結(jié)局存在顯著差異(P<0.000 1)。ROC分析獲得的AUC 3年和5年生存預(yù)測(cè)分別為0.749和0.721。這表明本CRC預(yù)后模型的特異性和敏感性相對(duì)較高。

本研究具有許多重要優(yōu)勢(shì)：第一，本研究納入了TCGA數(shù)據(jù)集中全部CRC相關(guān)的測(cè)序數(shù)據(jù)，使用的數(shù)據(jù)非常全面；第二，應(yīng)用基于隨機(jī)森林的機(jī)器學(xué)習(xí)方法和Cox回歸模型來(lái)構(gòu)建線粒體相關(guān)差異基因在CRC患者中的診斷及預(yù)后模型，應(yīng)用ROC曲線較為穩(wěn)健地評(píng)估效能。然而我們的研究仍然存在一些局限性：第一，本研究的研究樣本來(lái)自白種人、亞洲人、黑人或非裔美國(guó)人等不同人群，不同民族間基因表達(dá)譜可能存在一定差異。為了彌補(bǔ)這一缺陷，我們未來(lái)將繼續(xù)收集CRC患者臨床數(shù)據(jù)及腫瘤組織樣本，利用獨(dú)立的實(shí)驗(yàn)數(shù)據(jù)進(jìn)一步驗(yàn)證我們的預(yù)測(cè)模型。第二，線粒體相關(guān)差異基因的診斷效率和預(yù)后價(jià)值僅在TCGA數(shù)據(jù)集中進(jìn)行了分析和驗(yàn)證。因此，未來(lái)的研究中我們將進(jìn)一步在其他數(shù)據(jù)庫(kù)中驗(yàn)證本結(jié)果。第三，我們的研究方法主要基于生物信息學(xué)分析，需要進(jìn)一步的實(shí)驗(yàn)來(lái)驗(yàn)證基于腫瘤樣本和臨床數(shù)據(jù)的結(jié)果。此外，體內(nèi)和體外實(shí)驗(yàn)可以增強(qiáng)我們對(duì)CRC中線粒體相關(guān)差異基因功能作用的認(rèn)識(shí)。

綜上所述，本研究通過(guò)在大樣本中分析CRC的線粒體功能相關(guān)基因的表達(dá)，闡述了線粒體功能在CRC發(fā)生發(fā)展中的可能作用，結(jié)合臨床指標(biāo)構(gòu)建并評(píng)估了CRC診斷和預(yù)后模型，為CRC患者的預(yù)后判斷提供了有效的新方法。