帕金森病小腦樞紐基因及診斷模型

包成政，王功俊，羅雪蓮，王潔，藍學群，王喻，陳金梅，李雪斌,,4

帕金森?。≒arkinson’s disease，PD）是第二大退行性病變，臨床主要表現為靜止性震顫、運動遲緩、肌強直、姿勢步態異常等癥狀[1]。隨著疾病進展，患者可因治療效果差、并發癥多等導致生活質量嚴重下降，甚至死亡。深入研究PD 的發病機制，尋找更多早期診斷標志物及潛在干預靶點使PD得以早診斷、早治療是現階段的熱門研究方向。目前PD的研究主要集中在基底節區。有學者提出小腦有可能參與PD震顫的發生或發展過程[2]，也有學者認為小腦-丘腦-皮質回路有可能參與PD 震顫[3]。小腦的功能異常有可能是PD凍結步態的基礎[4]，使用神經調節技術可改善PD 患者運動癥狀[5]，刺激小腦可以降低PD 患者左旋多巴治療誘發的運動障礙的嚴重程度，小腦可能參與PD的疾病進程[6]。小腦在PD 的發病機制中具有重要地位，但目前相關研究較少，同時探索PD的診斷方式也是現行研究的熱點。

本研究基于生物信息學方法，從GEO數據庫獲取PD小腦總RNA芯片GSE20314、GSE28894 的數據，通過對該數據進行limma差異分析和加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA），并通過VennDiagram 尋找樞紐（Hub）基因。再使用邏輯回歸分析構建診斷模型，并通過2 個血液樣本數據集進行驗證，計算受試者工作特征曲線（receiver operator characteristic curve，ROC）的曲線下面積（area under the curve，AUC）評估其診斷性能。探索小腦在PD中的作用及構建具有應用價值的診斷模型。

1 資料與方法

1.1 數據采集

數據集來自NCBI 基因表達綜合公共數據庫（GEO）。GSE20314 和GSE28894 是在小腦中表達的RNA數據集，GPL96注釋的GSE20314包括4個對照組人群和4 個PD 患者的小腦RNA 注釋樣本；GPL6104注釋的GSE28894 包括14 個對照組人群和15 個PD 患者的小腦RNA注釋樣本。在分析前將2樣本去批次及歸一化，見圖1。GSE18838和GSE6613是外周血RNA表達的2個數據集，其中GPL5175注釋的GSE18838包括11個正常人樣本和17個PD患者樣本，GPL96注釋的GSE6613包括22個正常人樣本和50個PD患者樣本。

圖1 GSE20314、GSE28894箱線圖和Density圖

1.2 方法

1.2.1 差異表達分析 Limma是一種基于廣義線性模型的差異表達篩選方法[7]。在分析前將GSE20314 和GSE28894 利用“inSilicoMerging”包進行合并，再使用經典貝葉斯方法進行去除批次效應[8]。我們使用R 軟件limma包（version 3.40.6）進行差異分析。P＜0.05和差異倍數1.1 倍的基因被視為差異基因（DEGs）。DEGs的篩選及熱圖和火山圖繪制使用Sangerbox在線分析工具（http://vip.sangerbox.com/home.html）。

1.2.2 加權基因共表達網絡分析 基因之間共表達及基因與表型之間的關系探索使用在線分析工具Sangerbox（http://vip.sangerbox.com/home.html）。首先分別計算每個基因的MAD（median absolute deviation），剔除MAD 最小的前50%的基因。利用R 軟件包WGCNA的goodSamplesGenes方法去除離群的基因和樣本，進一步構建scale-free co-expression network。首先對所有配對基因分別采用Pearson相關矩陣和平均連鎖法，然后用冪函數A_mn=|C_mn|^β（C_mn=基因m和基因n之間的Pearson相關關系，A_mn=基因m和基因n 之間的鄰接關系）構造加權鄰接矩陣（β是一個能強調基因之間強相關性和懲罰弱相關性的軟閾值參數）。選擇適當的冪后，將鄰接關系轉化為一個拓撲重疊矩陣（TOM），該矩陣可以度量一個基因的網絡連通性，定義為該基因與所有其他基因的鄰接關系之和作為網絡基因比，并計算相應的不相似度（1-TOM）。為了將具有相似表達譜的基因劃分為基因模塊，根據基于tom的不相似度度量進行平均連鎖層次聚類。為進一步計算模塊特征基因的不同相似度，選擇一條切線作為模塊樹狀圖，合并一些模塊，并對每個模塊賦予不同的顏色描繪。模塊隸屬度（MM）是基因表達值與模塊ME 的相關系數，表示該基因與該模塊的相關性?；蝻@著性（GS）是基因表達值與表型之間的相關系數，代表基因與表型之間的關聯。

