包成政，王功俊，羅雪蓮，王潔，藍(lán)學(xué)群，王喻，陳金梅，李雪斌,,4

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

本研究基于生物信息學(xué)方法，從GEO數(shù)據(jù)庫(kù)獲取PD小腦總RNA芯片GSE20314、GSE28894 的數(shù)據(jù)，通過(guò)對(duì)該數(shù)據(jù)進(jìn)行l(wèi)imma差異分析和加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA），并通過(guò)VennDiagram 尋找樞紐（Hub）基因。再使用邏輯回歸分析構(gòu)建診斷模型，并通過(guò)2 個(gè)血液樣本數(shù)據(jù)集進(jìn)行驗(yàn)證，計(jì)算受試者工作特征曲線（receiver operator characteristic curve，ROC）的曲線下面積（area under the curve，AUC）評(píng)估其診斷性能。探索小腦在PD中的作用及構(gòu)建具有應(yīng)用價(jià)值的診斷模型。

1 資料與方法

1.1 數(shù)據(jù)采集

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

圖1 GSE20314、GSE28894箱線圖和Density圖

1.2 方法

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

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

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

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

1.3 統(tǒng)計(jì)分析

所有分析在Sangerbox、SPSS 23.0、Graphpad Prism 9.4.0中進(jìn)行，根據(jù)是否符合正態(tài)分布選擇t檢驗(yàn)或秩和檢驗(yàn)；P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 Limma差異基因

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

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

2.2 WGCNA結(jié)果

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

圖3 WGCNA結(jié)果

2.3 Hub基因

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

圖4 Veen圖

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

2.4 診斷模型的構(gòu)建和驗(yàn)證

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

圖6 SBNO1+VPS13C在各組相應(yīng)的ROC曲線和AUC值

3 討論

PD 主要為中腦黑質(zhì)多巴胺能神經(jīng)元缺失和紋狀體多巴胺神經(jīng)遞質(zhì)減少，導(dǎo)致出現(xiàn)運(yùn)動(dòng)遲緩、肌強(qiáng)直等癥狀，治療上常以補(bǔ)充多巴胺為主要治療手段[9]?，F(xiàn)有研究表明PD 震顫產(chǎn)生與小腦有較為明顯的關(guān)系，PD患者腹側(cè)被蓋區(qū)小腦連接性高于黑質(zhì)小腦連接性，在正常人群中這種現(xiàn)象則會(huì)相反[10]。Kawabata等[11]發(fā)現(xiàn)小腦和基底神經(jīng)節(jié)間的功能連接不僅與運(yùn)動(dòng)相關(guān)，與認(rèn)知表現(xiàn)也有聯(lián)系。小腦的功能或形態(tài)調(diào)節(jié)不僅與運(yùn)動(dòng)癥狀相關(guān)，與一些非運(yùn)動(dòng)癥狀也有不可忽視的聯(lián)系，在PD 中小腦發(fā)揮著代償?shù)淖饔肹12]。早期PD 患者中，小腦-丘腦-皮質(zhì)回路的募集與運(yùn)動(dòng)特征的進(jìn)展呈正相關(guān)，同樣也參與PD的補(bǔ)償機(jī)制[13]。甚至有學(xué)者提出小腦可能參與PD 的猝死機(jī)制[14]。Rusholt 等[15]的體視學(xué)研究，揭示了PD的白質(zhì)改變情況。基于現(xiàn)有小腦參與PD 發(fā)病機(jī)制的學(xué)說(shuō)，相關(guān)藥物的開發(fā)也逐步開展，如姜黃素膳食補(bǔ)充劑可能對(duì)PD 的小腦相關(guān)癥狀具有保護(hù)作用[16]。越來(lái)越多證據(jù)表明小腦在PD中的重要性。

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

鑒于PD 與小腦間存在著相應(yīng)聯(lián)系，我們認(rèn)為SBNO1、VPS13C等Hub基因有可能參與PD的發(fā)生、發(fā)展過(guò)程，未來(lái)或可作為PD的潛在干預(yù)靶點(diǎn)，且SBNO1+VPS13C 基因構(gòu)建的診斷模型具有較好的診斷性能。但其中的相關(guān)機(jī)制及應(yīng)用價(jià)值仍需進(jìn)一步探討及研究。