基于生物信息學分析潰瘍性結腸炎的差異表達基因和相關微小RNA

吳勝男, 蔣環宇, 陳浩然, 王欣瑤, 吳佳輝, 王璐琦

(山西醫科大學, 1. 管理學院, 2. 基礎醫學院, 山西 太原, 030000)

潰瘍性結腸炎(UC)是一種慢性非特異性炎癥,其特點是結直腸黏膜連續性、彌漫性發炎。UC作為炎癥性腸病(IBD)的一個重要類型,主要表現包括腹瀉、血便、腹痛和里急后重等癥狀,對患者的心理、經濟狀況和生活質量造成了嚴重的影響[1-3]。UC在工業化國家更為常見,特別是北美和西歐國家[4]。在IBD發病率迅速上升的亞太地區國家中,中國的發病率最高, UC是最常見的類型[5]。

目前UC的發病機制尚不明確,臨床很難根治。UC的發病主要是由于腸道黏膜屏障功能的損傷、黏膜免疫失衡以及腸道菌群紊亂,同時參與發病的因素還有遺傳易感性和環境因素[6]。免疫反應在UC的發生過程中扮演了重要角色,喪失免疫耐受會增加各種免疫炎癥因子,刺激抗原特異性效應物的增殖,從而導致局部和全身炎癥[7]。腸道損傷時,固有免疫細胞中性粒細胞能夠迅速聚集到受損部位。中性粒細胞起到限制外源性微生物入侵的關鍵作用,并具備重要的免疫調節和殺菌功能,但過度的免疫反應卻可能加劇UC早期腸道炎癥和黏膜損傷[8]。相關研究[9]通過流式細胞術檢測發現,與健康個體相比,UC患者的外周血T細胞表現出細胞激活標志物(如人類白細胞抗原-DR和β1整合素)的增加,而幼稚T細胞特征性抗原L選擇素(CD62L)的表達下降,提示在UC患者中T細胞的活化程度增加。微小RNA(miRNA)可參與調節各種生理過程[10-11]。超過60%的人類蛋白質編碼基因含有至少1個保守的miRNA結合位點[12]。miRNA與多種疾病(包括UC)具有顯著相關性。

本研究對GEO數據庫中的健康個體與UC患者的芯片原始數據進行加權基因共表達網絡分析和差異分析,從而獲取差異表達基因,并對其進行富集分析,根據關鍵的基因預測與DEG相關的潛在miRNA。

1 材料與方法

1.1 數據獲取

