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

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

1 材料和方法

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

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

1.2 WGCNA分析

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

1.3 關鍵基因的篩選

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

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

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

1.5 PD關鍵基因差異分析與共表達分析

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

1.6 PD相關基因的篩選

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

2 結果

2.1 芯片表達譜預處理及加權基因共表達網(wǎng)絡的構建

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

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

圖2 樣本聚類與對應的樣本信息

2.2 加權基因共表達網(wǎng)絡的構建

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

圖3 構建WGCNA分析過程中軟閾值篩選

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

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

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

圖5 WGCNA構建模塊（ME值）與樣本信息相關熱圖

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

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

2.4 PD表型相關關鍵基因識別

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

2.5 PD表型相關關鍵基因差異分析及相關性分析

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

圖8 10個PD相關關鍵基因差異分析和共表達分析

3 討論

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

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

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

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

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