基于生物信息學篩選炎癥相關結直腸癌基因及其特征分析

劉欣然 康悅 于守江 于永生 劉玉偉 王洪偉

結直腸癌作為最常見的消化道惡性腫瘤之一，嚴重影響著公眾健康[1]。結直腸癌的發病機制一直是腫瘤研究領域的重點問題，而其中炎癌轉化作為癌變潛在的原因越來越多地受到研究者關注[2-3]。炎癥和腫瘤之間存在千絲萬縷的聯系，既往研究已經證實炎癥性腸病患者患結直腸癌的風險明顯增加[4-5]，而口服非甾體抗炎藥能夠降低腸道腫瘤的發生率及炎癥性腸病癌變的風險，一定程度上驗證了炎癌轉化這一機制的存在[6-7]。但是由于炎癥性腸病發生癌變的過程較為隱匿且時間較長，也對這一過程的研究增加了一定難度[8]。

近年來，隨著大數據和多組學技術的發展，生物信息學在疾病發生與發展機制探究和疾病診斷和分型的應用中越來越廣泛[9]。包括表達譜芯片和轉錄組測序在內的技術方法能夠通過基因差異表達的篩選，確定不同病理狀態下的轉錄差異，找出全新生物標志物，從而更好地探究和理解疾病的發生發展機制，并進一步用于臨床的診斷及治療。與此同時，GEO，TCGA等數據庫將不同中心的數據進行匯總儲存，并且作為開源數據為研究者提供下載服務，成為當今生物科學，尤其是為腫瘤領域的研究開啟全新的格局。

本研究結合GEO和TCGA兩個公共數據庫，從生物信息學角度著手獲得炎癥性腸病相關的差異基因，然后結合結直腸癌表達譜和生存數據，篩選炎癥性腸癌相關的候選基因，為后續炎癌轉化機制和炎癥性腸癌預防及治療提供研究基礎和出發點。

資料與方法

一、在GEO數據庫中篩選炎癥性腸病相關基因

在GEO中篩選炎癥性腸病表達譜數據，經過全面篩選獲取RNA-seq數據集GSE95473和GSE107597。前者數據集中包含25例潰瘍型結腸炎（活動期）病灶標本，20例克隆恩病（活動期）病灶標本，還有29例健康對照組正常腸黏膜標本。后者則包括24例潰瘍型結腸炎的炎癥病灶和配對的非炎性的黏膜。在確定了候選數據集后，則進行下一步差異基因篩選，利用BioJupies分別對兩個數據集中的病灶和腸黏膜的數據進行處理。通過使用limma R包比較對照組和實驗組之間的基因表達水平來產生基因表達特征（表1）。

表1 來自GEO數據庫的數據集的相關信息

通過基于R語言limma包中的DESeq2法進行差異基因表達分析，對炎性病灶和正常黏膜樣品之間的差異表達基因（differentially expressed genes，DEGs）進行分析篩選。DEGs的閾值設定為P＜0.01和 |logFoldChange|＞ 2。

二、功能富集分析

通過基因本體論（gene ontology，GO）富集分析可以利用這些信息來識別生物過程，分子功能和細胞成分。通過從基因表達特征中分別提取具有最高和最低值的500個基因來產生上調和下調的基因組。使用Benjamini-Hochberg校正后，P＜0.1作為截止值。同時，生物學通路是細胞學行為中起關鍵作用的生物化合物之間相互作用序列。包括KEGG在內的數據庫包含了這些基因和通路的大量關聯信息，使用Benjamini-Hochberg校正后，P＜0.1作為截止值。

三、通過TCGA進一步篩選結直腸癌預后相關基因及數據庫內驗證

在TCGA數據庫中查找與結直腸癌預后相關的基因，篩選出的差異表達基因譜進行取交集，篩選出結直腸癌預后相關的基因。與此同時，利用GEPIA的工具（http://gepia.cancer-pku.cn/index.html）分析來自TCGA項目的RNA測序表達數據，用于本研究中的批量TCGA數據處理。差異分析的方法采用單向ANOVA，使用疾病狀態（腫瘤或正常）作為計算差異表達的變量。

四、統計學分析

采用SPSS 20.0統計軟件對數據進行處理。計量資料采用均值±標準差表示，計數資料采用頻率和百分比表示。采用獨立樣本t檢驗和單向ANOVA評估連續變量的統計學顯著性，通過Tukey檢驗進行多重比較。通過卡方檢驗或Fisher精確檢驗來分析等級資料之間的差異。采用Kaplan-Meier方法繪制存活曲線，并使用log-rank檢驗方法比較不同分組的結直腸癌患者總生存期。P＜0.05被認為差異具有統計學意義。

結 果

一、GEO數據庫中炎癥性腸病相關基因篩選結果

