劉丹 張敏 和申 彭代輝 肖澤萍

重性抑郁障礙（major depressive disorder，MDD）是常見的精神疾病，嚴(yán)重影響患者的生命健康[1]，但目前仍無單一機(jī)制能夠全面解釋該疾病的發(fā)生[2]。該疾病的診斷主要依靠癥狀學(xué)標(biāo)準(zhǔn)，尚無明確的用于診斷的生物標(biāo)志物存在。因此，探索MDD 生物標(biāo)志物至關(guān)重要。穩(wěn)健排序聚合（robust rank aggregation，RRA）算法是一種基于次序統(tǒng)計(jì)整合排序列表的算法[3]。由于其對(duì)高通量數(shù)據(jù)的異常值及噪音具有魯棒穩(wěn)健性，且不依賴于整合不同數(shù)據(jù)集的基因表達(dá)水平，而是整合各個(gè)數(shù)據(jù)集的排序列表，尤其適用于對(duì)不同平臺(tái)的測(cè)序數(shù)據(jù)進(jìn)行meta 分析[3]。本研究擬首次使用RRA 算法對(duì)公共數(shù)據(jù)庫中MDD 外周血轉(zhuǎn)錄組芯片數(shù)據(jù)進(jìn)行整合分析，篩選關(guān)鍵基因并建立診斷模型，為探索MDD 的生物標(biāo)志物提供參考。

1 資料與方法

1.1 數(shù)據(jù)來源從美國國家生物信息中心（NCBI）GEO數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo/）中檢索人類MDD 外周血轉(zhuǎn)錄組芯片數(shù)據(jù)集。納入標(biāo)準(zhǔn)：①研究設(shè)計(jì)為病例對(duì)照研究，病例組患者診斷須為MDD；②患者除MDD外不共患其他疾??；③數(shù)據(jù)集注釋后至少包含15000 個(gè)基因。篩選后共得到5 個(gè)MDD 外周血轉(zhuǎn)錄組芯片數(shù)據(jù)集，分別為GSE98793[4]、GSE76826[5]、GSE52790[6]、GSE38206[7]、GSE39653[8]。其中GSE98793 包含128 例MDD 患者及64 名健康對(duì)照樣本；GSE76826 包含10 例MDD患者及12 名健康對(duì)照樣本；GSE52790 包含10 例MDD 患者及12 名健康對(duì)照樣本；GSE38206 包含9例MDD 患者及9 名健康對(duì)照樣本；GSE39653 包含21例MDD患者及24名健康對(duì)照樣本。

1.2RRA 差異表達(dá)基因（DEG）篩選RRA 算法通過對(duì)比每一個(gè)項(xiàng)目在列表中的實(shí)際位置及零假設(shè)隨機(jī)排序的位置，對(duì)該項(xiàng)目進(jìn)行顯著性評(píng)分，用來重新排序并確定其顯著性。當(dāng)將該算法用于差異表達(dá)基因篩選時(shí)，其可通過對(duì)比每一個(gè)基因在每個(gè)數(shù)據(jù)集差異表達(dá)基因列表中的實(shí)際位置及零假設(shè)隨機(jī)排序的位置，對(duì)該基因進(jìn)行顯著性評(píng)分，確定其差異表達(dá)的統(tǒng)計(jì)學(xué)意義。使用R4.2.0中l(wèi)imma數(shù)據(jù)分析包分別對(duì)經(jīng)批次校正后的5 個(gè)數(shù)據(jù)集進(jìn)行差異表達(dá)基因篩選，設(shè)置|log2 fold change（FC）|＞0.5和P＜0.05 為篩選閾值，得到每個(gè)數(shù)據(jù)集分別的差異表達(dá)基因列表。使用RRA 包對(duì)5 個(gè)數(shù)據(jù)集的差異表達(dá)基因結(jié)果進(jìn)行整合分析，得出RRA 整合分析后的差異表達(dá)基因（P＜0.05）。

1.3 差異表達(dá)基因功能與信號(hào)通路富集分析使用在線工具M(jìn)etascape 對(duì)差異表達(dá)基因進(jìn)行基因本體（gene ontology，GO）功能信號(hào)通路富集分析，隨后使用R軟件進(jìn)行可視化分析。檢驗(yàn)水準(zhǔn)α=0.05。

