类风湿关节炎免疫特征基因及其与土茯苓黄酮类化合物的关系研究

杨 欣,黄 聪,姚血明

类风湿关节炎(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分析基因网络服从无尺度分布的假设,根据基因的表达相似性将其划分为不同的模块,并与表型进行关联分析,最后识别感兴趣的基因集。基于R语言“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、骨质疏松症在内的多种疾病的发病。滑膜T细胞中的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诊断标志物均有结合作用,其治疗机制可能与免疫相关,由于本研究尚未进行体内实验研究,具体机制有待进一步体内实验验证。