基于生物分子網絡分析的精神分裂癥功能模塊挖掘

廖苑君,陳應堅,趙小蕾,孫勝南,林 帆,覃繼恒,饒紹奇*

(廣東醫科大學a.公共衛生學院;b.醫學系統生物學研究所,中國廣東東莞523808)

精神分裂癥(schizophrenia,SCZ)是一種嚴重影響人類生命健康和生存質量的復雜性精神疾病,患病率約為1%,遺傳貢獻高達83%[1]。患者可出現幻覺、妄想、情感淡漠、社交退縮等典型精神癥狀,同時自知力喪失或不完整,造成嚴重的思維、情感、行為等多方面的障礙[2]。SCZ病癥具有較高的復發率和精神癥狀殘留,且具有不可預知性,給家庭和社會帶來巨大的經濟負擔,并占據了大量的醫療資源,導致資源利用價值下降。在此背景下,研究SCZ病因與病理機制對其預防及治療具有重要的意義。目前,大規模組學數據的發展為研究遺傳學機制提供了便利,然而分析方法主要采用的是單位點、單基因分析策略[3]。由于SCZ病因復雜,其發生、發展過程涉及多個生物分子(如基因、miRNA等),這些生物分子構成了一個復雜的網絡,因此與單位點、單基因分析方法相比,生物分子網絡分析能夠更加全面揭示SCZ的分子機制[4]。

在生物分子網絡分析中,節點代表生命體內的各種分子,邊(網絡連接)代表分子間的互作關系。根據對象的不同,生物分子網絡可分為基因調控網絡、蛋白質互作網絡、信號轉導網絡等[5]。迄今為止,生物分子網絡分析在生物醫學領域已被廣泛應用,具有成熟的分析方法,可用來篩選和挖掘網絡中關鍵的核心節點及功能模塊,從而提示具有重要生物學意義的生物聯系或功能結構。Xu等[6]發現基于功能注釋定義的功能模塊不僅可以用于癌癥的分型或預測,甚至能夠挖掘癌癥中潛在的亞型。Zhao等[4]通過開展冠心病風險通路內、通路間網絡研究,挖掘了冠心病與阿爾茨海默病共享的功能模塊。本研究借助網絡分析理論,深入挖掘與SCZ相關的功能模塊,以揭示潛在的分子致病機制,為SCZ的基礎研究提供新思路。

1 材料與方法

1.1 數據來源