1.4 關(guān)鍵基因篩選及診斷模型構(gòu)建Boruta 法篩選關(guān)鍵基因并采用多因素logistic 回歸分析構(gòu)建診斷模型。以樣本量最大的GSE98793 為訓(xùn)練集，使用R 中Boruta 包從差異表達(dá)基因中進(jìn)行關(guān)鍵基因篩選。Boruta 法是針對(duì)分類問題非常有效的特征選擇方法，其依托于隨機(jī)森林算法，能有效消除不重要的特征[9]。在篩選特征時(shí)，Boruta 法根據(jù)特征的重要性將特征分為“確認(rèn)”、“待定”及“拒絕”三類，選取前兩類基因作為關(guān)鍵基因，并使用rms 包采用多因素logistic 回歸構(gòu)建診斷模型。使用forestplot包繪制關(guān)鍵基因森林圖，分析單個(gè)關(guān)鍵基因的表達(dá)水平與MDD 的關(guān)系。使用roc.test法比較單個(gè)關(guān)鍵基因及診斷模型在訓(xùn)練集中的診斷價(jià)值。

1.5 診斷模型診斷價(jià)值分析采用Bootstrap 重抽樣法對(duì)GSE98793 重抽樣1000 次作為內(nèi)部驗(yàn)證。研究表明，相比于數(shù)據(jù)集拆分等其他方法，Bootstrap法在針對(duì)logistic 回歸模型的內(nèi)部驗(yàn)證方面具有更好的穩(wěn)定性和更小的偏倚[10]。將剩余4 個(gè)數(shù)據(jù)集（GSE76826、GSE52790、GSE38206、GSE39653）作為外部驗(yàn)證集，繪制診斷模型區(qū)分MDD 患者和健康對(duì)照的受試者工作特征（receiver operating characteristic,ROC）曲線并計(jì)算曲線下面積（area under curve,AUC），對(duì)診斷模型的診斷價(jià)值進(jìn)行外部驗(yàn)證。

本研究中數(shù)據(jù)分析處理及圖像生成均采用R4.2.0。檢驗(yàn)水準(zhǔn)α=0.05。

2 結(jié)果

2.1 差異表達(dá)基因分析結(jié)果使用RRA 聯(lián)合分析后共篩選出MDD 與健康對(duì)照差異表達(dá)基因31 個(gè)，其中上調(diào)基因20 個(gè)，下調(diào)基因11 個(gè)。差異表達(dá)基因熱圖見圖1。

圖1 MDD 與健康對(duì)照差異基因logFC 熱圖

2.2 差異表達(dá)基因的GO功能信號(hào)通路分析GO功能分析結(jié)果顯示，差異表達(dá)基因主要富集的生物過程包括白細(xì)胞激活調(diào)控、損傷反應(yīng)、刺激檢測(cè)、細(xì)胞因子生成負(fù)調(diào)控，細(xì)胞組分主要包括突觸前、特定顆粒腔及裂解空泡?？梢暬Y(jié)果見圖2。

圖2 差異基因GO 富集分析

2.3 關(guān)鍵基因篩選及診斷模型構(gòu)建使用GSE98793 作為模型訓(xùn)練集，使用Boruta 法對(duì)31 個(gè)差異表達(dá)基因進(jìn)行篩選。Boruta 法共篩選出3 個(gè)“確認(rèn)”基因（MMP8、TDRD9、FAM3B），4 個(gè)“待定”基因（LCN2、ARG1、NPTN、FANCF）。將上述7 個(gè)關(guān)鍵基因進(jìn)行多因素logistic 回歸分析并建立MDD 疾病診斷模型。關(guān)鍵基因森林圖表明：FANCF（OR=0.277,95%CI：0.100～0.744）及FAM3B（OR=0.763,95%CI：0.602～0.967）為MDD 的保護(hù)因素，其表達(dá)水平升高可以降低MDD 的發(fā)病風(fēng)險(xiǎn)；TDRD9（OR=2.362,95%CI：1.434～4.071）為MDD 的危險(xiǎn)因素，其表達(dá)水平升高可增大MDD 的發(fā)病風(fēng)險(xiǎn)（圖3）。繪制該模型在訓(xùn)練集中的ROC 曲線（圖4），AUC 為0.803（95%CI：0.740～0.867）。同時(shí)，各關(guān)鍵基因的ROC 曲線及AUC 值如圖4。其中MMP8的AUC 值最大，為0.702（95%CI：0.625～0.779）。使用roc.test法比較模型及MMP8的AUC 值，結(jié)果顯示模型的AUC 值高于MMP8的AUC 值，差異具有統(tǒng)計(jì)學(xué)意義（Z=2.793,P＜0.001）。

