帕金森病診斷基因的生物信息學(xué)分析

陳忠，連世忠

1 山西醫(yī)科大學(xué)第一臨床醫(yī)學(xué)院，太原 030000；2 山西醫(yī)科大學(xué)第一醫(yī)院神經(jīng)外科

帕金森病（PD）是與多巴胺能神經(jīng)元喪失相關(guān)的常見的神經(jīng)退行性疾病之一，其發(fā)病率隨著年齡的增長而穩(wěn)步上升［1］。PD 患者運(yùn)動遲緩并伴有至少一種以上表現(xiàn)（肌肉僵硬、靜止性震顫或姿勢不穩(wěn)）［2］。其主要特征是黑質(zhì)致密部多巴胺能神經(jīng)元早期變性死亡，細(xì)胞內(nèi)廣泛存在α-突觸核蛋白聚集［1，3］。迄今為止，PD 的病理機(jī)制尚未明確，目前研究 較 多 的 是SNCA、LRRK2、VPS35、EIF4G1、DNAJC13、CHCHD2基因突變，其治療仍以左旋多巴為首的藥物對癥治療［2，4］。2021 年10 月，本研究通過生物信息學(xué)分析，對不同數(shù)據(jù)集中的PD基因信息進(jìn)行整合，從而找到PD 的關(guān)鍵基因并建立診斷模型，預(yù)測靶向調(diào)控關(guān)鍵基因的微小RNA（miRNA）和轉(zhuǎn)錄因子（TF），為探索PD 的遺傳學(xué)病因和發(fā)病機(jī)制提供參考。

1 資料與方法

1.1 資料 從美國國家生物信息中心（NCBI）基因表達(dá)綜合數(shù)據(jù)庫（GEO，https：//www. ncbi. nlm. nih.gov/geo/）下載3個(gè)PD數(shù)據(jù)集（GSE20146、GSE20153、GSE20141）作為內(nèi)部訓(xùn)練集，3 個(gè)PD 數(shù)據(jù)集平臺注釋文件為Affymetrix 人基因表達(dá)陣列（GPL570）。GSE20146 數(shù)據(jù)集包括10 例PD 患者樣本和10 例正常樣本；GSE20153包括8例PD患者樣本、8例正常樣本；GSE20141包括10例PD患者樣本、8例正常樣本。此外從NCBI 基因表達(dá)綜合數(shù)據(jù)庫另外下載2 個(gè)PD數(shù)據(jù)集（GSE20291、GSE20292）作為外部驗(yàn)證集，2個(gè)PD 數(shù)據(jù)集平臺注釋文件是Affymetrix 人基因表達(dá)陣列（GPL96）。GSE20291數(shù)據(jù)集包括15例PD 患者的樣本、20例正常樣本；GSE20292包括11例PD患者的樣本、18例正常樣本。

1.2 差異表達(dá)基因（DEG）篩選 用R4.0.2 軟件（https：//www.R-project.org）將3 個(gè)PD 數(shù)據(jù)集合并作為內(nèi)部訓(xùn)練集，進(jìn)行預(yù)處理（背景校正、歸一化、log2 轉(zhuǎn)化）。當(dāng)多個(gè)探針對應(yīng)1 個(gè)共同基因時(shí)，取其平均值作為其表達(dá)值。用“sva”包消除3 個(gè)數(shù)據(jù)集之間的批次效應(yīng)。用“l(fā)imma”包篩選P＜0.05 的DEG，以|log2 Fold change（FC）|＞0.5 作為選擇DEG的截止點(diǎn)（FC為倍性變化）。

1.3 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（WGCNA）及模塊鑒定 為了探討基因間的相互作用，應(yīng)用系統(tǒng)生物學(xué)方法WGCNA 構(gòu)建基因共表達(dá)網(wǎng)絡(luò)。將樣本中有25%以上變異的基因整合的數(shù)據(jù)集導(dǎo)入WGCNA；為保證網(wǎng)絡(luò)構(gòu)建結(jié)果的可靠性，剔除離群樣本；用“pick-Soft-Threshold”函數(shù)，由共同表達(dá)式相似度得到的軟閾值功率β 計(jì)算鄰接度；將鄰接關(guān)系轉(zhuǎn)化為拓?fù)渲丿B矩陣（TOM），并計(jì)算相應(yīng)的不相似度（1-TOM）；通過分層聚類和動態(tài)樹切割函數(shù)對模塊進(jìn)行檢測；為將表達(dá)譜相似的基因分類到基因模塊中，對基因樹形圖進(jìn)行平均連鎖層次聚類，采用“TOMbased”差異測量方法，最小基因組為50；對與臨床屬性相關(guān)的模塊，計(jì)算模塊隸屬度（MM，特定基因與模塊特征基因之間的相關(guān)性）和基因顯著性（GS，特定基因與臨床變量之間的相關(guān)性）。最后，對特征基因網(wǎng)絡(luò)進(jìn)行可視化，進(jìn)一步分析模塊中的基因信息［5］。從內(nèi)部訓(xùn)練集提取的DEG 與WGCNA 重要模塊中的基因取交集得到交集基因。

