利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析方法識(shí)別帕金森病進(jìn)展過程中的關(guān)鍵基因

王樹文，樊冬梅，陳方正，謝玉平，呂曉彤，李依璘，林鴻程 （廣東醫(yī)科大學(xué)，廣東東莞 523808）

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)是探究基因復(fù)雜調(diào)控網(wǎng)絡(luò)最常用到的生物信息學(xué)分析方法[1-3]。帕金森?。≒D）作為一種涉及多基因調(diào)控的疾病，目前利用WGCNA 算法對(duì)PD 進(jìn)行研究的文獻(xiàn)報(bào)道不多。本研究通過構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)最終識(shí)別出PD 密切相關(guān)的關(guān)鍵基因模塊及關(guān)鍵基因，以期為闡明PD 發(fā)生發(fā)展機(jī)制及治療提供新的研究靶點(diǎn)分子。

1 材料和方法

1.1 PD數(shù)據(jù)集和樣本信息的獲取

PD相關(guān)表達(dá)譜芯片數(shù)據(jù)與樣本信息下載于GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22491)，芯片數(shù)據(jù)編號(hào)：GSE22491，芯片平 臺(tái)：GPL6480（Agilent-014850 Whole Human Genome Microarray 4x44K G4112F），檢測(cè)細(xì)胞類型:peripheral blood mononuclear cells。數(shù)據(jù)集樣本信息主要包括8例正常樣本和10例PD 樣本信息?；赗軟件中的“GEOquery”下載數(shù)據(jù)集矩陣文件和平臺(tái)注釋文件[4]，根據(jù)注釋文件對(duì)下載的矩陣文件進(jìn)行基因重注釋后，過濾掉注釋信息缺失的探針和低表達(dá)的探針，最后利用R 包“l(fā)imma”中的“normalizeBetweenArrays”函數(shù)對(duì)芯片表達(dá)譜數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理[5]。另外對(duì)矩陣文件中的樣本信息單獨(dú)整理出用于后續(xù)WGCNA 分析。為保證WGCNA 分析過程中運(yùn)行的順利性和可靠性，本研究選取了芯片表達(dá)譜中方差變化最大的前25%的基因用于下游分析。另外在正式進(jìn)行WGCNA分析前，還需要通過計(jì)算樣本間的相關(guān)性，繪制樣本的分層聚類圖，刪除離群樣本。

1.2 WGCNA分析

WGCNA 分析完成主要是基于R 軟件中的“WGCNA”包實(shí)現(xiàn)[6]，具體流程包括：根據(jù)表達(dá)譜計(jì)算出適宜的軟閾值，在計(jì)算出兩兩基因之間相互強(qiáng)度的基礎(chǔ)上，對(duì)其進(jìn)行加權(quán)計(jì)算并構(gòu)建出基于芯片表達(dá)譜的鄰接矩陣，隨后根據(jù)TOM 重疊計(jì)算公式將鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣(TOM)，根據(jù)聚類樹算法和1－TOM 計(jì)算公式，將表達(dá)相似基因合并在一起從而構(gòu)建出多個(gè)基因集模塊，同一模塊內(nèi)基因往往具有相似和相關(guān)的生物學(xué)功能。最后計(jì)算基因模塊的特征值(ME)，同時(shí)關(guān)聯(lián)樣本信息（正常與PD），使用“Pearson”相關(guān)系數(shù)計(jì)算模塊特征值與樣本性狀的相關(guān)程度，將與PD相關(guān)程度最強(qiáng)的模塊定義為PD關(guān)鍵基因模塊，隨后對(duì)該模塊內(nèi)基因進(jìn)行深入挖掘和功能分析。

1.3 關(guān)鍵基因的篩選

首先計(jì)算PD 關(guān)鍵基因模塊內(nèi)MM 值和GS 值，MM 值是基因表達(dá)水平和ME 值之間的相關(guān)系數(shù)，相關(guān)性越高，說明基因和該模塊的關(guān)聯(lián)性就越高，GS表示每個(gè)基因與性狀的相關(guān)程度，參考其他文獻(xiàn)我們將PD 關(guān)鍵基因模塊內(nèi)基因連接度最大的前10 個(gè)基因定義為PD 關(guān)鍵基因，后續(xù)通過文獻(xiàn)檢索，對(duì)PD 關(guān)鍵基因進(jìn)行進(jìn)一步篩選，識(shí)別出目前在PD 研究中還未涉及到的關(guān)鍵基因。

1.4 基因功能(GO)和京都基因與基因組百科全書(KEGG)富集分析

