類風(fēng)濕關(guān)節(jié)炎免疫特征基因及其與土茯苓黃酮類化合物的關(guān)系研究

楊 欣,黃 聰,姚血明

類風(fēng)濕關(guān)節(jié)炎(rheumatoid arthritis,RA)是一種慢性的、高度致殘的自身免疫性疾病。RA患者關(guān)節(jié)滑膜液中存在大量的炎性細胞浸潤,如巨噬細胞、樹突狀細胞、T細胞等[1]。土茯苓是百合科植物光葉菝葜的干燥根莖,其在抗炎、免疫系統(tǒng)和腫瘤等方面具有明顯的藥理活性。已有研究[2]表明土茯苓具有免疫調(diào)節(jié)作用,其主要成分落新婦苷具有抗痛風(fēng)性關(guān)節(jié)炎的作用。目前,生物信息學(xué)技術(shù)與高通量測序快速發(fā)展,能快速獲取疾病相關(guān)的生物分子信息,并對相關(guān)的生物分子進行功能研究。同時,基因共表達網(wǎng)絡(luò)分析(weighted correlation network analysis,WGCNA)在尋找RA潛在生物標(biāo)志物、疾病治療特征基因等中發(fā)揮重要作用[3]。該研究利用了GEO數(shù)據(jù)集結(jié)合最小絕對收縮和選擇算法(least absolute shrinkage and selection operator,LASSO)、單樣本基因集富集分析(single sample gene set enrichment analysis,ssGSEA)和 WGCNA 等篩選與特定免疫相關(guān)的RA診斷標(biāo)志物。并通過超高效液相色譜-四極桿-靜電場軌道阱高分辨質(zhì)譜 (UHPLC-Q-Exactive Orbitrap MS)技術(shù)鑒定土茯苓黃酮類化合物,分子對接篩選土茯苓黃酮類化合物與靶點的結(jié)合情況,為本課題組后面的實驗研究以及未來相關(guān)的科學(xué)研究提供依據(jù)。

1 材料與方法

1.1 材料

1.1.1藥材 2022年4月土茯苓藥材采自貴州省安順市,由貴州中醫(yī)藥大學(xué)魏升華教授鑒定為百合科植物光葉菝葜( Smilax glabra Roxb.) 的干燥根莖,植物標(biāo)本存放在貴州中醫(yī)藥大學(xué)棟垣樓417(標(biāo)本號:TFL2022.04)。

1.1.2主要試劑 甲醇、乙腈和甲酸(貨號:34860-1L-R、34851-1L、F0507-1L,美國sigma公司),L-2-氯苯丙氨酸(貨號:C2001,上海恒柏生物科技有限公司)。

1.1.3主要儀器 離心機(型號:Heraeus Fresco 17)、高分辨質(zhì)譜儀(型號:Q Exactive Focus)和超高效液相色譜(型號:Vanquish)均購自賽默飛世爾科技中國有限公司,電子天平(型號:BSA124S-CW,德國賽多利斯科學(xué)儀器有限公司),研磨儀(型號:JXFSTPRP-24,上海凈信科技有限公司),超聲儀(型號:YM-080S,深圳市方奧微電子有限公司),色譜柱(型號:ACQUITY UPLC BEH C18,沃特世科技上海有限公司)。

1.2 方法

1.2.1基因芯片處理 借助美國國立生物技術(shù)信息中心的Gene Expression Omnibus(GEO)數(shù)據(jù)庫檢索 “rheumatoid arthritis” 樣本。將R 4.3.0軟件的“surrogate variable analysis (SVA)”軟件包用于消除批間差異后使用二維主成分分析(principal component analysis,PCA)群集圖顯示樣本間校正的效果[4]。正常對照組和RA組差異基因識別采用(Linear Models for Microarray Data,Limma)程序包,差異基因篩選條件為|log2 Fold change|>1及校正P<0.05 被認(rèn)為差異有統(tǒng)計學(xué)意義。

1.2.2加權(quán)共表達網(wǎng)絡(luò)構(gòu)建和模塊識別 WGCNA分析基因網(wǎng)絡(luò)服從無尺度分布的假設(shè),根據(jù)基因的表達相似性將其劃分為不同的模塊,并與表型進行關(guān)聯(lián)分析,最后識別感興趣的基因集。基于R語言“WGCNA”軟件包構(gòu)建了基于差異基因表達情況的加權(quán)共表達網(wǎng)絡(luò)[5]。分別定義了基因顯著性(gene significance,GS)及基因與模塊的相關(guān)性( module membership,MM)。GS和 MM的絕對值越大,代表基因與其所屬模塊以及臨床表型的相關(guān)性越大。