作為SCZ遺傳公共數據庫之一,SZDB數據庫[7](http：//www.szdb.org/)收錄了大量的數據集和候選基因。對SZDB數據庫進行挖掘,納入標準：1)基因表達數據集;2)數據集的患者病例數超過10例,經篩選得到5套完整數據集[8～12]。利用該數據庫提供的在線處理分析工具篩選差異表達基因。首先對數據進行標準化(歸一化)預處理,隨后利用R和Bioconductor平臺中的Limma包進行差異表達分析。考慮到基因表達豐度低的mRNA的差異表達倍數不一定足夠大,本研究不設置差異倍數閾值。為了更好地控制差異表達結果和假陽性率,將在3套或以上數據集中具有顯著性差異表達(P＜0.05)的基因納入后續分析。最終共挖掘了944個與SCZ顯著相關的基因,并定義為初始種子基因集。

HPRD 數據庫[13](http：//www.hprd.org/)是人類蛋白質參考數據庫,含有人類蛋白質互作注釋信息。迄今為止,該數據庫儲存了41 327對蛋白質互作(protein-protein interaction,PPI)對子數據,數據可信度高,證據可靠。本研究提取該數據庫的PPI對子數據,構建致病基因網絡。

1.2 致病基因網絡的構建

本研究構建的網絡節點集合包含了初始種子基因集及其在PPI的一級鄰居基因,對于集合中的任意兩個不同節點,假設其產物在PPI網絡中存在互作關系即視為網絡中連接。因此,構建的網絡可視為SCZ特異性的風險基因網絡,其中節點代表基因,節點間的邊代表基因間相互作用。

1.3 網絡功能模塊的挖掘

生物分子傾向于協同執行一定的生物過程或發揮特定生物功能。同樣地,在生物分子網絡中,其對應的節點連接高度緊密,聚合為節點簇,構成網絡的模塊化(modularity)。因此,與同樣規模的隨機網絡相比,模塊化網絡連接更為緊密。從功能模塊的角度分析網絡,能夠將生物分子進行功能分區的劃歸,對預測其功能、理解生命活動的層級和結構具有關鍵提示作用。功能模塊的結構劃分評價可看作是網絡結構分解的優化問題。基于研究假設,本研究通過Newman算法[14～15]對構建的網絡進行分解,從而識別網絡中高度模塊化的功能模塊。該算法將網絡分解問題轉變為二次型優化問題,并將之簡化為特征向量提取問題。Newman算法擬優化的目標函數為其中,當Aij=1時說明節點i和節點j之間存在互作關系,反之取值為0;m代表網絡中邊的總數;ki和kj分別代表節點i和j的連通度;s代表有n個節點指示變量的列向量,向量內元素的取值為1或-1,si和sj分別代表節點i和j所屬的功能模塊,即si=即兩節點位于同一功能模塊時,Q才能取得最大值。為了方便,將等式轉換為二次型的形式：Q=為s的倒置,由列向量變為行向量;B代表著一個內部元素為Bij=Aij-kikj/2m的對稱矩陣,當s取B的最大特征值對應的特征向量時,Q值最大。因為s只可取值為±1,所以s只保留矩陣特征向量對應的正負號,通過符號方向將網絡分解為兩部分。通過以上方法重復網絡一分為二的過程,直到網絡不能分解為止。最終將這些稠密、不可分割的子網提取出來,并定義為功能模塊。本研究設定的功能模塊節點數大于50,上述算法通過R語言“igraph”軟件包完成。

1.4 網絡拓撲屬性分析

對網絡的拓撲屬性進行評價時,評價指標包括以下五方面的內容。

1)網絡直徑(network diameter)：即網絡中任意兩個連通節點間最短路徑的最大值,可反映網絡的大小。

2)連通度(degree)：某個節點的連通度是指網絡中直接與該節點相連的邊數,一般用符號k表示,能夠衡量該節點在維持網絡完整性中的作用。在隨機網絡中,連通度的分布符合泊松分布[16],通過檢驗某個節點的連通度是否顯著大于一般節點,可進一步篩選和挖掘核心節點。基于理論假設,某個節點的連通度大于或等于t的概率為：其中 λ=NP,λ 代表節點的期望連通度,N代表節點數目,P代表任意兩個節點發生連接的概率。通過Bonferroni校正,以P＜0.01為差異有統計學意義。當P值小于預設檢驗水準時,該節點即定義為核心節點。

3)介數(betweenness)：描述網絡中某節點的路徑占其他節點間最短路徑的比例。介數越高,表示其在維持網絡緊密連接性方面起更加重要的作用。

4)聚類系數(clustering coefficient)：衡量網絡中某節點與周圍節點的密集連接程度。

5)無標度網絡(scale-free network)：指由少數高連通度節點及多數低連通度節點組成的網絡。一般認為,無標度網絡中連通度的分布符合冪律分布,即P(k)∝k-α,其中k代表接通度,α代表指數參數,使用極大似然估計法[17]。本研究采用KS(K-olmogorov-Smirnov)擬合優度檢驗評判無標度特性。

1.5 功能富集分析

通過R語言(版本3.5.0)“clusterProfiler”軟件包[18]對每個功能模塊進行KEGG(Kyoto Encyclopedia of Genes and Genomes)功能富集分析。通過FDR校正,以P＜0.001為差異有統計學意義,并針對每個功能模塊最顯著的前5條富集通路開展后續分析。根據富集到的KEGG通路類別評估功能模塊之間的相互作用,利用Cytoscape軟件(版本3.6.1)構建功能模塊與通路分類的關系。

本研究遵循上述流程開展,技術路線見圖1。

2 結果

2.1 網絡的構建與拓撲屬性

通過PPI引導,構建的SCZ特異性基因網絡含有3 207個節點和5 207個互作對子。由于網絡中的節點之間并非全部緊密連接,故將孤立作用的節點(TRIM63、PDK4、DLAT 等)去除掉,只對最大子網進行后續分析,最終生成的網絡包含3 062個節點和5 120個互作對子。通過網絡拓撲屬性分析得到網絡直徑為12,聚類系數為0.034,網絡平均路徑長度為4.972。

2.2 網絡功能模塊的挖掘與拓撲屬性分析

