基于多層基因網絡的關鍵基因識別算法

魏丕靜,劉晶晶,趙永敏,蘇延森,鄭春厚

(1．安徽大學 物質科學與信息技術研究院,合肥 230601;2.安徽大學 計算機科學與技術學院, 合肥 230601;3.安徽大學 人工智能學院,合肥 230601)

疾病的關鍵基因是指在人體系統中與某種疾病密切相關的一組基因,其在人類生理過程和疾病發生過程中具有不可忽視的調控作用,了解關鍵基因在疾病中的功能和作用,對研究疾病的調控方式、復雜通路、治療和預后等具有重要意義。關鍵基因可用于診斷疾病、判斷疾病分期、預測和評價新藥或新療法的有效性等。藥物靶向治療的關鍵是藥物能特異性作用于疾病相關基因位點,故識別與疾病緊密相關的關鍵基因十分重要。但由于基因數量龐大,僅通過生物實驗的方法測定基因功能將會耗費巨大的時間成本和經濟成本。因此,基于計算模型識別疾病關鍵基因的預測算法亟待開發。目前有很多研究致力于發現疾病關鍵基因,此方面研究有助于探索人類復雜疾病的內部發病機制、研究疾病細胞存活所需的最小基因集和后續對疾病的治療方式及治療藥物的研究[1]。

研究表明,基因并非獨立的發揮生物作用,基因之間的相互作用普遍存在,并通過相互作用共同維持著生物內部整體環境的穩定性[2],故基于基因相互作用網絡的基因排序技術得到了廣泛的應用。Wang等提出了在蛋白質-蛋白質相互作用網絡上基于邊緣聚類的關鍵蛋白識別方法,該方法認為節點的重要性由節點與相鄰節點之間的相互作用的邊緣系數與聚類系數之和決定[3]。Fan等人提出了關鍵蛋白質預測方法,該方法將亞細胞室信息與基因表達信息相結合,并運用修改后的PageRank算法獲得加權蛋白質-蛋白質相互作用網絡,實驗結果表明其有更好的關鍵蛋白質預測性能[4]。由此可以看出,將網絡拓撲信息和生物學信息結合為研究關鍵基因提供了很好的思路。然而,雖然目前有多種技術可以用來識別疾病基因,但是大部分方法往往都是通過整合多個樣本構建基因共表達網絡,弱化了疾病樣本與正常樣本之間的差異信息,忽略了疾病樣本的個體特異性。此外,個體特異性網絡構建思想在揭示疾病的個體特征方面已經得到有效的驗證[5-6]。

本文以基因間表達相似性為基礎構建基因網絡,并用來篩選有價值的生物標志物或關鍵基因,探索基因和疾病之間復雜關系。具體來說,首先利用正常樣本的基因表達數據構建參考基因共表達網絡,然后依次將每個疾病樣本的基因表達數據與正常樣本組合,構建疾病樣本擾動網絡,根據此擾動網絡和參考網絡,得到每一個疾病樣本的個體特異性網絡。然后將個體特異性網絡作為單層網絡,并將單層網絡之間的基因聯系起來,從而得到多層基因網絡,這樣既保留了疾病樣本的特異性又將多個疾病樣本聯系在一起。最后,利用Wu等[7]提出的基于張量的多層網絡中心性的計算方法,對多層網絡中的基因節點中心性進行打分,從而得到關鍵基因集。與其他經典算法的對比分析表明該方法在預測藥物靶標基因上具有一定的優勢,功能和通路富集分析證明關鍵基因集與疾病聯系緊密。

1 數據與方法

1.1 數據集