通過NCBI的GEO數據庫(https://www.ncbi.nlm.nih.gov/geo/)以“ulcerative colitis” “ulcerative colitis clinic”為關鍵詞進行檢索,限制Entry type為“Series”、Top Organisms為“Homo sapiens”, 篩選之后獲得GSE66407、GSE65114基因芯片數據作為研究對象,所在平臺分別為GPL19833以及GPL16686。GSE66407中包含368個樣本,通過篩選疾病以及是否具有表型數據,最終納入的數據為對照組99個樣本, UC組為156個樣本; GSE65114中包含28個樣本,其中對照組為12個, UC組為16個,全部納入研究。

1.2 數據預處理及篩選差異表達基因

通過bioDBnet(https://biodbnet-abcc.ncifcrf.gov/)對GSE66407數據進行id轉化,用于后續差異分析; 通過R軟件中的dplyr包對GSE65114數據進行id轉化,最終得到二者相對應的基因表達矩陣文件。對于GSE65114的基因表達矩陣文件,采用R軟件中的limma包進行差異基因篩選。差異基因的篩選標準為|logFC|>1以及經過調整后P值<0.05。

1.3 基因集富集分析(GSEA)

本研究使用GSEA軟件(version 4.2.3)對GSE65114數據進行分組。為評估相關的途徑和分子機制,從MolecularSignatures Database(https://www.gsea-msigdb.org/gsea/msigdb)下載了c2.cp.kegg.v2023.1.Hs.symbols.gmt子集合。設置基因集的大小范圍,最小基因集為15, 最大基因集為500, 并進行1 000次重復抽樣。設定P值<0.05和FDR<0.25的閾值。通過基因表達譜數據和表型數據進行分析,得出了相關結果。

1.4 構建加權基因共表達網絡(WGCNA)

當基因在生物學上顯示高度相關時,使用WGCNA將其聚類到不同的模塊中。此外, WGCNA還能夠探索模塊與外部性狀之間的關聯性,進而有助于識別候選的生物標志物和潛在的治療靶點[13]。本研究的目的在于整合GSE66407的基因表達矩陣文件和臨床表型文件,并利用R軟件中提供的WGCNA包來構建WGCNA。本研究首先對數據進行聚類,以便檢測異常值,并確定適當的軟閾值來建立一個無尺度的網絡。然后,采用層次聚類樹和拓撲重疊矩陣的方法來鑒定和檢測基因模塊。為了確定各個模塊與表型特征的相關性,本研究計算了皮爾遜相關系數,并找出與表型特征相關性最強的基因模塊。

1.5 交集基因的功能富集分析

首先將差異基因與基因模塊中的基因導入venny 2.1.0中進行交集分析,以獲取交集基因。然后使用R軟件中的ClusterProfiler包對這些交集基因進行功能富集分析,包括基因本體論(GO)和京都基因與基因組百科全書(KEGG)通路分析。GO功能分析包括生物學過程(BP)、細胞組分(CC)以及分子功能(MF)。設定P值<0.05作為篩選標準,用于確定顯著富集的基因集條目。

1.6 蛋白互作網絡的構建

為了揭示UC相關蛋白質的作用及其相互關系,將交集基因導入在線數據庫STRING (https://cn.string-db.org/)中,以構建蛋白質相互作用(PPI)網絡。在構建過程中,設定confidence>0.4作為篩選標準。隨后,利用Cytoscape軟件對構建好的PPI網絡進行可視化操作。通過使用CytoNCA插件中Degree值的排序,選出排名前10位的基因作為樞紐基因(Hub基因)。

1.7 與Hub基因相關的miRNA

使用NetworkAnalyst 3.0(https://www.networkanalyst.ca/)在線平臺[14], 將得到的10個Hub基因映射到各自的miRNA, 用來識別基因調控網絡中miRNA-基因之間的相互作用。最后,通過Cytoscape軟件對miRNA-基因網絡圖進行可視化。

2 結 果

2.1 差異基因的篩選結果

基于GSE65114的基因表達數據,將對照組和UC組進行比較。為了篩選差異基因并進行可視化,使用R軟件中的ggplot2包。篩選標準為: |logFC|>1且P<0.05。最終,經過數據處理,得到了277個差異表達基因,其中有200個基因上調,77個基因下調。詳細的基因熱圖和火山圖見圖1。在熱圖中,紅色代表基因上調,綠色代表基因下調,而黑色則表示中間表達水平。基因表達越顯著,顏色越深。

2.2 GSEA分析結果

通過對GSE65114數據中所有基因的表達矩陣進行GSEA分析,結果顯示富集通路主要涉及神經活性配體受體相互作用通路、利什曼原蟲感染、朊病毒病害以及心電圖受體相互作用通路等,見圖2。

2.3 基因共表達網絡

本研究使用R軟件中的WGCNA包對GSE66407的基因表達數據進行分析,首先利用goodSamplesGenes方法檢測缺失值,并利用平均值計算基因的平均表達量,設置過濾標準為FPKM>0.5。通過WGCNA算法計算合適加權系數β(軟閾值),得到β為11時構建的網絡更符合無尺度特征(圖3A)。接著通過層次聚類法,使用動態剪切樹進行模塊的識別和合并,最終獲得16個共表達模塊(圖3B)。然后隨機選取1 000個基因做TOM熱圖(圖3C)發現了明顯的顏色聚集,用來描述基因之間的拓撲重疊矩陣,確定每個模塊獨立存在。

2.4 模塊與性狀關聯分析

本研究將樣本臨床特征中的年齡、炎癥、組織納入表型數據,將每個模塊與樣本的性狀特征進行關聯分析,發現brown模塊中的基因與UC樣本中發炎結腸的相關性最強(cor=0.69,P<0.01)(圖4A)。

對統計棕色模塊中基因的基因與性狀關系(GS)和基因與模塊關系(MM)進行分析,可以得出棕色模塊內基因與顯著性關系的散點圖(圖4C), 由散點圖中可發現,基因的GS與MM在棕色模型中基本上都是呈現線性排布,存在著相關性與一致性,說明該模塊與UC炎癥的發生相關性最大,該模塊中可能存在調控UC炎癥反應的關鍵基因,因此之后將對棕色模塊中的基因進一步分析。

2.5 交集基因的富集分析結果

利用GEO數據庫中的GSE66407和GSE65114基因表達數據集進行分析,將獲得的277個差異基因與基因模塊中的2 027個基因導入venny 2.1.0中取交集,獲得123個交集基因(圖5A), 其中上調基因114個,下調基因9個。利用R軟件中的ClusterProfiler包對所獲得的交集基因進行GO功能和KEGG通路富集分析,設定P<0.05為富集條目的篩選標準。

GO富集結果(圖5B)顯示,與BP相關的條目主要富集在白細胞遷移、正向調節對外部刺激的反應、白細胞趨化性、髓系白細胞遷移、單核細胞遷移、中性粒細胞趨化性等;與CC相關的條目主要富集在膠原蛋白、細胞外基質分泌顆粒腔、細胞質囊泡腔、等離子體的囊泡腔外側、細胞膜頂端部分等;與MF相關的條目主要富集在結構成分趨化因子活性、肝素結合、趨化因子受體結合、Toll樣受體結合、細胞因子活性等。KEGG富集分析結果(圖5C)顯示,前12條通路主要富集在細胞因子受體相互作用、磷脂酰肌醇-3激酶/蛋白激酶B(PI3K-AKT)信號通路、趨化因子信號通路、核轉錄因子kappa B(NF-κB)信號通路、Toll樣受體信號傳導通路、白細胞介素-17(IL-17)信號通路等相關過程。

2.6 PPI網絡及Hub基因的篩選結果

將交集基因導入String在線數據庫,從中獲取PPI網絡(見圖6A)。然后,利用Cytoscape軟件對所獲得的PPI網絡進行可視化操作(見圖6B), 使用CytoNCA插件計算節點的Degree值, Degree值越大,節點越大。選出排名前10位的基因作為Hub基因,分別是C-X-C趨化因子配體8(CXCL8)、Toll樣受體2(TLR2)、細胞間黏附分子1(ICAM1)、選擇素L(SELL)、趨化因子受體4(CXCR4)、細胞毒性T淋巴細胞相關抗原(CTLA4)、細胞分化抗原69(CD69)、雙糖鏈蛋白多糖(BGN)、C-X-C趨化因子配體13(CXCL13)、金屬蛋白酶抑制劑1(TIMP1)。將這10個關鍵基因建立Hub基因網絡(見圖6C)。

2.7 建立miRNAs-Hub基因調控網絡

miRNA在基因表達的調控中起著多種作用。基于NetworkAnalyst數據庫,使用Cytoscape構建miRNA-Hub基因調控網絡。Hub基因及其相應的調控miRNA見圖7。其中, miRNA (hsa-miR-335-5p)有5個靶基因(CXCL8、TLR2、ICAM1、CXCR4以及CD69); miRNA(hsa-miR-146a-5p)有4個靶基因(CXCL8、TLR2、ICAM1、CXCR4); miRNA(hsa-miR-92a-3p)有3個靶基因(ICAM1、CTLA4、CD69); miRNA(hsa-miR-155-5p)有3個靶基因(CXCL8、ICAM1、CTLA4); miRNA(hsa-miR-26b-5p)有3個靶基因(CTLA4、CXCL13、TIMP1); 另外,CXCL8、ICAM1是3個miRNA的共同靶標(hsa-miR-4426、hsa-miR-4462b和hsa-miR-4647);ICAM1、CD69是2個miRNA的共同靶標(hsa-miR-32-5p和hsa-miR-92b-3p);CXCL8、ICAM1是2個miRNA的共同靶標(hsa-miR-98-5p和hsa-miR-93-5p)。

3 討 論

UC是IBD中的一種重要類型,其主要的特征是結直腸黏膜連續受累和彌漫型炎癥。據認為, UC的發生是易感個體在環境因素作用下產生的先天性和適應性免疫功能障礙,導致的慢性、非受限制的炎癥過程。這種炎癥過程受遺傳因素影響,并以免疫介導的方式持續存在。UC患者的病理學特點主要表現為結腸黏膜連續彌漫性慢性活動性炎癥[15]。因此,探究發炎黏膜中的潛在生物標志物,對于獲取新的治療靶點尤為重要,同時對于后續臨床診斷及治療UC具有重要意義。

本研究利用GEO數據庫中的GSE66407和GSE65114基因表達數據集進行了分析。首先,通過應用WGCNA和差異分析得到的差異基因取交集,獲得了123個交集基因。其中, 117個基因表達上調,9個基因表達下調。然后對這些交集基因進行了功能富集分析,并構建了一個蛋白互作網絡。這個過程篩選出了10個關鍵基因,包括CXCL8、TLR2、ICAM1、SELL和CXCL4等,這些基因在UC的發生和發展中可能起到重要的作用。

GO和KEGG功能富集分析結果顯示,核心基因與白細胞遷移、白細胞介素、趨化因子活性以及膠原蛋白密切相關。涉及的信號通路包括PI3K-AKT信號通路、NF-κB信號通路、Toll樣受體信號傳導通路以及IL-17信號通路等,這些通路與炎癥相關信號的激活和促炎因子的釋放密切相關,提示其在UC的發病和加重過程中起著重要的作用。白細胞遷移參與IBD的進展,因此阻斷白細胞向腸道遷移是控制疾病和緩解癥狀的主要策略[16]。研究[17]表明, PI3K-AKT細胞信號通路的異常被認為與UC的發病有關,對UC誘導的細胞凋亡和炎癥以及長期UC誘導的結腸癌有很大影響。GO和KEGG富集通路分析結果表明,本研究中鑒定的關鍵基因可能通過上述途徑參與UC的進展。

炎性趨化因子CXCL8是最早和研究最深入的趨化因子之一,可由多種細胞類型釋放, CXCL8-CXCR1/2軸在炎癥反應及腫瘤生長中起多效作用,包括中性粒細胞和粒細胞募集和活化到炎癥部位,有助于消除病原體,在各種炎癥性疾病中起重要作用。與正常對照組相比,活動性UC結腸黏膜中CXCL8表達增加已得到廣泛證實, CXCL8表達上調與炎癥程度相關[18-19]。TLR2作為模式識別受體(PRR)家族中的一員,能夠感知病原相關分子模式(PAMP)和損傷相關分子模式。在機體免疫微環境中, TLR2信號通路對巨噬細胞、NK細胞、T細胞、肥大細胞等細胞的表型和功能產生影響,進而影響機體對病原體的清除能力和抗腫瘤能力[20]。ICAM-1是一種細胞表面糖蛋白和黏附受體,在許多細胞類型中充當黏附分子和信號受體,以引發炎癥反應,啟動炎癥和愈合的消退,并調節腫瘤細胞的存活和傳播[21]。SELL是一種細胞表面黏附分子,屬于黏附歸巢受體家族。該基因產物在白細胞與內皮細胞之間的結合和隨后的滾動過程中起著必要的作用,從而促進其遷移到次級淋巴器官和炎癥部位,同時還可以限制中性粒細胞的激活[22]。CXCL4可以改變單核細胞分化的軌跡,誘導包括Ⅱ類主要組織相容性復合體轉錄激活因子(CIITA)在內的關鍵轉錄調節因子介導的新型促炎和促纖維化表型[23]。

10個關鍵基因中,已有9個被報道[24-25]。與UC密切相關的基因包括CXCL8、TLR2、ICAM1、CXCR4、CTLA4、CD69、CXCL13、TIMP。此外,相關研究[26]通過隨機森林算法和Limma分析得出, BGN可能是UC加重的基因特征。SELL與UC相關的報道較少見。研究[27]發現,SELL是牙周炎與非酒精性脂肪性肝病的共同關鍵基因。SELL還被發現是卵巢癌三級淋巴結構的相關特征基因,可作為預測卵巢癌預后和指導免疫治療的生物標志物[28]。因此,進一步研究上述基因在UC中的作用,可以作為認識UC關鍵調控基因的新角度。

miRNA能夠調節人類基因的表達,在自身免疫性疾病的發病機制中起著至關重要的作用[29]。糖皮質激素(GC)是治療中度至重度活動性UC最有效的方法,研究[30]發現, miR-32-5p在GC耐藥組中顯著下調,熱休克蛋白90β家族成員1(HSP90B1)受到了miR-642a-5p和miR-150-5p 的共同調控。相關研究[31]發現并驗證了長鏈非編碼RNA(lncRNA)MEG3的升高能夠抑制UC大鼠中miR-98-5p的上調,并且通過調控miR-98-5p可以促進白細胞介素-10(IL-10)的表達。CHEN Z P等[32]認為,miR-146a-5p通過上調環指蛋白8(RNF8)和抑制Notch1/mTORC1通路減少腸上皮屏障損傷。PATHAK S等[33]發現,炎癥介質誘導UC患者腸道肌成纖維細胞(IMF)中miR-155的表達,通過下調細胞因子信號抑制蛋白1(SOCS1)的表達。WU C P等[34]發現, Toll樣受體4(TLR4)是hsa-miR-375的直接下游,其水平由has-miR-375負介導。hsa-miR-146a-5p、has-miR-155-5p、hsa-miR-335-5p、hsa-miR-32-5p、hsa-miR-98-5p、hsa-miR-93-5p及其靶基因已有與UC發生相關的報道[35], 本研究所得出的其余miRNA包括hsa-mir-92a-3p、hsa-mir-26b-5p、hsa-mir-4426、hsa-mir-4462b、hsa-mir-4647、hsa-mir-92b-3p及其靶基因可能在UC中發揮著關鍵的作用。

本研究使用WGCNA篩選差異基因,并且結合GSEA基因集富集方法和建立基因-miRNA調控網絡,以探索UC潛在的調控機制,識別新的潛在生物標志物和通路,具有一定的創新性。基因-miRNA調控網絡在UC病理生理學中起著重要作用。本研究的不足之處在于僅通過生物信息學分析微陣列表達譜,而沒有通過相關實驗進行驗證以及探索Hub基因和miRNA調節UC的詳細機制,因此,還需擴大臨床研究樣本和進一步實驗以闡明機制。