1.2.3篩選特征基因 基于R語言軟件中“glmnet”包對差異表達基因進行LASSO回歸分析,篩選RA的特征基因。lambdas增加,自由度和殘差減少,因此選取最小lambda值。然后通過glmnet函數(shù)進行交叉檢驗,選擇均方誤差最小時的 lambda值并輸出圖形。將lambda帶入模型,篩選最終輸出的變量作為RA的特征基因[6]。通過R語言軟件在數(shù)據(jù)集中繪制受試者工作特征曲線(receiveroperator characteristic curve,ROC)曲線,為了進一步驗證這5個特征基因,計算ROC曲線下的面積(area under curve,AUC)評估特征基因用于RA診斷的價值,且使用標(biāo)準(zhǔn)正態(tài)分布計算了AUC的95%置信區(qū)間(confidence interval,CI),AUC>0.80 且P<0.05的基因作為后續(xù)分析的重點。

1.2.4特征基因與免疫細胞的相關(guān)性分析 在TISIDB(http://cis.hku.hk/TISIDB/downl)下載免疫細胞的數(shù)據(jù)集,包括28 種免疫細胞,分為適應(yīng)性免疫細胞類和先天性免疫細胞類。通過R軟件相關(guān)程序輸入免疫細胞浸潤文件、4個GEO數(shù)據(jù)集合并文件,繪制免疫特征基因與TISIDB數(shù)據(jù)庫中免疫細胞的可視化圖。

1.2.5UHPLC-Q-Exactive Orbitrap MS鑒定分析土茯苓黃酮類化合物 按照文獻中的方法進行操作[7],超高效液相色譜條件:UPLC BEH C18 色譜柱(100 mm×2.1 mm,1.7 μm)。進樣體積為5 μl。流速單位為μl/min,A和B相中均加入0.1%甲酸。洗脫程序如下:0 min,5%B相;3.5 min,15%B相;6 min,30%B相;12 min,70%B相;18 min,100%B相;26 min,5%B相;30 min,5%B相,流速為0.3 ml/min,進樣量為1 μl,柱溫為40 ℃。質(zhì)譜條件:采用電噴霧離子源(ESI),正離子與負離子兩 種掃描模式,掃描模式為全掃描/數(shù)據(jù)依賴的二級掃描(full/ddMS2)。質(zhì)量掃描范圍m/z 為100～1 200,毛細管溫度350 ℃,鞘氣流速30 arb,輔助氣流速10 arb,MS2采用低、中、高3種碰撞能。

1.2.6分子對接分析 采用分子對接軟件(SYBYL 2.1.1)中的 Surflex-Docking 模塊分析土茯苓黃酮類化合物與5個特征基因的結(jié)合情況,結(jié)果以總分值(total_ score,T_Score)>5為閾值進行分析[8]。從蛋白質(zhì)晶體結(jié)構(gòu)數(shù)據(jù)庫RCSB (http://www.rcsb.org/pdb)中獲取載脂蛋白D(apolipoprotein D,APOD,PDB ID:2HZQ),含鋅指和 BTB 域 16 (zinc finger and BTB domain containing 16,ZBTB16,PDB ID:1BUO),趨化因子C-C亞族受體5(C-C chemokine receptor type 5,CCR5,PDB ID:7F1R),基質(zhì)金屬蛋白酶-1(matrix metalloproteinase 1,MMP1,PDB ID:3SHI)和冠蛋白1A(coronin-1A,CORO1A,PDB ID:2AQF)晶體結(jié)構(gòu)。

2 結(jié)果

2.1 GEO數(shù)據(jù)集處理及差異基因分析4個數(shù)據(jù)集合并作為訓(xùn)練數(shù)據(jù)集(表1),正常對照組32個樣本,RA組44個樣本。以PCA聚類圖的形式展示消除批次之間的差異前和后,PCA結(jié)果顯示,圖1A為矯正前的GEO數(shù)據(jù)集的分布情況,圖1B為矯正后的GEO數(shù)據(jù)集的分布情況,在數(shù)據(jù)矯正前4個數(shù)據(jù)集有明顯差異,而批次矯正后基本上消除了4個數(shù)據(jù)集的批處理效應(yīng)。

圖1 批次矯正前后PCA結(jié)果

表1 RA基因芯片基本信息

2.2 關(guān)鍵模塊與RA的相關(guān)性借助 R 軟件的 “WGCNA”軟件包對 4個GEO數(shù)據(jù)集的表達譜數(shù)據(jù)構(gòu)建加權(quán)共表達網(wǎng)絡(luò),通過 WGCNA 分析將差異表達的基因構(gòu)建4個不同的模塊。圖2A 顯示了藍色模塊、藍綠色模塊、棕色模塊、黃色模塊與RA的相關(guān)性,其中藍色模塊和藍綠色模塊與RA的相關(guān)性較高(r=0.75、0.73),且藍色模塊和藍綠色模塊中基因與RA表型呈正相關(guān)(r=0.87、0.91,P<0.05)(圖2B、C),由于藍色模塊中基因與RA表型相關(guān)性低于藍綠色模塊,說明這些基因在藍綠色模塊中具有重要角色,因此本研究選擇藍綠色模塊中的基因進行深入研究。

圖2 基因模塊與類風(fēng)濕關(guān)節(jié)炎的相關(guān)性分析圖

2.3 類風(fēng)濕關(guān)節(jié)炎的特征基因?qū)⑺{綠色模塊中連接性最高的69個基因作為候選中心基因。4個數(shù)據(jù)集合并處理后鑒定差異表達基因共63個,圖3A顯示有9個交集基因。圖3B為LASSO回歸的10折交叉驗證圖,顯示基于LASSO回歸算法對9個基因進行篩選,圖中的每一條曲線代表了每個自變量系數(shù)的變化軌跡,縱坐標(biāo)是系數(shù)的值,下橫坐標(biāo)是L1范數(shù),上橫坐標(biāo)是此時模型中非零系數(shù)的個數(shù)。參數(shù)系數(shù)絕對值隨著L1范數(shù)值變小而變大。當(dāng)L1范數(shù)達到一定值以后,一部分不重要的變量(特征基因)將被壓縮為0,代表該變量已被剔除。圖3C為最佳調(diào)諧參數(shù)選擇圖,得到交叉驗證曲線和最小化平均交叉驗證誤差的lambda的值。共篩選出5個對RA疾病具有診斷價值的特征基因,分別為APOD、ZBTB16、CCR5、MMP1和CORO1A。

圖3 類風(fēng)濕關(guān)節(jié)炎特征生物標(biāo)志物的鑒定和LASSO回歸圖

2.4 特征基因的表達情況基因CCR5、APOD、ZBTB16、MMP1和CORO1A在數(shù)據(jù)集中均存在差異表達(圖4)。與正常對照組相比,RA組CCR5、MMP1和CORO1A表達增高 (P<0.001 )。APOD和ZBTB16表達降低(P<0.001 )。

圖4 正常對照組和RA組特征基因表達的比較

2.5 特征基因的 ROC曲線分析結(jié)果從標(biāo)準(zhǔn)化表達矩陣中獲取CCR5、APOD、ZBTB16、MMP1和CORO1A基因的表達情況,結(jié)果顯示5個特征基因的AUC均大于0.85(圖5)。

圖5 特征基因的ROC曲線

2.6 特征基因與免疫細胞浸潤之間的關(guān)系通過 R 語言繪制5個特征基因與免疫細胞的相關(guān)性熱圖,縱坐標(biāo)代表免疫細胞,橫坐標(biāo)代表5個特征基因,模塊中顏色代表5個特征基因與免疫細胞間的相關(guān)性,模塊中顏色深淺代表該免疫特征基因與免疫細胞間的相關(guān)性的強弱,藍色代表負相關(guān),紅色代表正相關(guān)(圖6)。發(fā)現(xiàn)CCR5、MMP1和CORO1A與多種免疫細胞呈正相關(guān),其中CCR5、MMP1和CORO1A均與伽馬三角洲T細胞(γδT細胞)呈正相關(guān)(均P<0.05)。CCR5和CORO1A與巨噬細胞、調(diào)節(jié)性 T 細胞、單核細胞、骨髓源性抑制細胞呈正相關(guān)(均P<0.05)。而APOD基因與多種免疫細胞呈負相關(guān)(P<0.05)。

圖6 特征基因與免疫細胞的相關(guān)性熱圖

2.7 土茯苓黃酮類化合物的基本信息檢測數(shù)據(jù)經(jīng)峰匹配與校準(zhǔn)后共得到133個化合物離子用于后續(xù)分析。通過與標(biāo)準(zhǔn)品、保留時間、精確分子量和二級質(zhì)譜比對共結(jié)構(gòu)鑒定出20個化合物(表2)。

表2 鑒定出的土茯苓黃酮類化合物

2.8 土茯苓黃酮類化合物與特征基因的結(jié)合情況通過分子對接發(fā)現(xiàn)14個化合物與MMP1有較好的結(jié)合(T_Score>5);16個化合物與CCR5有較好的結(jié)合(T_Score>5);17個化合物與APOD有較好的結(jié)合(T_Score>5);8個化合物與CORO1A有較好的結(jié)合(T_Score>5);3個化合物與ZBTB16有較好的結(jié)合(T_Score>5);同時發(fā)現(xiàn),Mulberrin和Neobavaisoflavone化合物能同時與5個特征基因結(jié)合(表3),有進一步研究的價值。

表3 土茯苓黃酮類化合物與特征基因結(jié)合總分

3 討論

RA是一種常見的自身免疫性疾病,能反復(fù)和持續(xù)激活先天性和獲得性免疫系統(tǒng),導(dǎo)致免疫耐受失敗、大量炎性細胞因子產(chǎn)生、致殘性軟骨和骨骼損傷,也會出現(xiàn)肺部和心臟的全身性炎癥反應(yīng)。本研究通過生物信息學(xué)方法,對GEO數(shù)據(jù)庫中RA相關(guān)數(shù)據(jù)深度挖掘,發(fā)現(xiàn)5個RA特征基因(CCR5,APOD,ZBTB16,MMP1,CORO1A),且5個RA特征基因與多種免疫細胞相關(guān)。RA的發(fā)病過程中趨化因子通過調(diào)節(jié)免疫和炎癥反應(yīng)發(fā)揮介質(zhì)作用。CCR5是一種G蛋白偶聯(lián)受體。已有研究[9]發(fā)現(xiàn)CCR5參與了包括RA、骨質(zhì)疏松癥在內(nèi)的多種疾病的發(fā)病。滑膜T細胞中的CCR5水平升高,并與RA的發(fā)生發(fā)展有關(guān),CCR5在外周血單核細胞和 T 細胞中的表達水平低,但在體內(nèi)外經(jīng)炎癥刺激后表達明顯增強;APOD是多功能蛋白質(zhì),涉及炎癥、抗氧化等。在多種組學(xué)中APOD 表達影響糖尿病的發(fā)生和發(fā)展[10]。但其在RA的表達情況尚未明確,值得進一步深入研究;在基質(zhì)金屬蛋白酶(matrix metalloproteinases,MMPs)家族中,MMP1 屬于膠原酶和間質(zhì)溶解素類,并且有研究[11]證明MMP1的水平與 RA 疾病活動度相關(guān),并且可以預(yù)測RA 的功能和影像學(xué)結(jié)果,是預(yù)測 RA 早期骨破壞的指標(biāo)。CORO1A基因參與自身免疫病的發(fā)生發(fā)展,其是特異性調(diào)節(jié)因子,參與維持外周初始 T 細胞穩(wěn)態(tài)[12]。以上提示這些特征基因在RA的發(fā)生和發(fā)展中發(fā)揮重要作用。

RA發(fā)病機制復(fù)雜,多種免疫細胞參與RA的發(fā)病及進展過程,通過探索RA特征基因與免疫細胞之間的關(guān)系,可以更好地了解RA的免疫學(xué)機制。本研究篩選特征基因與免疫細胞的相關(guān)性發(fā)現(xiàn)CCR5、MMP1和CORO1A與多種免疫細胞呈正相關(guān),其中CCR5、MMP1和CORO1A均與γδT細胞呈正相關(guān)。CCR5和CORO1A與巨噬細胞、調(diào)節(jié)性 T 細胞、單核細胞呈正相關(guān)。單核細胞和巨噬細胞在 RA的發(fā)病機制中起著核心作用,其能通過產(chǎn)生不同種類的趨化因子、細胞因子和生長因子來調(diào)節(jié)免疫和炎性反應(yīng),導(dǎo)致骨和軟骨的破壞,誘導(dǎo)RA的發(fā)生[13];調(diào)節(jié)性T細胞(Treg)在多種疾病中發(fā)揮作用。是輔助性T細胞的特殊亞群之一,Treg能有效防止自身免疫性疾病的發(fā)展。RA患者的Treg細胞以功能受損和表型改變?yōu)樘卣鱗14];γδT細胞是免疫系統(tǒng)中的一類特殊細胞,已有研究[15]顯示γδT細胞在RA發(fā)病中具有啟動和保護雙重作用。RA作為自身免疫性疾病,其自身反應(yīng)性 CD4+T細胞、巨噬細胞、炎癥細胞因子、趨化因子等水平異常升高。而本研究發(fā)現(xiàn)特征基因與免疫細胞的相關(guān)性,對未來相關(guān)的科學(xué)研究提供理論依據(jù)。

綜上所述,基于差異表達數(shù)據(jù)利用LASSO回歸篩選出5個RA特征基因,在 ROC曲線分析中,CCR5,APOD,ZBTB16,MMP1,CORO1A的AUC值均大于0.85,這5個RA特征基因作為RA診斷標(biāo)志物具有潛在價值。此外,本研究選發(fā)現(xiàn)土茯苓黃酮類化合物與RA診斷標(biāo)志物均有結(jié)合作用,其治療機制可能與免疫相關(guān),由于本研究尚未進行體內(nèi)實驗研究,具體機制有待進一步體內(nèi)實驗驗證。