基因表達數據集來源于基因表達綜合數據庫GEO(https：//www.ncbi.nlm.nih.gov/geo/)。本文主要考慮樣本量偏少的數據集,因此從GEO數據庫中獲取哮喘疾病的基因表達數據集GSE31773和GSE43696。在哮喘疾病樣本選取的過程中,由于mRNA在CD8+T細胞中的表達差異性大于在CD4+T細胞中,因此選擇的疾病樣本為CD8+類型的。此外,根據控制變量的原則,盡量使得正常樣本和異常樣本的其他生物信息如年齡,性別等保持一致。因此,在GSE31773中選取了8個正常樣本和6個疾病樣本,每個樣本包含8 789個基因。同理,在GSE43696中選取20個正常樣本和6個疾病樣本,每個樣本包含9 194個基因。

疾病相關的基因來源于DisGeNet(https：//www.disgenet.org/)和Phenopedia(https：//phgkb.cdc.gov/PHGKB/startPagePhenoPedia.action)數據庫。從兩個數據庫中獲取與哮喘相關的2 712個基因,并與GSE31773和GSE43696數據集中的數據進行整合,分別得到2 522個基因和2 478個基因的表達數據。

此外,從TTD(http：//db.idrblab.net/ttd/)數據庫獲取11個針對哮喘已獲批準的藥物靶標。

1.2 方法

1.2.1 多層基因網絡構建

多層基因網絡構建主要分為四步,具體構建過程如圖1所示。

第一步是獲取疾病相關基因的表達數據。首先從GEO數據庫獲取正常樣本和疾病樣本的基因表達數據,從疾病基因相關數據庫獲取所有與所要研究的疾病潛在相關的基因,從正常樣本和疾病樣本的表達數據中篩選出疾病相關基因的表達數據。

第二步是利用所有正常樣本構建參考基因網絡[5]。設參考網絡為Gref(V,E,W),其中點集V是由與疾病相關的基因所構成,邊集E表示基因對之間的邊集,W表示邊權,即基因對間的皮爾遜相關系數,其計算方式如式(1)。

ω(ij)=

(1)

其中,Cik表示基因i在第k個正常樣本中的表達值,n為基因節點的總數。

第三步是針對每個疾病樣本構建個體特異性網絡[5]。個體特異性網絡的構建參考Liu等[5]提出的方法。具體而言,首先在所有正常樣本的表達數據中加入一個疾病樣本的表達數據,根據第二步的公式(1),求新的表達數據中基因之間的皮爾遜相關系數,構建一個新的基因網絡,并將其看作是加入該疾病樣本后的擾動網絡[5]。然后根據參考網絡和擾動網絡構建個體特異性網絡,其中邊權值為擾動網絡和參考網絡的邊權值的差值絕對值[5]。接著利用拐點分析法設置閾值對網絡中的邊進行選擇,刪除一些不顯著的邊。這種方法考慮到了每個樣本的個體特異性,體現了參考網絡受到疾病樣本的干擾程度,有效衡量了基因間相互作用關系與疾病的相關程度。

圖1 多層基因網絡構建示意圖Fig.1 Schematic diagram of multilayer network construction

注：多層基因網絡構建分成四步：第一步是數據收集,即在GEO數據庫中選擇正常樣本和疾病樣本,在疾病基因數據庫中選擇出與哮喘相關的基因;第二步是利用正常樣本構建參考網絡;第三步是構建疾病樣本的個體特異性網絡,首先利用疾病樣本構建擾動網絡,再用擾動網絡減去參考網構建個體特異性網絡,然后根據閾值去除部分異常邊權值后得到最終的個體特異性網絡。第四步是整合所有個體特異性網絡得到多層基因相互作用網絡,層間邊連接每層網絡中的相同基因,邊的權值為1.

第四步是整合單層基因網絡得到多層基因網絡。將得到的個體特異性網絡作為多層基因網絡的每一層,依次連接每兩個單層網絡中的相同節點構建層與層之間的邊,邊的權值為1,得到多層復用基因網絡。

1.2.2 基因節點中心性計算

(2)

根據單層網絡中PageRank算法的冪法求解過程,交互張量H相當于轉移概率矩陣,求解張量方程HΦ=λΦ得到中心性二階張量Φ,其中Φiα表示當前迭代中第α層的第i個基因節點的中心性值,λ表示特征系數,λ這里取值為1,保證二階張量Φ的存在性和唯一性。迭代結束后,將每個基因節點在所有層中的中心性均值作為該基因的最終中心性值,降序排序后選取排名靠前的基因作為關鍵基因,分值越高說明基因在疾病中發揮的作用越重要。

2 結果與討論

2.1 網絡邊閾值的選取

通過設置皮爾遜相關系數的閾值得到多層網絡。具體而言,針對GSE31773和GSE43696兩個數據集,分別利用拐點分析法選擇拐點,并將其作為篩選邊的閾值。根據圖2可以發現,GSE31773數據集拐點示意圖中,當邊權值大于1時,趨勢不再有明顯上升,因此構建網絡的閾值選擇為1。同理對于數據集GSE43696閾值選擇為0.6。確定數據集GSE31773和GSE43696構建6層網絡的閾值分別為1和0.6。

圖2 拐點分析圖Fig.2 Analysis diagram of inflection point

2.2 多層基因網絡有效性分析

利用本文的方法,針對兩個獨立數據集GSE31773和GSE43696分別構建多層網絡,其信息如表1所示,其中層間邊連接每層的相同基因,例如數據集GSE43696,其中層間的邊數是每層節點連接其他五層中相同節點,即總計37 170條邊。以數據集GSE43696構建的多層網絡為例,將其可視化后如圖3所示。

表1 多層網絡信息Table 1 Information of multilayer network

圖3 GSE43696： 6層基因網絡示意圖Fig.3 GSE43696： Diagram of 6-layers gene network

在多層網絡構建過程中,多層基因網絡的層數有多種選擇。為了驗證本文構建6層網絡的有效性,本文在數據集GSE31773中隨機選擇不同數量的疾病樣本構建了不同層數的多層基因網絡,并對比已知哮喘藥物靶標基因在不同層數的網絡中的排名結果,如表2所示。其中在選擇哮喘藥物靶標基因時,首先選擇有效治療哮喘的藥物,并在數據庫中尋找藥物關鍵基因靶標,最終選擇包含在本文數據集中的11個靶標基因。表2中“排名1”和“排名2”指隨機選擇了兩次相同數量樣本的結果。從表中看出,在六層基因網絡中,有5-LOX、IL17、CCR4、IL5RA、ROS等5個哮喘的藥物靶標基因排名更靠前;在五層基因網絡中,有H1R、IL5、JAK-1等3個基因排名更優,在四層基因網絡中,只有基因CAMP有更好的排名;在三層基因網絡中,有2個基因JAK-2、IL4R排名更優。綜上所述,在識別關鍵基因集時,構建六層網絡的效果更好。

表2 哮喘靶標在不同層網絡排名情況Table 2 Ranking of asthma targets in different layers of networks

2.3 算法對比分析

根據哮喘基因數據集,利用本文提出的方法,可以得到哮喘相關的基因。為了進一步評估已知的疾病特異性通路或基因是否在預測的關鍵基因上具有顯著的優先級,本研究利用GSEA軟件的GSEAPreranked工具對結果進行分析。GSEA富集分析主要是用來評估一個預先定義的基因集在與表型相關的基因排序列表中的分布趨勢,它不需要進行基因過濾,輸入數據主要包括兩部分,一種是預先定義的基因集,一種是給定的基因排序列表。本文中,預先定義的基因集是KEGG通路數據庫中的哮喘特異性相關的基因集,基因排序列表是本文預測的所有基因排序結果。通過GSEA富集分析揭示我們的模型結果和KEGG通路數據庫中哮喘特異性相關的基因集之間的關聯,以GSE31773數據集為例,根據其所有基因排名和KEGG通路數據庫中哮喘特異性相關的基因集進行加權K-S檢驗得到p值,結果如圖4所示。結果表明,與其他預測關鍵基因的排序方法MI[8],t-Test[9],PCC[10],SCC,FC[11],NetRank[12],MarkRank[13]相比,本研究中的算法在對疾病關鍵特異性基因進行優先排序時具有顯著的p值。

圖4 不同對比算法的哮喘通路富集分析Fig.4 Enrichment analysis of asthma pathways withdifferent comparison algorithms

2.4 哮喘關鍵基因集分析

為了驗證本方法所識別的疾病相關基因的重要性,針對GSE31773和GSE43696兩個數據集,分別選擇排名前10的關鍵基因(見表3),分析是否已有研究證實其為哮喘關鍵基因。研究發現TP53、MAP3K1、COL18A1、DACT1、CD40LG、ANKRD55、CD4以及TNFSF18、AFM、NKX2-1、SCGB1A1、RAG1、FRAS1、HSD11B2、GSTO2、SOAT1、IL19等基因在哮喘發生發展過程中起重要作用。例如,Yuan等[14]的研究表明,與遲發性哮喘臨床表型相關的TP53差異甲基化位點是早期篩選的有效生物標志物。Zhang等的研究證明DACT1可能是治療哮喘的潛在靶點[15]。對于CD40LG,有研究表明CD86和CD40LG之間的相互作用會促進過敏性哮喘的發展[16]。CD4T細胞淋巴細胞活化在嚴重哮喘發病機制中起重要作用[17]。SCGB1A1是肺重要的防御分子,防止SCGB1A1被抑制可有效的改善哮喘[18]。有研究表明GSTO2是哮喘易感基因,GSTO2基因的多態性和哮喘有關[19]。此外,有研究證實,IL-19基因在哮喘中高度表達,在變應性疾病中起著重要作用[20]。研究還發現,嗜酸性粒細胞的凋亡在支氣管哮喘病理生理中發揮至關重要的作用[21],并且PUS10基因能夠調節Trail誘導的細胞凋亡過程[22]。軸突或突觸結構調控哮喘的激發機制[23],且F5蛋白在膜-細胞骨架相互作用和突觸結構或功能的動態方面發揮重要作用[24]。此外,“下丘腦-垂體-腎上腺”軸功能與肺功能改善程度相關[25],并且CRHBP調節促腎上腺皮質激素控制“下丘腦-垂體-腎上腺”軸功能[26]。由此推斷,PUS10、F5及CRHBP等基因也與哮喘發生發展緊密相關。

2.5 GO功能富集分析

為了分析本算法預測的疾病關鍵基因的功能相關性,利用本算法分別在兩個獨立數據集上選擇排名前100的基因,使用基因功能分析工具DAVID對其作GO功能富集分析?；贒AVID分析工具,得到與前100個基因顯著相關的基因本體,圖5展示了排名前10的基因本體。圖的縱坐標展示了GO的功能注釋,橫坐標上的值表示GO在關鍵基因集中的富集顯著性值-log(p)。由圖可以發現,在排名前10的基因本體中,免疫反應、調控T細胞增殖、T細胞刺激以及細胞因子活性均被證實與哮喘有密切聯系[38]。具體來說,哮喘是由免疫系統對環境因子和不同的基因表達的聯合反應引起的呼吸系統疾病。T細胞是哮喘中過敏性氣道炎癥的關鍵介質[39],T細胞的增殖會引起免疫球蛋白水平增加和支氣管高反應性即哮喘發作,細胞因子也會輔助T細胞增殖的反饋控制。此外,炎癥反應也與哮喘相關,在哮喘惡化過程中伴隨著循環嗜酸性粒細胞、嗜堿性粒細胞及其前體細胞的變化等各種炎癥反應[40-41]。除上述機制外還有幾種潛在的新機制,例如藥物反應,內皮細胞分化,蛋白質磷酸化調控,信號調控,應對缺氧,轉錄調控等在哮喘發展過程中都起著重要的作用。

表3 排名前10的關鍵基因集Table 3 Top 10 critical gene sets

圖5 關鍵基因富集Gene OntologyFig.5 Key gene enrichment Gene Qntology

2.6 通路富集分析

為了定位關鍵通路的關鍵基因,本文基于DAVID平臺對兩個獨立數據集排名前100的基因進行通路富集分析,得到與100個基因顯著相關(p_val≤0.05)的通路,表4和表5展示了顯著相關的通路。由上述通路富集分析結果可知,細胞因子受體相互作用、趨化因子信號通路、T細胞受體信號通路、原發性免疫不全四條通路都與哮喘緊密相關。腫瘤壞死因子(TNF)信號通路、TGF-beta信號通路、Th1/Th2分化等通路也被證明與哮喘有關。TNF信號通路的壞死因子α是免疫和炎癥反應的有效調節劑,可以引起包括哮喘在內的多種自身免疫性疾病[42]。哮喘會通過TGF-beta信號通路促進小鼠脈絡膜血管新生[43]。T淋巴細胞介導的對過敏原的免疫應答是哮喘發病機制的早期關鍵因素,而Th1/Th2平衡是哮喘發病機制的核心[44]。此外,還有若干個與哮喘潛在相關的通路,包括黏著連接、焦點粘連、鞘脂類信號通路等。

表4 GSE31773：關鍵基因通路富集分析Table 4 GSE31773： Pathways enrichment analysisof critical genes

表5 GSE43696：關鍵基因通路富集分析Table 5 GSE43696： Pathways enrichmentanalysis of critical genes

3 結 論

1)復雜疾病的發生發展本質上與基因和生物功能過程的改變密切相關,疾病關鍵基因的識別對于研究疾病機理尤其是藥物靶向治療具有重要意義。哮喘作為全球范圍內發病率最高的慢性呼吸道疾病之一,其發病率在逐年上升。識別出與哮喘成因緊密相關的基因有助于提高治療效果。然而臨床研究中由于疾病樣本數較少,通常導致疾病相關基因識別困難。針對上述問題,本研究提出基于少數樣本構建多層網絡,進而利用多層網絡隨機游走識別疾病相關的關鍵基因的方法。該方法有助于挖掘樣本數量受限條件下的疾病相關基因,加深對疾病致病機理的理解。

2)構建的多層網絡對識別小樣本疾病的致病基因可行且有效。本文利用皮爾遜相關系數計算出每條邊的權值;為增強網絡結構的穩定性,采用拐點分析法尋找最佳閾值,保留擾動程度較大的邊;通過比較對已知疾病關鍵基因的排序選取最優的網絡層數。例如針對數據集GSE31773的實驗分析表明,構建六層基因網絡效果最佳。

3)與其他方法相比,本算法識別的哮喘相關基因的排名更具顯著性。利用本算法分別在GSE31773和GSE43696數據集中挖掘排名前10的關鍵基因,研究發現TP53、MAP3K1、COL18A1、DACT1、CD40LG、ANKRD55、CD4以及TNFSF18、AFM、NKX2-1、SCGB1A1、RAG1、FRAS1、HSD11B2、GSTO2、SOAT1、IL19等基因在哮喘發生發展過程中起重要作用,并推斷PUS10、F5及CRHBP等基因也與哮喘發生發展緊密相關。

4)對分別從GSE31773和GSE43696兩個數據集中所得排名前100的關鍵基因進行通路富集分析和GO功能富集分析,分析結果表明所識別的基因能夠顯著富集到與哮喘相關的通路和功能中。