利用R 軟件中的“clusterprofiler”包[7]對(duì)PD 關(guān)鍵模塊內(nèi)基因進(jìn)行進(jìn)行GO 功能富集分析和KEGG 通路富集分析和可視化，分析結(jié)果中P＜0.05 表示富集結(jié)果的差異具有統(tǒng)計(jì)學(xué)意義。其富集結(jié)果以氣泡圖形式呈現(xiàn)。

1.5 PD關(guān)鍵基因差異分析與共表達(dá)分析

利用R 軟件中的“l(fā)imma”包進(jìn)行篩選出的10 個(gè)PD 關(guān)鍵基因在正常組織樣本與PD 樣本差異性，差異條件設(shè)置為：差異倍數(shù)絕對(duì)值≥1.2 且較正后的P＜0.001。最后利用R軟件中內(nèi)置的“cor”函數(shù)進(jìn)行10個(gè)PD關(guān)鍵基因的相關(guān)性分析，分析方法選用“Pearson”。

1.6 PD相關(guān)基因的篩選

基于Genecards 數(shù)據(jù)庫（https://www.genecards.org/）篩選目前比較明確的PD 相關(guān)基因，從而對(duì)本研究篩選出的10個(gè)PD關(guān)鍵基因進(jìn)行初步驗(yàn)證。

2 結(jié)果

2.1 芯片表達(dá)譜預(yù)處理及加權(quán)基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建

整理GSE22491 芯片表達(dá)譜，過濾其中探針信息注釋不全和重復(fù)的基因，最終獲得16 576 個(gè)基因，其中包括16 327個(gè)編碼基因和249個(gè)長(zhǎng)非編碼RNA，選取方差變化最大的前4 082 個(gè)基因進(jìn)行后續(xù)的WGCNA 分析。通過計(jì)算樣本間的相關(guān)性，繪制樣本的分層聚類圖，將聚類高度設(shè)置為90（紅線標(biāo)記），見圖1。刪除異常樣本GSM558687，剩余17 例樣本信息及相應(yīng)的4 082 個(gè)基因芯片表達(dá)譜數(shù)據(jù)合并用于后續(xù)WGCNA構(gòu)建模塊，見圖2。

