類風濕關節炎免疫特征基因及其與土茯苓黃酮類化合物的關系研究

楊 欣,黃 聰,姚血明

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

1 材料與方法

1.1 材料

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

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

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

1.2 方法

1.2.1基因芯片處理 借助美國國立生物技術信息中心的Gene Expression Omnibus(GEO)數據庫檢索 “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 被認為差異有統計學意義。

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

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

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

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 ℃。質譜條件:采用電噴霧離子源(ESI),正離子與負離子兩 種掃描模式,掃描模式為全掃描/數據依賴的二級掃描(full/ddMS2)。質量掃描范圍m/z 為100～1 200,毛細管溫度350 ℃,鞘氣流速30 arb,輔助氣流速10 arb,MS2采用低、中、高3種碰撞能。

1.2.6分子對接分析 采用分子對接軟件(SYBYL 2.1.1)中的 Surflex-Docking 模塊分析土茯苓黃酮類化合物與5個特征基因的結合情況,結果以總分值(total_ score,T_Score)>5為閾值進行分析[8]。從蛋白質晶體結構數據庫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),基質金屬蛋白酶-1(matrix metalloproteinase 1,MMP1,PDB ID:3SHI)和冠蛋白1A(coronin-1A,CORO1A,PDB ID:2AQF)晶體結構。

2 結果

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

圖1 批次矯正前后PCA結果

表1 RA基因芯片基本信息

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

圖2 基因模塊與類風濕關節炎的相關性分析圖

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

圖3 類風濕關節炎特征生物標志物的鑒定和LASSO回歸圖

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

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

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

圖5 特征基因的ROC曲線

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

圖6 特征基因與免疫細胞的相關性熱圖

2.7 土茯苓黃酮類化合物的基本信息檢測數據經峰匹配與校準后共得到133個化合物離子用于后續分析。通過與標準品、保留時間、精確分子量和二級質譜比對共結構鑒定出20個化合物(表2)。

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

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

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

3 討論

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

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

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