圖3 關(guān)鍵基因森林圖

2.4 診斷模型內(nèi)部及外部驗(yàn)證內(nèi)部驗(yàn)證結(jié)果顯示AUC 為0.804（95%CI：0.757～0.851）。內(nèi)部驗(yàn)證的ROC 曲線見圖5，模型的校準(zhǔn)曲線見圖6。同時(shí)，該模型在外部驗(yàn)證集中表現(xiàn)出較好的診斷性能。在GSE76826 中的AUC 為0.781（圖7A），在GSE38206中的AUC 為0.901（圖7B），在GSE39653 中的AUC為0.722（圖7C），在GSE52790 中的AUC 為0.725（圖7D）。

圖5 模型內(nèi)部驗(yàn)證的ROC 曲線

圖6 模型內(nèi)部驗(yàn)證校正曲線圖

圖7 模型在各外部驗(yàn)證集中的ROC曲線

3 討論

本研究使用RRA 的方法整合5 個(gè)MDD 外周血轉(zhuǎn)錄組數(shù)據(jù)集，得到31 個(gè)差異表達(dá)基因。將包含最大樣本量共192 例樣本的數(shù)據(jù)集GSE98793 作為訓(xùn)練集，從31個(gè)差異表達(dá)基因中通過Boruta法篩選出7 個(gè)基因作為關(guān)鍵基因建立診斷模型，該模型在訓(xùn)練集及為外部驗(yàn)證集中均表現(xiàn)出具有較好的預(yù)測(cè)能力。

既往基于轉(zhuǎn)錄組數(shù)據(jù)篩選MDD 生物標(biāo)志物及構(gòu)建診斷模型的研究存在一定局限，如HE 等[11]使用GSE98793 數(shù)據(jù)集建立了4 基因診斷模型，其AUC 為0.799，但其數(shù)據(jù)來源單一，僅使用了兩個(gè)轉(zhuǎn)錄組數(shù)據(jù)集，未驗(yàn)證其模型在其他來源數(shù)據(jù)上的診斷價(jià)值。在轉(zhuǎn)錄組芯片數(shù)據(jù)集的meta分析方面，目前僅有1 項(xiàng)研究整合MDD 相關(guān)數(shù)據(jù)集構(gòu)建了包含70個(gè)關(guān)鍵基因的診斷模型，該模型在外部驗(yàn)證集中AUC 為0.78[12]。而本研究構(gòu)建的診斷模型僅包含7個(gè)基因，同時(shí)在外部驗(yàn)證集上表現(xiàn)較好，因此在可行性上明顯優(yōu)于上述模型。

對(duì)RRA 整合分析后得到的差異表達(dá)基因進(jìn)行GO 功能分析，結(jié)果顯示差異表達(dá)基因主要富集的生物過程包括白細(xì)胞激活調(diào)控、損傷反應(yīng)、刺激檢測(cè)及細(xì)胞因子生成負(fù)調(diào)控，這些通路均與免疫及免疫炎癥相關(guān)。目前已有大量文獻(xiàn)表明，免疫失調(diào)及免疫炎癥反應(yīng)與MDD 的病理生理機(jī)制及預(yù)后相關(guān)[13-14]。差異表達(dá)基因主要富集的GO 細(xì)胞組分包括突觸前、特定顆粒腔及裂解空泡，其均與神經(jīng)遞質(zhì)傳遞及神經(jīng)元活動(dòng)相關(guān)。