我們在GEO中篩選炎癥性腸病表達譜數據，經過全面篩選獲取RNA-seq數據集GSE95473和GSE107597。前者數據集中包含25例潰瘍型結腸炎（活動期）病灶標本，20例克隆恩病（活動期）病灶標本，還有29例健康對照組正常腸黏膜標本。后者則包括24例潰瘍型結腸炎的炎癥病灶和配對的非炎性的黏膜。

采用將P＜0.01和|logFoldChange|＞2作為閾值，其中GSE95473中通過比較炎癥性腸病病灶和正常對照組腸黏膜進行差異表法分析，獲得114個DEGs，同樣地通過比較GSE107597中炎癥病灶和非炎癥黏膜，篩選到了709個DEGs。兩個數據集篩選到的DEGs按照|logFoldChange|由大到小的順序進行排列，前30的基因展示在表2和表3中。

二、炎癥性腸病相關差異表達基因的功能富集分析

將上述兩個數據集篩選出的基因取交集，獲得100個DEGs。為了進一步探索其腫瘤生物學行為，我們對其進行了功能富集分析。100個DEGs按照GO細胞組分、分子功能和生物學過程進行注釋。并富集到包括KEGG在內信號通路中。結果證實，在細胞組分上與顆粒腔、囊泡腔以及各種顆粒有關。在這分子功能上，與趨化因子受體結合、趨化因子的激活以及各種金屬離子的結合有關。在生物過程中，則與炎性反應、粒細胞調控、中性粒細胞介導的免疫反應和趨化因子介導的一系列信號通路有關。在信號通路上與IL-17信號通路，趨化因子-趨化因子受體相互作用、趨化因子信號通路等關系密切，見圖1。

三、通過TCGA篩選出提示結直腸癌預后不良的基因

為了建立炎癥性腸病和結直腸癌之間的關聯，我們將從上述兩個炎癥性腸病數據集篩選到的100個DEGs在TCGA腫瘤數據庫中進行驗證。基于上述差異基因在TCGA數據庫中對結腸腫瘤及正常黏膜、直腸腫瘤及正常黏膜表達譜數據進行主成分分析。從結果可以看出，數據結果能夠被明顯聚成兩類，腫瘤跟正常組織能夠明顯區分，進而提示這100個DEGs跟結直腸癌的表達特征存在內部關聯，見圖2。

與此同時，為了進一步探究上述基因跟腫瘤預后的關系，我們篩選結直腸癌數據（COAD＋READ）中總生存率LogrankP＜0.05的差異表達基因，共獲得10個。分別為AQP8，ALDOB，SPINK4，REG4，IL1B，C2CD4A，TNIP3，CXCL8，NOS2，CXCL3，見圖 3。

四、在TCGA數據庫中探索候選基因在結直腸癌中的表達模式

為了進一步研究上述基因在結直腸癌中的表達模式，我們在TCGA數據庫中進一步探究其表達。結果顯示AQP8在結腸癌數據集中表達下降，TNIP3在結腸和直腸兩個數據集中腫瘤灶與正常黏膜之間的表達差異無統計學意義，而其他8個基因和正常黏膜相比腫瘤中的表達明顯升高，差異有統計學意義，見圖4。

討 論

隨著組學技術的應用和發展，有學者根據結直腸癌的病理特征和特征性表達譜將其分為杯狀細胞型、腸上皮細胞型、干細胞型、炎癥型和轉化放大型[10]。同時發現其中的炎癥型腸癌患者對放療相對敏感，并且可以通過口服阿司匹林等抗炎治療能夠達到預防腸癌發生的目的[11]。炎癥性腸病作為結直腸癌發生的高危因素，慢性炎癥腸病患結直腸癌的風險明顯高于普通人群[12]。與此同時，也有文獻證實對于結直腸癌的影響不僅限于炎癥性腸病，結直腸癌的發生發展與某些隱匿性炎癥有關[13]。炎癥性腸病相關結直腸癌的研究對于炎癌轉化和炎癥相關結直腸癌的發生與發展具有重要的價值[14]。

表2 GSE95473差異表達基因（|logFoldChange|由大到小的順序排名前30的基因）

本研究首先在GEO數據庫中找到兩個炎癥性腸病相關的RNA測序數據，通過常規的差異表達分析流程，獲得了100個共有的DEGs。這些基因能夠反映出炎癥性腸病的表達特征，在此基礎上我們提出假設，即這些基因中包含腸癌發生與發展過程相關的基因。為了驗證這一想法，我們對TCGA中的結直腸癌的腫瘤病灶和正常黏膜中這100個基因進行主成分分析，從結果可以明顯看出這些基因能夠很好區分腫瘤和正常組織。此外，如果一個基因的表達能夠影響腫瘤的預后，我們可以認為該基因在腫瘤的病理生理過程中發揮一定的作用，基于此點，本研究中利用TCGA數據庫中生存相關的數據信息，在結直腸癌中篩選有預后意義的基因。同時，為了進一步了解這些基因在腫瘤和正常黏膜中表達情況，我們再一次在TCGA的數據集中進行驗證，發現10個基因中有8個在結腸癌和直腸癌中的表達均顯著升高，而另外兩個基因AQP8和TNIP3則沒有這個趨勢，這可能跟這兩種基因本身的表達量基線較低有關。這也從另一個角度說明通過預后來篩選腫瘤差異基因的科學性和可行性。