圖1 樣本聚類并篩選異常樣本(紅線表示高度=90）

圖2 樣本聚類與對(duì)應(yīng)的樣本信息

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建

基于文獻(xiàn)中WGCNA 分析流程，借助R 軟件中“WGCNA”包實(shí)現(xiàn)本研究所需的加權(quán)基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建。在WGCNA分析過程中，軟閾值的選取對(duì)于構(gòu)建符合生物學(xué)意義的無尺度網(wǎng)絡(luò)具有決定性的作用，本研究基于“WGCNA”包中的“pickSoftThreshold”函數(shù)進(jìn)行了軟閾值的篩選，見圖3。當(dāng)選軟閾值即β=7 時(shí)，可以保證本研究構(gòu)建出的網(wǎng)絡(luò)既達(dá)到無尺度標(biāo)準(zhǔn)，也使網(wǎng)絡(luò)的平均連接程度不至于較小而缺乏足夠的信息。隨后計(jì)算基因間的相關(guān)性矩陣和鄰接矩陣，根據(jù)TOM 重疊計(jì)算公式將鄰接矩陣轉(zhuǎn)換為TOM 重疊矩陣，根據(jù)1－TOM 公式構(gòu)建基因間分層聚類樹，同時(shí)使用動(dòng)態(tài)剪切樹的算法把基因分成15個(gè)模塊，根據(jù)模塊相似度大于0.75標(biāo)準(zhǔn)合并相似模塊（圖4），最終基于表達(dá)譜數(shù)據(jù)可以識(shí)別出15個(gè)基因模塊，包括:Black 模塊(190 個(gè)基因)、Blue 模塊(576 個(gè)基因)、Brown模塊(430個(gè)基因)、cyan模塊(130個(gè)基因）、Green 模塊(257 個(gè)基因)、Greenyellow 模塊(149 個(gè)基因）、Grey 模塊(109 個(gè))、Magenta 模塊(158 基因）、pink模塊（162個(gè)基因）、purple模塊（150個(gè)基因）、Red模塊(230 個(gè)基因)、samon 模塊(140 個(gè)基因）、tan 模塊(149個(gè)基因)、Turquoise 模塊(937 個(gè)基因)和Yellow 模塊(315個(gè)基因)。鑒于Grey模塊內(nèi)包含的基因是指不能歸類于其他任何模塊的基因，后續(xù)研究中一般將Grey模塊進(jìn)行舍棄。

圖3 構(gòu)建WGCNA分析過程中軟閾值篩選

圖4 基因共表達(dá)網(wǎng)絡(luò)分層聚類樹與共表達(dá)模塊

2.3 與PD 樣本表型的相關(guān)基因模塊識(shí)別及該模塊內(nèi)基因功能和通路富集分析

基于Pearson相關(guān)性分析方法，通過計(jì)算各個(gè)基因模塊與樣本表型之間的關(guān)系，繪制出了基因模塊與樣本表型熱圖，可以看出模塊Turquoise(ME-turquoise)與PD 樣本顯著相關(guān)（r=0.97，P=3e-10)，見圖5、6。隨后對(duì)模塊Turquoise 內(nèi)937 個(gè)基因進(jìn)行功能注釋和通路富集分析，結(jié)果顯示GO 功能注釋BP 條目顯著富集在氯離子跨膜轉(zhuǎn)運(yùn)過程、氧運(yùn)輸、陰離子轉(zhuǎn)運(yùn)和中性粒細(xì)胞生長(zhǎng)及代謝等生物學(xué)過程，見圖7。KEGG 通路富集分析顯示937 個(gè)基因顯著富集在神經(jīng)活動(dòng)配體與受體相關(guān)作用和細(xì)胞因子-細(xì)胞因子受體相互作用等通路途徑。

圖5 WGCNA構(gòu)建模塊（ME值）與樣本信息相關(guān)熱圖

圖6 turquoise 模塊內(nèi)基因MM 與GS 散點(diǎn)圖(紅線表示GS=0.9，MM=0.9)

圖7 turquoise模塊內(nèi)基因GO生物學(xué)過程和KEGG通路分析

2.4 PD表型相關(guān)關(guān)鍵基因識(shí)別

基于上述相關(guān)性分析識(shí)別出與PD顯著相關(guān)基因模塊，參考其他文獻(xiàn)我們將PD 表型最相關(guān)的Turquoise 基因模塊內(nèi)基因基因連接度最大的前10 個(gè)基因定義為PD 關(guān)鍵基因，包括IFITM5、MYPOP、SIK1、NPAS3、SMARCC2、FBXO17、SLC34A3、AMN、TNRC18 和FBRSL1，其中基因IFITM5 和FBRSL1 在Genecards 數(shù)據(jù)庫中已經(jīng)證實(shí)參與PD 進(jìn)程，而剩余8個(gè)基因與PD相關(guān)性的研究尚未見文獻(xiàn)報(bào)道。

2.5 PD表型相關(guān)關(guān)鍵基因差異分析及相關(guān)性分析

利用R 軟件中的“l(fā)imma”包對(duì)上述篩選出的10個(gè)PD 關(guān)鍵基因在正常組織樣本與PD 樣本差異性進(jìn)行分析，結(jié)果顯示，與正常組織相比，篩選出的10 個(gè)PD 關(guān)鍵基因在PD 樣本顯著上升（差異倍數(shù)大于1.2且矯正后P＜0.001）。相關(guān)性分析顯示10 個(gè)PD 關(guān)鍵基因具有強(qiáng)的相關(guān)程度（r＞0.9），見圖8。

圖8 10個(gè)PD相關(guān)關(guān)鍵基因差異分析和共表達(dá)分析

3 討論

PD 是1 種多發(fā)于中老年人、以運(yùn)動(dòng)障礙為主要表現(xiàn)的神經(jīng)系統(tǒng)變性疾病[8]。隨著社會(huì)老齡化的日益加劇，PD 已成為我國(guó)僅次于阿爾茨海默病的神經(jīng)系統(tǒng)性疾病[9]。目前，PD的發(fā)病機(jī)制仍存在諸多未知且臨床上缺乏有效的治療藥物。因此，開展PD 分子機(jī)制的相關(guān)研究對(duì)改善目前PD現(xiàn)狀具有積極意義。

PD 最主要的病理改變比較明確，首先由于中腦黑質(zhì)多巴胺（dopamine,DA）能神經(jīng)元的變性死亡，繼而引起紋狀體DA 含量顯著性減少而導(dǎo)致患者致病[10]。目前關(guān)于PD 臨床特征和發(fā)病機(jī)制并未明了，積極探索PD發(fā)生發(fā)展機(jī)制對(duì)于人們掌握其主要臨床特征，篩選有效的藥物治療靶點(diǎn)具有重要意義。

多基因參與了疾病的發(fā)生發(fā)展過程，PD 也不例外。近年來出現(xiàn)的高通量測(cè)序技術(shù)為人們研究多基因在疾病進(jìn)程中的作用提供了必要條件，相較于傳統(tǒng)的單個(gè)分子的生物學(xué)研究，基于高通量的生物信息學(xué)研究從系統(tǒng)生物學(xué)的角度實(shí)現(xiàn)了對(duì)疾病相關(guān)多個(gè)分子靶點(diǎn)及其功能的探究。WGCNA 分析算法是目前進(jìn)行基因共表達(dá)網(wǎng)絡(luò)研究最常用的生物信息學(xué)分析方法，廣泛用于癌癥研究領(lǐng)域，其穩(wěn)定性和有效性在多項(xiàng)研究中得到證實(shí)，該算法的優(yōu)勢(shì)是在計(jì)算兩兩基因相互關(guān)系的基礎(chǔ)上引入加權(quán)值，從而構(gòu)建出了更符合生物學(xué)意義的基因共表達(dá)網(wǎng)絡(luò)。在本研究中，我們首先將WGCNA 算法用于PD 轉(zhuǎn)錄組基因分析，試圖利用先進(jìn)的生物信息學(xué)分析方法，通過對(duì)PD 高通量芯片數(shù)據(jù)的挖掘，為人們深入探究PD 的發(fā)病機(jī)制提供新的研究思路和必要的候選靶點(diǎn)分子。

通過構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)，我們識(shí)別出Turquoise 模塊內(nèi)基因與PD 具有非常強(qiáng)的相關(guān)性，通過對(duì)模塊內(nèi)基因進(jìn)行功能注釋和通路富集分析，結(jié)果顯示該模塊內(nèi)基因顯著富集在氯離子跨膜轉(zhuǎn)運(yùn)過程、氧運(yùn)輸、陰離子轉(zhuǎn)運(yùn)和中性粒細(xì)胞生長(zhǎng)及代謝相關(guān)生物學(xué)過程和神經(jīng)活動(dòng)配體與受體相關(guān)作用與細(xì)胞因子-細(xì)胞因子受體相互作用等通路途徑，其中已有文獻(xiàn)報(bào)道[11-12]，功能和通路中相關(guān)分子參與了PD 進(jìn)展的過程，本研究結(jié)論與前述文獻(xiàn)報(bào)道較一致。我們隨后對(duì)Turquoise 模塊內(nèi)基因進(jìn)行挖掘，最終識(shí)別出10 個(gè)基因，包括IFITM5、MYPOP、SIK1、NPAS3、SMARCC2、FBXO17、SLC34A3、AMN、TNRC18 和FBRSL1，在Turquoise模塊內(nèi)具有MM值和GS均大于0.9，且網(wǎng)絡(luò)連接度較高的特征。我們將上述10 個(gè)基因定義為PD 相關(guān)的關(guān)鍵基因。與正常樣本相比，10 個(gè)基因在PD 組織中的表達(dá)均顯著升高。進(jìn)一步文獻(xiàn)檢索發(fā)現(xiàn)，此10 個(gè)基因目前在與PD 相關(guān)的研究中均較少出現(xiàn)。利用genecards 數(shù)據(jù)庫檢索此10 個(gè)基因與PD 的相關(guān)性得知，IFITM5 基因和FBRSL1 在genecards 數(shù)據(jù)庫與PD均具有一定的相關(guān)性，其中IFITM5基因即干擾素誘導(dǎo)跨膜蛋白5，作為干擾素誘導(dǎo)跨膜蛋白家族成員之一，在成骨細(xì)胞礦化過程中的具有重要作用。鑒于有學(xué)者進(jìn)行骨髓基質(zhì)細(xì)胞向神經(jīng)細(xì)胞分化及對(duì)PD 的治療研究[13]，我們推測(cè)IFITM5 可能通過影響骨髓基質(zhì)細(xì)胞向成骨細(xì)胞和神經(jīng)細(xì)胞的轉(zhuǎn)化過程環(huán)節(jié)參與了PD 的發(fā)生、發(fā)展過程，盡管在genecards數(shù)據(jù)庫中檢索到FBRSL1 基因?qū)儆赑D 相關(guān)基因，但其具體機(jī)制并不明確。考慮到WGCNA 構(gòu)建同一模塊內(nèi)基因具有相似性的思路，我們對(duì)此10 個(gè)基因進(jìn)行了相關(guān)性分析，結(jié)果顯示，9 個(gè)基因間均與IFITM5具有非常強(qiáng)的相關(guān)性，其參與PD 進(jìn)程的具體機(jī)制值得進(jìn)一步探索。

綜上所述，本研究利用WGCNA分析方法識(shí)別出了與PD具有強(qiáng)相關(guān)性的基因模塊。通過對(duì)該模塊內(nèi)基因進(jìn)行挖掘，識(shí)別出10 個(gè)與PD 強(qiáng)相關(guān)性的關(guān)鍵基因，但目前相關(guān)研究仍非常缺乏。因此，本研究的完成將會(huì)為進(jìn)行PD 致病機(jī)理的相關(guān)研究及PD 藥物靶點(diǎn)研究提供重要的候選靶點(diǎn)分子。