1.4 關(guān)鍵基因篩選及診斷模型構(gòu)建 逐步回歸法選擇關(guān)鍵基因并采用多因素邏輯回歸分析構(gòu)建診斷模型。對所有樣本中交集基因的表達(dá)量取中位數(shù)，基因表達(dá)量高于中位數(shù)則定義為高表達(dá)，反之則為低表達(dá)。用SPSS23.0 統(tǒng)計(jì)軟件對數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，用逐步回歸篩選出關(guān)鍵基因；用“rms”包采用多因素邏輯回歸分析用關(guān)鍵基因構(gòu)建診斷模型并繪制列線圖。

1.5 診斷模型的校準(zhǔn) 使用“rms”包繪制校準(zhǔn)曲線以評估PD 診斷模型列線圖，用“rms”包測量C 指數(shù)量化列線圖的辨別性能。用R 軟件對2 個(gè)PD 數(shù)據(jù)集（GSE20291、GSE20292）進(jìn)行合并，作為外部驗(yàn)證集，預(yù)處理方法同1.2。用外部驗(yàn)證集驗(yàn)證診斷模型，計(jì)算C 指數(shù)，繪制受試者工作特征（ROC）曲線并計(jì)算曲線下面積。

1.6 miRNA-TF-mRNA 調(diào)控網(wǎng)絡(luò)的構(gòu)建 利用miRTarBase、Starbase 和Targetscan 數(shù)據(jù)庫預(yù)測關(guān)鍵基因的靶向miRNA。為了提高預(yù)測的準(zhǔn)確性，只保留3 個(gè)數(shù)據(jù)庫預(yù)測的miRNA。使用rich 數(shù)據(jù)庫（http：//amp. pharm. mssm. edu/Enrichr/）預(yù)測關(guān)鍵基因的靶向TF。選取P＜0.05 的結(jié)果作為截?cái)嘀怠Ｔ讷@得miRNA-TF-mRNA 調(diào)控關(guān)系后，使用Cytoscape3.7.2 軟件對miRNA-TF-mRNA 調(diào)控網(wǎng)絡(luò)進(jìn)行可視化。

1.7 統(tǒng)計(jì)學(xué)方法 采用SPSS23.0 統(tǒng)計(jì)軟件和R4.0.1 軟件進(jìn)行數(shù)據(jù)處理，用R4.0.1 軟件和Cytoscape3.7.2軟件進(jìn)行圖像生成處理。P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 DEG 篩選結(jié)果 DEG 共405個(gè)，其中表達(dá)上調(diào)基因191個(gè)，表達(dá)下調(diào)基因214個(gè)。表達(dá)上調(diào)基因前5 位為KDM6B、GRINA、IFI35、TRIM36、NME4，表達(dá)下 調(diào) 基 因 前5 位 為BASP1、NCKAP1、CNTNAP2、ZNF536、SYT1。

2.2 WGCNA 分析和模塊鑒定結(jié)果 所有基因按照方差由大到小進(jìn)行排序，選取了方差前25%（3 100 個(gè)）的基因進(jìn)行分析。用“flashClust”工具包進(jìn)行聚類分析，將閾值設(shè)置為100，檢測并刪除3 個(gè)離群樣本，保留了51 個(gè)樣本。用“WGCNA”包中的“pickSoftThreshold”函數(shù)篩選出1~20 的功率參數(shù)。選擇β=12 的冪作為軟閾值，以保證網(wǎng)絡(luò)的無標(biāo)度。將閾值設(shè)置為0.20，合并集群樹中的類似模塊，共獲得8 個(gè)模塊，其中基因具有相似的共表達(dá)性狀。tan 模塊的模塊特征基因（ME）與PD 的正相關(guān)性最高（r=0.34，P＜0.05），red 模塊的ME 與PD 的負(fù)相關(guān)性最高（r=-0.31，P＜0.05）。因此，包含100 個(gè)基因的tan 模塊和包含180 個(gè)基因的red 模塊被鑒定為PD的重要模塊。在tan模塊中獲得基因的MM與GS呈正相關(guān)（cor=0.28，P＜0.05）；在red 模塊中獲得基因的MM 與GS 呈正相關(guān)（cor=0.24，P＜0.05）。從內(nèi)部訓(xùn)練集篩選的DEG 分別與tan 模塊和red 模塊中的基因取交集得到19個(gè)交集基因。