通過Newman分解算法對最終生成的網絡進行分解,得到14個節點數大于50的功能模塊。14個功能模塊的拓撲屬性和無標度檢測結果如表1所示。從表1可知,功能模塊內的節點數差別很大。例如：最大的功能模塊(M9)由522個節點和860條邊組成,每個節點平均含有1.6條邊,而最小的功能模塊(M14)由51個節點和52條邊組成,每個節點平均含有1條邊。從直徑的分布結果來看,功能模塊尺寸越大,其直徑反而呈下降的趨勢,這表明隨著功能模塊尺寸的增加,節點間聚集的程度會越發緊密。這些功能模塊的拓撲性質指標接近,直徑數值在8～15波動,特征路徑長度均小于6,呈現“小世界”的特性,符合六度分離理論,支持了生物學上無標度網絡的常見表現[19]。無標度檢驗結果顯示,各功能模塊指數參數估計值均為2～3,連通度均近似符合冪律分布(KS檢驗,P＞0.05),支持了功能模塊的無標度性質,可視為無標度網絡。

圖1 數據獲取及分析流程圖Fig.1 Flow chart of data acquisition and analysis

2.3 核心基因的識別

利用泊松分布檢驗篩選出102個核心基因(P＜0.01)。通過檢索PubMed數據庫,得到57個核心基因在已有文獻中被報道與SCZ相關,其余的45個核心基因暫無文獻支持其與該疾病相關(表2)。在14個功能模塊中,M9含有的核心基因最多,多達21個;M7含有的核心基因最少,僅有1個。圖2是我們利用R語言軟件展示的部分功能模塊與核心基因分布圖。

表1 各功能模塊的基本特征和拓撲屬性Table 1 Basic characteristics and topological properties of each functional module

根據連通度對基因進行排序,列出連通度較大的前10個基因的拓撲參數(表3),結果顯示EGFR基因具有最高的連通度,維系著網絡中15.29%的基因間最短路徑的連通路徑,提示其在網絡中發揮重要作用。

2.4 模塊的功能富集分析

KEGG通路富集分析結果顯示,模塊內的基因顯著富集于細胞凋亡(hsa04210)、ErbB信號通路(hsa04012)、細胞周期(hsa04110)、磷脂酶D信號通路(hsa04072)、PI3K-Akt信號通路(hsa04151)等多個有統計學意義的信號通路(表4)。根據富集到的KEGG通路類別評估這14個功能模塊之間的相互作用,利用Cytoscape軟件構建功能模塊與通路分類的關系,結果如圖3所示。圖3信息表明,提取的SCZ功能模塊主要參與折疊、分類和降解(folding,sorting and degradation)、細胞通訊(cellular communication)、免 疫 系 統(immune system)、傳染病(infectious diseases)、癌癥(cancers)及內分泌系統(endocrine system)等多個功能類別。

3 討論

表2 每個功能模塊的核心基因Table 2 Hub genes of each functional module

表3 連通度排名前10的網絡節點的拓撲屬性Table 3 Topological properties of some network nodes

近年來隨著基因組學和蛋白質組學的發展,眾多研究機構利用組學實驗手段產生了大量與疾病相關的實驗數據。利用公共數據庫資源或整合實驗數據,能夠幫助研究者抽取有重要生物學意義的生物信息,重構不同層次的生物分子網絡。從功能模塊的角度分析網絡的功能學特征,能夠挖掘有意義的生物分子或功能結構。為此,本研究基于網絡模型,利用SZDB和HPRD公共數據庫構建SCZ特異性的風險基因網絡,并進一步通過Newman算法對構建的網絡進行分解,挖掘了緊密聯系的功能模塊和核心基因,分析了功能模塊行使的生物學功能。

圖2 部分功能模塊與核心基因分布圖含有標簽的節點代表對應功能模塊的核心節點,節點的大小代表該基因的連通度大小。Fig.2 Part of functional modules and hub gene distributionThe node containing the tag represents the hub node of the corresponding functional module,and the size of the node represents the degree of the gene.