1.2.3 樞紐（Hub）基因鑒定 使用“VennDiagram”包獲取DEGs 和WGCNA 獲得的關鍵模塊基因的交集，進而獲得與PD較為密切相關的Hub基因，Hub基因在對照組人群和PD患者小腦中的表達使用箱線圖表示。

1.2.4 邏輯回歸模型 邏輯回歸是一種廣義線性回歸分析模型，可用于疾病的自動診斷。本研究采用2 個響應變量的邏輯回歸，響應變量為1 時表示PD 樣本，響應變量為0時表示正常樣本。逐步回歸分析最初用于消除對響應變量不顯著的因素，僅保留顯著的因素以簡化模型。逐步回歸從模型中迭代地添加或刪除變量，直到保留最顯著因素。之后使用邏輯回歸來擬合這些顯著因素與響應變量之間的關系。最后使用受試者工作特征曲線（ROC）和ROC 曲線下面積（AUC）評估模型的診斷效果。

1.3 統計分析

所有分析在Sangerbox、SPSS 23.0、Graphpad Prism 9.4.0中進行，根據是否符合正態分布選擇t檢驗或秩和檢驗；P＜0.05為差異有統計學意義。

2 結果

2.1 Limma差異基因

DEGs以P＜0.05、差異倍數1.1倍的基因作為篩選條件，共得到116 個DEGs，使用火山圖及熱圖（顯示t值前10位及最后10位基因）表示，見圖2。

圖2 Limma差異分析火山圖及熱圖

2.2 WGCNA結果

從2個數據集合集選取基因并用于構建加權基因共表達網絡。當軟閾值功率設置為16時，尺度獨立性達到0.89，平均連接值為0.29。當最小模塊大小設置為30時，合并距離＜0.25的模塊，最終獲得13個共表達模塊；值得注意的是grey模塊被認為無法被分配給任何模塊的基因集合。然后對每個模塊與臨床特征進行相關性分析。purple 模塊與PD 的正相關性最高（r=0.35，P=0.03）。計算與基因的表達相關性獲得GS，同時計算模塊特征向量與基因的表達相關性獲得MM，基于截止標準（|MM|＞0.6及|GS|＞0.1），purple模塊中具有高連通性的188個重要基因被提取出來用于進一步分析，見圖3。

圖3 WGCNA結果

2.3 Hub基因

使用VennDiagram 獲得到7 個Hub 基因，使用Veen圖表示，見圖4；使用箱線圖表示各基因在各組的表達，見圖5。

圖4 Veen圖

圖5 GSE20314和GSE28894小腦組織中的Hub基因

2.4 診斷模型的構建和驗證

使用邏輯回歸算法構建基因預測模型，逐步回歸分析提示SBNO1、VPS13C 對PD 的診斷影響最顯著。在小腦樣本數據集GSE20314、GSE28894 中，SBNO1+VPS13C ROC曲線AUC值為0.9386，見圖6A。在血液樣本數據集GSE18838、GSE6613中進一步驗證該模型提示具有良好的診斷性能，AUC 值分別為0.8663、0.6564，見圖6B、C。

圖6 SBNO1+VPS13C在各組相應的ROC曲線和AUC值

3 討論