本研究診斷模型納入MMP8、TDRD9、FAM3B、LCN2、ARG1、NPTN、FANCF等7 個(gè)基因。其中，NPTN基因編碼兩種亞型的跨膜糖蛋白，Np55 及Np65。Np65 是僅表達(dá)于大腦中的細(xì)胞黏附分子，主要表達(dá)于杏仁核、海馬等與記憶及情緒調(diào)節(jié)相關(guān)的腦區(qū)中。在分子層面，其能與γ-氨基丁酸（γ-aminobutyric acid，GABA）受體結(jié)合，影響GABA能神經(jīng)元的信息傳遞[15]。GABA 能系統(tǒng)參與包括抑郁癥在內(nèi)的多種神經(jīng)精神疾病的發(fā)生[16-17]，作為新型抗抑郁藥物的可能靶標(biāo)之一，相關(guān)研究表明提升GABA 神經(jīng)遞質(zhì)水平在抑郁癥模型小鼠中具有抗抑郁效果[18]。同時(shí)，動(dòng)物實(shí)驗(yàn)表明NPTN基因失活影響5-HT水平，導(dǎo)致小鼠皮質(zhì)酮水平增高，增加小鼠抑郁樣行為[19]。LCN2編碼Lipocalin-2 蛋白，面對(duì)炎癥刺激時(shí)誘導(dǎo)中樞神經(jīng)系統(tǒng)產(chǎn)生趨化因子，參與先天性免疫反應(yīng)，調(diào)節(jié)神經(jīng)元炎癥及神經(jīng)元降解[20]。同時(shí)，LCN2還參與調(diào)節(jié)多種生物行為反應(yīng)，如認(rèn)知功能、疼痛敏感性、抑郁和焦慮[20-22]。研究表明，LCN2敲除小鼠出現(xiàn)焦慮及抑郁樣行為，另有研究表明出現(xiàn)海馬腹側(cè)顆粒細(xì)胞和錐體細(xì)胞肥大，而這兩種細(xì)胞的功能與情緒調(diào)節(jié)相關(guān)[22]。一項(xiàng)以老年MDD 患者為研究對(duì)象的研究發(fā)現(xiàn)老年MDD 患者的血漿Lipocalin-2 蛋白水平較健康對(duì)照顯著上升[23]，另有研究表明血漿Lipocalin-2 蛋白水平與心衰患者的抑郁癥狀嚴(yán)重程度相關(guān)[24-25]。動(dòng)物研究表明，LCN2參與調(diào)節(jié)炎癥性腸病小鼠模型的抑郁癥狀[26]。兩項(xiàng)全基因組關(guān)聯(lián)分析研究表明MMP8與MDD 及圍產(chǎn)期抑郁障礙的發(fā)生相關(guān)[27-28]。2021 年的一項(xiàng)研究表明在慢性輕度應(yīng)激小鼠模型的海馬中過表達(dá)白細(xì)胞介素4 可以上調(diào)ARG1表達(dá)，進(jìn)而促進(jìn)海馬神經(jīng)元生成并降低應(yīng)激相關(guān)抑郁的風(fēng)險(xiǎn)[29]。MMP8、FAM3B及ARG1均參與免疫炎癥反應(yīng)及細(xì)胞因子調(diào)節(jié)，而如前文所述，免疫炎癥反應(yīng)與MDD 的發(fā)病機(jī)制密切相關(guān)[13,30]，因此，這3 個(gè)基因也可能通過調(diào)控免疫炎癥反應(yīng)介導(dǎo)抑郁的發(fā)生。目前尚未有研究探索FANCF或TDRD9與MDD 之間的關(guān)系，而由于本研究為基于生物信息學(xué)的研究，無法探究其作用機(jī)制，故仍需進(jìn)一步的基礎(chǔ)實(shí)驗(yàn)探索可能的作用機(jī)制。

綜上所述，本研究通過對(duì)現(xiàn)有MDD 公共數(shù)據(jù)庫數(shù)據(jù)進(jìn)行整合分析，篩選出7 個(gè)MDD 關(guān)鍵基因，并構(gòu)建出在內(nèi)部及外部驗(yàn)證集中均表現(xiàn)出較好診斷能力的診斷模型。同時(shí)，7 個(gè)關(guān)鍵基因中2 個(gè)基因與MDD 的關(guān)系目前尚未見研究報(bào)告。因此，本研究為進(jìn)一步探究MDD 的發(fā)病機(jī)制提供了新的可能方向，為進(jìn)一步明確其潛在病理機(jī)制及基于生物標(biāo)志物的診斷提供了依據(jù)。