2.3 關(guān)鍵基因及診斷模型 用逐步回歸法從19 個(gè)交集基因中最終篩選出5 個(gè)關(guān)鍵基因，包括PNMA3、AEBP1、PABPC1、GCA、GSTM2。將5 個(gè)關(guān)鍵基因進(jìn)行多因素邏輯回歸分析，見表1，建立包含上述5 個(gè)關(guān)鍵基因的診斷模型，并將其表示為列線圖。該診斷模型的校正曲線顯示一致性良好，計(jì)算C 指數(shù)為0.890，并進(jìn)一步繪制了ROC 曲線，曲線下面積為0.890。

表1 關(guān)鍵基因的多因素邏輯回歸分析

2.4 miRNA-TF-mRNA 調(diào)控網(wǎng)絡(luò)分析結(jié)果 得到24 個(gè)miRNA-TF-mRNA 調(diào)控關(guān)系。該調(diào)控網(wǎng)絡(luò)由Cytoscape3.7.2軟件構(gòu)建，包括8個(gè)miRNA、9個(gè)TF、5 個(gè)mRNA，見 圖1。8 個(gè)miRNA 分 別 為hsa-miR-532-3p、hsa-miR-423-3p、hsa-miR-485-5p、hsa-miR-149-5p、hsa-miR-520g-3p、hsa-miR-429、hsa-miR-200c-3p、hsa-miR-145-5p；9 個(gè)TF 分 別 為SNAI2、ATF4、PLAU、SP3、JUN、NFIA、LTF、TCF3、SNAI1；5個(gè)mRNA 分別為PNMA3、AEBP1、PABPC1、GCA 和GSTM2。

圖1 miRNA-TF-mRNA調(diào)控網(wǎng)絡(luò)

2.5 診斷模型驗(yàn)證結(jié)果 外部驗(yàn)證集驗(yàn)證診斷模型，計(jì)算C 指數(shù)為0.752，并繪制了ROC 曲線，曲線下面積為0.752。

3 討論

PD 是與多巴胺能神經(jīng)元喪失相關(guān)的最常見神經(jīng)退行性疾病之一，其發(fā)生機(jī)制尚未完全明確。在本研究中，我們整合了3 個(gè)PD 數(shù)據(jù)集作為內(nèi)部訓(xùn)練集，使用兩種不同的方法（WGCNA 分析和差異表達(dá)分析）得到交集基因。通過差異表達(dá)分析，我們在PD 樣本和正常對照組之間共鑒定出405 個(gè)DEG，其中表達(dá)上調(diào)基因191 個(gè)，表達(dá)下調(diào)基因214 個(gè)。前5位表達(dá)上調(diào)基因?yàn)镵DM6B、GRINA、IFI35、TRIM36、NME4，而前5 位表達(dá)下調(diào)基因?yàn)锽ASP1、NCKAP1、CNTNAP2、ZNF536、SYT1。有文獻(xiàn)報(bào)道，CNTNAP2可作為參與特發(fā)性和攜帶LRRK2 基因G2019S 突變的PD 患者的潛在候選基因［6］。也有研究發(fā)現(xiàn)，在PD 小鼠模型中miR-34-5p 可以靶向調(diào)控SYT1 從而發(fā)揮神經(jīng)保護(hù)作用［7］。本研究WGCNA 分析發(fā)現(xiàn)，tan模塊的ME與PD的正相關(guān)性最高，red模塊的ME與PD 的負(fù)相關(guān)性最高，因此包含100 個(gè)基因的tan模塊和包含180 個(gè)基因的red 模塊被鑒定為PD 的重要模塊。從內(nèi)部訓(xùn)練集提取的DEG與WGCNA重要模塊中的基因取交集得到了19個(gè)交集基因。