表3 GSE107597差異表達基因（|logFoldChange|由大到小的順序排名前30的基因）

通過上述流程，我們最終獲得了8個候選基因。在這些基因中，已有一部分已在多個研究中證實其介導炎癥反應，并同時在結直腸癌發生發展甚至是遠處轉移過程中發揮重要作用，如IL1B，CXCL8，NOS2，CXCL3[15-18]。既往研究已經證明，IL-1β作為一個促炎因子，炎癥性腸病活動期往往伴隨高水平的IL-1β，與此同時IL-1β能夠促進腫瘤侵襲，并且在多種消化道腫瘤的發生與進展中發揮作用[15，19-21]。同時在體外實驗中證實IL-1β可以調節VEGF-α和COX-2的表達，參與腫瘤血管的形成[22-23]。已有文獻證實CXCL8與結直腸癌關系密切，CXCL8能夠誘導結腸癌上皮間質轉化、同時抵抗失巢凋亡，促進腫瘤的侵襲、遷移，還可介導腫瘤細胞-免疫細胞之間的相互作用，負向調節機體免疫功能[16，24-25]。而誘導型一氧化氮合酶（nitric oxide synthase 2，NOS2）在結直腸癌中的表達有助于判斷結直腸癌的惡性程度，并且抑制NOS2在結直腸癌中表達可以抑制腫瘤血管生成[26-28]。有另研究表明，CXCL3的表達水平在結直腸腫瘤組織中顯著提高，并且具有潛在的診斷價值，CXCL3可能通過調控DNA修復，細胞凋亡及TP53的表達水平影響患者預后[18，29]。本研究還發現了ALDOB，SPINK4，REG4，C2CD4A等四個候選基因，這些基因既往關于結直腸癌的炎癌轉化相關研究較少，而本研究發現這些基因在結直腸腫瘤中高表達，并且這些基因的高表達往往提示患者預后較好。而這些相關轉化機制和人群隊列研究有待于后續進一步的探索和驗證。

圖1 GO富集和KEGG通路富集分析結果。該圖包含交互式條形圖，顯示使用Enrichr生成的GO富集分析和KEGG信號通路的結果。條形圖的長度表示每個條目的富集分數

圖2 主成分析結果。基于兩個隊列重合的100個差異基因在TCGA數據庫中對結腸腫瘤及正常黏膜、直腸腫瘤及正常黏膜表達譜數據進行主成分分析。2A：是降維過程中每一個主成分中的比例。2B和2C：顯示了結直腸癌病灶及正常黏膜測序數據的散點圖（二維和三維）。每個點代表RNA-seq樣品。具有相似基因表達譜的樣品在三維空間中更接近

圖3 基于TCGA數據集候選基因在結直腸癌中的預后意義。根據TCGA的數據，候選基因的表達情況對結直腸癌的總生存期的影響，采用對數秩檢驗進行計算

關于尚存在以下局限有待于后續研究，首先，本課題篩選到了8個候選基因，尤其有4個既往研究報道較少的基因，有進一步探索的價值。同時本研究假設的前提是由炎癥性腸病發展而來的結直腸癌和散發性結直腸癌基因特征相近。因此在后續研究中我們會納入炎癥性腸病發展而來的結直腸癌的相關意義基因。其次，從當前篩選的方法來看，只能夠說明這些基因跟炎癥性腸病和結直腸癌的表達及預后有關，而具體這些基因是如何通過調節炎癥反應，誘導腫瘤的發生與發展，最終影響患者預后的，仍有待于進一步機制探究和模型驗證。第三，本研究應用公共數據來篩選和探究結直腸癌的表達模式，缺乏大樣本的隊列的驗證和臨床因素對于這些標志物表達水平相關性的研究；只是研究了炎癥性腸病的差異表達基因，以及候選基因的表達模式和預后意義。此外，本研究發現一個有趣的結果，篩選出的生存相關的差異表達基因，高水平表達提示預后良好，可能由于這些基因在腫瘤癌變之后，發揮炎癥反應的作用，抑制腫瘤細胞的惡性生物學行為，相對好的預后，這也是本研究發現的一個重要的結果，其生物學作用仍有待于進一步探索。

圖4 基于TCGA數據集候選基因在結直腸癌中的表達模式。結直腸癌原發灶和正常黏膜中的候選基因的表達情況，腫瘤灶用紅色表示，正常組織用綠色表示，*表示P＜0.05