本研究共挖掘了14個與SCZ相關的功能模塊,并通過拓撲學指標評價各模塊的拓撲性質,結果顯示所有功能模塊的拓撲屬性接近,符合六度分離理論,具有小世界性,支持了生物學上無標度網絡的常見表現。KS擬合優度檢驗結果顯示,所得功能模塊的連通度均近似服從冪律分布,符合無標度性。無標度性是指功能模塊包含了極少數連通度較高的核心基因,它們位于功能模塊的中樞,發揮“信息樞紐”的作用,能夠維持其功能的穩定性,而這些核心基因一旦受到破壞將會影響功能模塊的功能,進而導致疾病的發生發展。利用泊松分布檢驗篩選得到102個核心基因,通過檢索PubMed數據庫,發現共有57個核心基因已被報道與SCZ相關,例如EGFR、HAX1、IL1-R1、RALGDS、NSF等基因。其中,EGFR基因在網絡中具有最高的連通度,維系著網絡中15.29%的基因間最短路徑的連通路徑。EGFR是具有酪氨酸激酶活性的表皮生長因子受體家族(epithelial growth factor receptor family)成員。相關研究發現,EGFR能夠抑制多巴胺D3樣受體的信號傳導,是SCZ潛在的治療靶點[20]。除此之外,一些尚未有報道證實與SCZ相關的核心基因也在研究過程中被發現,例如 SVIL、DNAJA1、RABAC1、STX6 等基因,雖然這些基因仍未經過完善的研究和驗證,但從本研究的結果來看,它們可能會成為治療SCZ的重要靶標。

表4 各功能模塊的KEGG通路分析結果Table 4 KEGG pathways of each functional module

通過對所得功能模塊進行KEGG信號通路的富集分析,發現大部分通路已有文獻支持與SCZ的發病機制相關,例如：TGF-β1(hsa04350)和Hippo信號通路(hsa04390)在誘導Th17細胞的分化(hsa04659)過程中起著積極作用,進而介導組織的炎癥反應[21～23];血管緊張素轉換酶是腎素-血管緊張素系統(hsa04614)的關鍵酶,可以催化神經肽的降解并調節多巴胺能和5-羥色胺能神經傳遞[24];長時程抑制(hsa04730)作為突觸可塑性的主要形式之一,在學習和記憶形成的機制中發揮重要作用[25];FcγR介導的吞噬(hsa04666)過程能夠誘導肌動蛋白細胞骨架的組裝(hsa04810),從而影響細胞的運動、分裂、黏附及吞噬過程,目前研究發現精神疾病的發病與神經突觸被過度吞噬密切相關[26];腫瘤抑制因子PTEN對PI3K-Akt信號通路(hsa04151)起到負性調節作用,導致SCZ患者的認知功能受到影響[27～28]。同時,我們也得到了未明確證實與SCZ發病機制相關的通路,例如Ras信號通路(hsa04014)。據報道,Ras基因的失調能夠導致5-羥色胺系統和多巴胺系統受損[29];Ras信號通路(hsa04014)在神經細胞的生長、分化等多個生物學過程中具有重要作用[25],推測其與精神疾病的發病密切相關。

本研究根據富集到的KEGG通路類別探尋這14個功能模塊之間的相互作用,得到了更為深入的功能聯系。結果顯示,得到的功能模塊主要參與折疊、分類和降解、細胞通訊、免疫、傳染病、癌癥及內分泌系統等多個功能類別。其中,信號轉導(M2,M3,M6,M9,M12)、細胞通訊(M1,M9,M11,M12)以及折疊、分類和降解(M2,M6,M13,M14)這3種功能類別的連通度較高,提示這些功能與SCZ的病理機制關系更為密切。除此之外,功能模塊M4、M7、M11共同參與免疫炎癥反應;功能模塊M3、M5、M7共同調控細胞周期;功能模塊M10主要富集于冠心病通路,提示功能模塊M10可能成為解釋SCZ和冠心病共享的關鍵功能模塊。功能模塊分析顯示,大部分功能模塊不是單獨發揮作用的,它們共同影響SCZ的發生發展,具有共享的病理機制。

圖3 功能模塊與通路分類的關系綠色圓圈表示功能模塊,紅色圓圈表示生物學通路對應的功能類別。Fig.3 Relationship between functional modules and pathway classificationThe green circle represents the functional modules,and the red circle represents the functional category corresponding to the biological pathway.

總的來講,本研究從功能模塊的角度分析了SCZ網絡的功能學特征,為疾病的基礎研究提供了新的思路。此外,研究過程采用的方法具有普適性和推廣性,能夠進一步推廣到其他疾病研究。需要指出的是,研究的預測結果在一定程度上受到數據庫的信息量完整性和正確性的影響,存在一定的局限性,在今后的研究中仍需要整合更多的數據信息來加以完善。