通過對54 個(gè)樣本中19 個(gè)交集基因的表達(dá)量取中位數(shù)，則基因表達(dá)量高于中位數(shù)則定義為高表達(dá)，反之則為低表達(dá)。用逐步回歸選擇關(guān)鍵基因，篩選出5 個(gè)關(guān)鍵基因，包括PNMA3、AEBP1、PABPC1、GCA 和GSTM2，并采用多因素邏輯回歸分析構(gòu)建診斷模型，該診斷模型的校正曲線顯示了良好的一致性，計(jì)算C 指數(shù)為0.890，并進(jìn)一步繪制了ROC 曲線，曲線下面積為0.890，說明了模型的良好分辨力。外部數(shù)據(jù)集也再一次驗(yàn)證了其良好的診斷能力，計(jì)算C 指數(shù)為0.752，曲線下面積為0.752。然而暫無文獻(xiàn)證實(shí)這5 個(gè)基因?qū)D 的診斷有重要意義。但有研究表明，AEBP1 為膠質(zhì)瘤的潛在致癌驅(qū)動基因，對其治療和干預(yù)具有潛在影響［8］。也有文獻(xiàn)報(bào)道，PABPC1 過表達(dá)可抑制膠質(zhì)母細(xì)胞瘤的惡性進(jìn)展［9］，而且GCA 對于檢測銅綠假單胞菌和其他假單胞菌屬物種是否存在有重要意義［10］。研究還表明，GSTM2 是胰腺癌化療優(yōu)化和預(yù)后生物標(biāo)志物的潛在靶標(biāo)［11］。這5 個(gè)基因?qū)τ赑D 的診斷價(jià)值需要進(jìn)一步驗(yàn)證。

通過對關(guān)鍵基因進(jìn)行miRNA-TF-mRNA 調(diào)控網(wǎng)絡(luò)分析得到24 個(gè)miRNA-TF-mRNA 調(diào)控關(guān)系，包括8 個(gè)miRNA 和9 個(gè)TF。研究發(fā)現(xiàn)，抑制miR-429 可以減弱小膠質(zhì)細(xì)胞的炎癥反應(yīng)和創(chuàng)傷性腦損傷介導(dǎo)的腦損傷［12］。LTF 在腦出血后血腫解毒中也有重要作用［13］。NFIA 則被發(fā)現(xiàn)是一種膠質(zhì)的生成開關(guān)，能夠從多能干細(xì)胞中快速衍生出功能性星形膠質(zhì)細(xì)胞［14］。也有研究表明，MALAT1與miR-200c-3p結(jié)合可以上調(diào)SIRT1 表達(dá)，從而誘導(dǎo)腦微血管內(nèi)皮細(xì)胞自噬并保護(hù)腦微血管內(nèi)皮細(xì)胞免受氧和葡萄糖侵害［15］。抑制miR-145-5p 表達(dá)可減小大鼠急性腦缺血中的梗死體積［16］。并且ATF4 的表達(dá)可加深原發(fā)性腦腫瘤的惡性程度，促進(jìn)腫瘤細(xì)胞增殖和腫瘤血管生成［17］。然而以上miRNA 和TF 如何對關(guān)鍵基因進(jìn)行調(diào)控，則需進(jìn)一步驗(yàn)證。

綜上所述，通過綜合生物信息學(xué)分析得到5 個(gè)關(guān) 鍵 基 因，即PNMA3、AEBP1、PABPC1、GCA 和GSTM2，進(jìn)而建立了診斷PD 的模型，并且證實(shí)該模型具有良好的診斷能力。通過對關(guān)鍵基因進(jìn)行miRNA-TF-mRNA 調(diào)控網(wǎng)絡(luò)分析得到24 個(gè)miRNATF-mRNA 調(diào)控關(guān)系，包括8 個(gè)miRNA 和9 個(gè)TF。然而關(guān)鍵基因?qū)τ赑D 的診斷及miRNA 和TF 如何對關(guān)鍵基因進(jìn)行調(diào)控，進(jìn)而影響PD 的進(jìn)展，則需進(jìn)一步研究。以上關(guān)鍵基因、miRNA 和TF 也許會成為PD 的生物標(biāo)志物，同時(shí)也為今后探明PD 的發(fā)病機(jī)制提供了更多方向。