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

陳忠，連世忠

1 山西醫科大學第一臨床醫學院，太原 030000；2 山西醫科大學第一醫院神經外科

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

1 資料與方法

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

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

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

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

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

1.6 miRNA-TF-mRNA 調控網絡的構建 利用miRTarBase、Starbase 和Targetscan 數據庫預測關鍵基因的靶向miRNA。為了提高預測的準確性，只保留3 個數據庫預測的miRNA。使用rich 數據庫（http：//amp. pharm. mssm. edu/Enrichr/）預測關鍵基因的靶向TF。選取P＜0.05 的結果作為截斷值。在獲得miRNA-TF-mRNA 調控關系后，使用Cytoscape3.7.2 軟件對miRNA-TF-mRNA 調控網絡進行可視化。

1.7 統計學方法 采用SPSS23.0 統計軟件和R4.0.1 軟件進行數據處理，用R4.0.1 軟件和Cytoscape3.7.2軟件進行圖像生成處理。P＜0.05為差異有統計學意義。

2 結果

2.1 DEG 篩選結果 DEG 共405個，其中表達上調基因191個，表達下調基因214個。表達上調基因前5 位為KDM6B、GRINA、IFI35、TRIM36、NME4，表達下 調 基 因 前5 位 為BASP1、NCKAP1、CNTNAP2、ZNF536、SYT1。

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

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

表1 關鍵基因的多因素邏輯回歸分析

2.4 miRNA-TF-mRNA 調控網絡分析結果 得到24 個miRNA-TF-mRNA 調控關系。該調控網絡由Cytoscape3.7.2軟件構建，包括8個miRNA、9個TF、5 個mRNA，見 圖1。8 個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 個TF 分 別 為SNAI2、ATF4、PLAU、SP3、JUN、NFIA、LTF、TCF3、SNAI1；5個mRNA 分別為PNMA3、AEBP1、PABPC1、GCA 和GSTM2。

圖1 miRNA-TF-mRNA調控網絡

2.5 診斷模型驗證結果 外部驗證集驗證診斷模型，計算C 指數為0.752，并繪制了ROC 曲線，曲線下面積為0.752。

3 討論

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

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

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

綜上所述，通過綜合生物信息學分析得到5 個關 鍵 基 因，即PNMA3、AEBP1、PABPC1、GCA 和GSTM2，進而建立了診斷PD 的模型，并且證實該模型具有良好的診斷能力。通過對關鍵基因進行miRNA-TF-mRNA 調控網絡分析得到24 個miRNATF-mRNA 調控關系，包括8 個miRNA 和9 個TF。然而關鍵基因對于PD 的診斷及miRNA 和TF 如何對關鍵基因進行調控，進而影響PD 的進展，則需進一步研究。以上關鍵基因、miRNA 和TF 也許會成為PD 的生物標志物，同時也為今后探明PD 的發病機制提供了更多方向。