PD 主要為中腦黑質多巴胺能神經元缺失和紋狀體多巴胺神經遞質減少，導致出現運動遲緩、肌強直等癥狀，治療上常以補充多巴胺為主要治療手段[9]?，F有研究表明PD 震顫產生與小腦有較為明顯的關系，PD患者腹側被蓋區小腦連接性高于黑質小腦連接性，在正常人群中這種現象則會相反[10]。Kawabata等[11]發現小腦和基底神經節間的功能連接不僅與運動相關，與認知表現也有聯系。小腦的功能或形態調節不僅與運動癥狀相關，與一些非運動癥狀也有不可忽視的聯系，在PD 中小腦發揮著代償的作用[12]。早期PD 患者中，小腦-丘腦-皮質回路的募集與運動特征的進展呈正相關，同樣也參與PD的補償機制[13]。甚至有學者提出小腦可能參與PD 的猝死機制[14]。Rusholt 等[15]的體視學研究，揭示了PD的白質改變情況?；诂F有小腦參與PD 發病機制的學說，相關藥物的開發也逐步開展，如姜黃素膳食補充劑可能對PD 的小腦相關癥狀具有保護作用[16]。越來越多證據表明小腦在PD中的重要性。

人類VPS13基因的突變是造成神經發育和神經退行性疾病的原因，VPS13可能在細胞器之間的脂質交換中起作用。作為VPS13之一的液泡蛋白13C（VPS13C）與PD 有著較為密切的聯系，而VPS13A 與舞蹈癥棘皮細胞增多癥存在聯系[17]。通過全外顯子組測序已成功識別導致家族性PD 的基因，發現VPS13C 位列其中[18]。已有研究提示VPS13C 變異與PD 相關[19]。PD 和路易體癡呆均為路易體病，VPS13C 與兩者都有聯系。VPS13C中罕見的錯義突變與路易體病相關，隱性復合雜合錯義突變可能對VPS13C的表達和功能產生不同的影響。罕見的隱性復合雜合錯義突變的組合可降低VPS13C 表達并導致路易體病風險增加[20]。由此可見VPS13C與PD間具有密切關系。在PD細胞模型中，可見VPS13C 部分定位于線粒體的外膜。VPS13C 沉默與較低的線粒體膜電位、線粒體片段化、呼吸頻率增加、PINK1/Parkin 依賴的線粒體惡化以及響應線粒體損傷的PARK2轉錄上調有關。VPS13C功能喪失會增加線粒體對壓力的脆弱性，從而激活PINK1/Parkin 依賴性線粒體質量控制途徑[21]，或許由此參與PD 的發病過程。不同PD 患者群之間VPS13C存在著些許區別。VPS13C 基因的單核苷酸多態性rs2414739 與等位基因模型的PD 易感性存在顯著聯系，但這不包括中國人群[22,23]。但在另外一項研究中篩選中國早發性PD患者隊列中的VPS13C突變，結果發現PD 患者中仍有攜帶VPS13C 的復合雜合突變[24]。VPS13C也與東亞人群有著顯著聯系[25]。在伊朗人中，VPS13C rs2414739 在PD 與正常人群中也存在顯著差異[26]。VPS13C在不同人群中是否有著區別，需進一步擴大樣本去驗證。綜上研究可見VPS13C 與PD 間密切的關系，這與我們的研究結果相符。我們的研究結果進一步發現VPS13C 與SBNO1 基因共同構建基于小腦的PD 診斷模型具有良好的診斷性能。在外周血液中，其診斷性能較腦組織有所下降，這或許與腦組織樣本往往比血液樣本更能代替PD 的病理改變有關。綜合本研究結果及現有研究報道，我們推測該模型對PD患者的診斷或許存在一定程度上的指導意義。

鑒于PD 與小腦間存在著相應聯系，我們認為SBNO1、VPS13C等Hub基因有可能參與PD的發生、發展過程，未來或可作為PD的潛在干預靶點，且SBNO1+VPS13C 基因構建的診斷模型具有較好的診斷性能。但其中的相關機制及應用價值仍需進一步探討及研究。