联合机器学习和网络药理学探究小柴胡汤治疗类风湿关节炎的作用机制*

陈丝梦,孟祥文,王金,程芷洛,宣自华,贾晓益

1.安徽中医药大学,安徽 合肥 230012; 2.安徽省中药复方重点实验室,安徽 合肥 230012;3.安徽省中药研究与开发重点实验室,安徽 合肥 230012

类风湿关节炎(rheumatoid arthritis,RA)是一种以滑膜增生和免疫细胞浸润为特征的自身免疫性疾病,临床特征为进行性关节肿胀、疼痛和关节畸形,并逐渐出现关节功能障碍,最终导致残疾[1]。全球疾病负担研究报告发现,RA作为一种顽固性疾病,在所有关节疾病中致残率最高,其全球发病率、患病率和发病年龄均有所增加,患病率已达0.27%[2]。目前,临床治疗RA的常用药物存在一定的不良反应,部分药物作用机制单一,不能从根本上阻止RA的发展进程[3-4]。因此,寻找疗效显著且不良反应小的药物尤为重要。

根据中医理论,RA可归于“痹证”范畴[5]。本病活动期基本病机为脏腑内蕴湿热毒,外以风、寒、湿等邪从化,最终导致湿热毒邪痹阻经络、筋脉骨节,湿热为主要病理因素[6]。中药具有多组分、多途径、多靶点作用的特点,不良反应少,在防治RA方面可发挥重要作用。小柴胡汤记载于《伤寒论》,由柴胡、黄芩、人参、半夏、甘草、生姜和大枣七味药材组成,广泛用于治疗多种疾病,研究表明,其在增强免疫、抗血管生成和抗感染与解热等方面疗效显著[7-9]。本研究应用机器学习和网络药理学筛选和预测小柴胡汤活性成分在治疗RA中的潜在靶点和信号通路,为探讨小柴胡汤治疗RA的作用机制提供分子基础,并为小柴胡汤抗RA药物的进一步发掘提供了理论依据。

1 资料与方法

1.1 药物成分获取与成分靶点搜集基于中药系统药理学数据库和分析平台(Traditional Chinese Medicine Systems Pharmaclogy Database and Analysis Platform,TCMSP)获取小柴胡汤的化学成分,认为口服生物利用度(oral bioavailability,OB)≥30%且类药性(drug likeness,DL)≥0.18的成分是活性成分。利用Swiss Target Prediction和HERB数据库预测活性成分靶点并合并两个数据平台靶点。

1.2 GEO芯片下载与差异基因分析GEO(Gene Expression Omnibus)平台储存着大量测序芯片数据。我们利用R软件包GEOquery获取RA相关数据集GSE93272的表达谱和临床信息。运用R软件包limma对数据进行差异分析,|logFC|>0.3、P<0.05的基因被认为是疾病的差异靶点基因。通过ggplot2绘制差异分析后的火山图。

1.3 支持向量机分析支持向量机(support vector machine,SVM)是一类对数据进行二元分类的广义线性分类器,属于机器学习中的监督学习。我们利用SVM进一步预测靶点基因对于疾病和正常状态的影响。Rank值大于1 000的差异靶点基因被定义为疾病靶点基因。疾病靶点基因与药物靶点的交集靶点为小柴胡汤治疗RA的潜在作用靶点,结果采用韦恩图展示。

1.4 富集分析利用clusterProfiler包进行交集靶点的基因本体(Gene Ontology,GO)功能富集分析和WikiPathway通路富集分析,富集方法为超几何分布。enrichplot包被用来绘制富集分析柱状图。其中GO功能分析包括生物过程(biological process,BP)、细胞组分(cell component,CC)和分子功能(molecular function,MF)。

1.5 “药物-活性成分-潜在靶点-疾病”网络构建将筛选得到的活性成分和潜在作用靶点导入Cytoscape 3.8.2软件,构建“药物-活性成分-潜在靶点-疾病”网络图。图中节点(node)分别为药物、活性成分、潜在作用靶点,边(edge)表示各节点之间的相互关系。利用插件CytoNCA对网络图进行分析,获取小柴胡汤治疗RA的核心成分。

1.6 LASSO回归分析LASSO回归在拟合广义线性模型的同时进行变量筛选,同时利用λ参数控制模型的惩罚力度。利用LASSO回归分析筛选小柴胡汤治疗RA的核心靶点基因。

1.7 核心靶点表达量与诊断分析将核心靶点基因在样本中的表达量提取后,通过箱线图可视化核心靶点在正常样本和RA患者样本中的表达差异。核心靶点基因在疾病中的诊断效能通过受试者工作特征曲线(receiver operating characteristic curve,ROC)呈现。

1.8 免疫浸润分析ssGSEA算法可以分析28种免疫细胞在样本中的浸润情况,利用ssGSEA算法分析RA中免疫细胞的浸润情况和核心靶点基因与免疫细胞的关系。

1.9 分析流程全文整体分析流程见图1。

图1 研究流程图

2 结果

2.1 小柴胡汤活性成分与相关靶点基于TCMSP收集小柴胡汤193个活性成分,其中半夏13个、柴胡17个、大枣29个、甘草92个、黄芩36个、人参22个、生姜5个,其中共同成分11个。基于Swiss Target Prediction和HERB数据库共收集到2 719个小柴胡汤活性成分的作用靶点。

2.2 差异基因分析通过limma包差异分析共筛选出1 361个差异基因,其中上调基因1 155个,下调基因206个。差异基因分析火山图如图2A所示。

注:A:差异基因火山图;B:SVM准确度曲线;C:SVM错误率曲线;D:药物靶点基因与疾病靶点基因交集韦恩图。

2.3 支持向量机分析为进一步明确RA的致病靶点基因,SVM-RFE被用来进一步分析。基于SVM构建的模型,选择Rank>1 000的基因作为疾病的靶点基因,共得到357个基因。模型的准确度和错误率展示在图2B—2C。小柴胡汤活性成分靶点和SVM筛选后RA靶点基因的交集被认为是药物治疗RA的潜在作用靶点,共得到62个交集靶点基因,见图2D。

2.4 富集分析采用clusterProfiler包对62个交集靶点进行GO富集分析和信号通路富集分析,以P<0.05为条件进行筛选后共得到489个GO条目(BP:383个;CC:54个;MF:52个),选择每项前5个GO条目绘制气泡图,见图3A。分析发现,BP主要富集在核染色体分离、长链脂肪酸生物合成过程、肽交联、前列腺素代谢过程、其他生物细胞杀伤正调节等;CC主要富集在黏着斑、细胞-基质连接、转移酶复合物,转移含磷基团、整合素复合物、参与细胞黏附的蛋白质复合物等;MF主要富集在生长因子结合、成纤维细胞生长因子结合、纤连蛋白结合、γ-微管蛋白结合、趋化因子结合等。根据P<0.05筛选出70条Wiki信号通路,选取前13条通路绘制气泡图,见图3B。分析发现,交集靶点主要富集在磷脂酰肌醇3-激酶(phosphatidylinositol 3-kinases,PI3K)-蛋白激酶B(protein kinase B,AKT)信号通路、细胞周期、血小板源生长因子(platelet-derived growth factor,PDGF)信号通路等。

注:A:交集靶点基因GO功能富集分析;B:交集靶点基因通路分析。

2.5 “药物-活性成分-潜在靶点-疾病”网络构建将筛选得到的活性成分和潜在作用靶点导入Cytoscape 3.8.2软件构建“药物-活性成分-潜在靶点-疾病”网络,见图4。网络图中,紫色代表中药,黄色代表活性成分,绿色代表潜在靶点,红色代表疾病。为进一步探究小柴胡汤治疗RA的核心成分,利用Cytoscape中CytoNCA插件进行度值排序,结果发现,度值排序最靠前的活性成分为香豆雌酚(度值为38),这提示香豆雌酚可能是小柴胡汤治疗RA时发挥作用的核心成分。

图4 “药物-活性成分-潜在靶点-疾病”网络图

2.6 LASSO回归分析为进一步明确小柴胡汤治疗RA的核心作用靶点,我们运用LASSO回归构建模型筛选关键作用靶点。62个基因的权重系数被展示在图5A;62个基因的模型可解释偏差被展示在图5B;62个基因的LASSO系数分布见图5C,其中模型惩罚力度设置为最小,即λ为最小λ值0.154。模型建立后的ROC曲线见图5D,曲线下面积(area under curve,AUC)为0.884,表明模型可靠。通过LASSO回归筛选到2个核心靶点基因,分别为趋化素样因子(chemokine like factor,CKLF)和胱抑素 A(cystatin A,CSTA),提示这两个靶点可能是小柴胡汤治疗RA的核心作用靶点。

注:A:调整参数选择的交叉验证 (权重系数);B:调整参数选择的交叉验证 (可解释偏差);C:交集靶点基因的LASSO系数分布;D:SVM的ROC预测。

2.7 核心靶点表达和诊断分析提取CKLF和CSTA在GSE93272中的表达数据和分组信息用于分析组间差异,结果显示,与正常样本组比较,RA患者样本中CKLF和CSTA均呈现高表达趋势(P<0.000 1)。通过ROC分析展示CKLF和CSTA对于RA的诊断效能,结果显示,CKLF和CSTA在RA中具有良好的诊断效能,AUC分别为0.876和0.884。见图6。

注:A:CKLF在正常样本和RA样本中的表达差异;B:CSTA在正常样本和RA样本中的表达差异;C:CKLF的诊断ROC曲线;D:CSTA的诊断ROC曲线。

2.8 免疫浸润分析RA的病变过程与免疫细胞的变化密不可分,我们通过ssGSEA算法探究核心靶点基因与28种免疫细胞的关系。免疫细胞在每个样本中的评分热图见图7A。正常样本和RA患者样本中的免疫细胞差异被展示在图7B,结果显示,活化的CD8+T细胞(activated CD8+T cell)、中央记忆型CD8+T细胞(central memory CD8+T cell)、效应记忆型CD8+T细胞(effector memeory CD8+T cell)、活化的CD4+T细胞(activated CD4+T cell)、中央记忆型CD4+T细胞(central memory CD4+T cell)、效应记忆型CD4+T细胞(effector memeory CD4+T cell)、γδ T细胞(gamma delta T cell)、Th17细胞(type 17 T helper cell)、Th2细胞(type 2 T helper cell)、记忆B细胞(memory B cell)、CD56bright自然杀伤细胞(CD56bright natural killer cell)、自然杀伤性T细胞(natural killer T cell)、活化的树突状细胞(activated dendritic cell)、浆细胞样树突状细胞(plasmacytoid dendritic cell)、嗜酸性粒细胞(eosinophil)、肥大细胞(mast cell)和中性粒细胞(neutrophil)具有显著差异,提示RA的发生发展可能与这些免疫细胞失调相关。在此基础上,进一步将RA患者数据分别再分为CKLF高、低表达组和CSTA高、低表达组。结果发现,CKLF的高、低表达与活化的CD8+T细胞、中央记忆型CD8+T细胞、效应记忆型CD8+T细胞、自然杀伤性T细胞、浆细胞样树突状细胞、嗜酸性粒细胞、肥大细胞和中性粒细胞等免疫细胞浸润显著相关,见图7C;而CSTA的高、低表达与活化的CD8+T细胞、中央记忆型CD8+T细胞、活化的CD4+T细胞、中央记忆型CD4+T细胞、γδ T细胞、效应记忆型CD4+T细胞等免疫细胞浸润显著相关,见图7D。对CKLF和CSTA与免疫细胞的相关性分析发现,与CKLF最相关的免疫细胞是嗜酸性粒细胞(r=0.59)和中性粒细胞(r=0.52);与CSTA最相关的免疫细胞是效应记忆型CD4+T细胞(r=0.64)和rδ T细胞(r=0.62),见图8。以上结果表明,CKLF和CSTA可通过调控多种免疫细胞的浸润进而缓解RA的进程。

注:A:28种免疫细胞在样本中的浸润评分热图;B:28种免疫细胞在正常样本和RA患者样本中的浸润差异;C:28种免疫细胞在CKLF高、低表达组间的浸润差异;D:28种免疫细胞在CSTA高、低表达组间的浸润差异。

注:A:CKLF和28种免疫细胞的相关性;B:CSTA和28种免疫细胞的相关性。

3 讨论

本研究利用机器学习和网络药理学探讨多组分、多途径、多靶点的中药复方小柴胡汤治疗RA的可能分子机制。首先,基于“药物-活性成分-潜在靶点-疾病”网络发现,香豆雌酚具有较高的度值,表现出显著的抗RA作用。香豆雌酚是一种天然存在的强效植物雌激素,属于Coumestans家族,广泛存在于大豆和球芽甘蓝中[10]。越来越多的研究表明,香豆雌酚具有强大的抗炎、抗氧化特性。香豆雌酚作为Coumestans家族的成分之一,其可显著减少脂多糖诱导的RAW264.7细胞中一氧化氮(nitric oxide,NO)、前列腺素E2(prostaglandin E2,PGE2)和活性氧(reactive oxygen species,ROS)等炎症介质的产生,从而发挥抗炎作用[11]。同时,香豆雌酚可降低促炎细胞因子IL-1β刺激的原代大鼠软骨细胞中基质降解酶如基质金属蛋白酶-13(matrix metalloproteinase-13,MMP-13)、MMP-3和MMP-1的表达及一氧化氮合酶(nitric oxide synthase,NOS)、环氧合酶-2(cyclooxygenase-2,COX-2)、PGE2和炎性细胞因子等的水平,其机制可能是通过拮抗原代培养大鼠软骨细胞中的IL-1β来实现的[12]。Jin等[13]用甲状腺球蛋白(thyroglobulin,TG)和弗氏佐剂诱导实验性自身免疫性甲状腺炎(experimental autoimmune thyroiditis,EAT)小鼠模型并给予香豆雌酚干预,结果发现,香豆雌酚治疗后小鼠血清抗TG-IgG2a、IgG3、IgG1滴度、IgG2a/IgG1比值及脾细胞中辅助性T细胞1(helper T cell 1,Th1)百分率均显著降低,同时小鼠脾组织中INF-γ的水平显著降低,表明香豆雌酚可以通过抑制Th1反应来减少EAT发生过程中甲状腺自身抗体的产生并在一定程度上延缓了该疾病的进程。然而,香豆雌酚在RA中的作用尚未被研究。RA是一种经典的自身免疫性疾病,且RA的发生发展与过度分泌的炎性因子和软骨细胞功能异常密切相关,这突出了香豆雌酚作为RA治疗重要成分的潜力。

其次,筛选得到62个潜在治疗靶点,并根据LASSO回归进一步证明CKLF和CSTA是小柴胡汤治疗RA的核心靶点。CKLF是第一个被发现的CMTM家族成员,其至少包括CKLF1、CKLF2、CKLF3、CKLF4四种形式。目前,已有大量研究表明,CKLFs在炎症性疾病和自身免疫性疾病中发挥重要作用。早期有学者通过反转录-聚合酶链反应和免疫组织化学反应发现CKLF1在RA患者滑膜中高表达[14]。He等[15]通过R软件limma包进行差异分析发现,CKLF是PBMC中可能与RA相关的差异表达基因,且体内实验证实了CKLF在RA患者PBMC中显著上调。CKLF1可与其功能受体CCR4结合,进而抑制软骨细胞增殖及促进内皮细胞增殖、炎性介质分泌和Th细胞等的趋化作用,从而参与包括RA在内的炎症反应[16]。有研究称 CKLF1 和CCR4在RA患者的滑膜中表达上调,且CKLF1在RA患者中的表达水平与C反应蛋白及血沉相关[17]。CSTA是具有多种生物活性的相关蛋白大家族Cystatins中的一员,是一种天然存在的单链组织蛋白酶抑制物,可在不同的细胞类型、组织和肿瘤中表达,同时也是研究最少的半胱氨酸氨基转移酶。Komura等[18]观察到胰腺导管腺癌患者外周血中CD4+T细胞中CSTA的基因水平高于健康志愿者,且患者血清中CSTA浓度也高于健康志愿者,提示CSTA在肿瘤细胞和肿瘤浸润性免疫细胞中表达,改变了局部肿瘤组织的炎症状态。然而,CSTA在RA中的表达和功能尚不清楚,但我们的结果为中医药防治RA提供新的分子机制。在本研究中,我们进一步使用ssGSEA算法分析了28个免疫细胞在RA样本中的渗透情况。与正常样本比较,RA患者样本PBMC中活化的CD8+T细胞、中央记忆型CD8+T细胞、效应记忆型CD8+T细胞、活化的CD4+T细胞、中央记忆型CD4+T细胞、效应记忆型CD4+T细胞、γδ T细胞、Th17、Th2、记忆B细胞、CD56bright自然杀伤细胞、自然杀伤性T细胞、活化的树突状细胞、浆细胞样树突状细胞、嗜酸粒细胞、肥大细胞和中性粒细胞水平发生显著变化。有研究显示,CKLF1拮抗肽C19在肝脏缺血再灌注期抑制CKLF1表达,可减轻肝脏缺血再灌注损伤,减少中性粒细胞浸润,减轻炎症反应[19]。胰腺导管腺癌患者外周血CD4+T细胞中CSTA表达上调[18]。即CKLF和CSTA的表达与许多免疫细胞的水平有关,这突显了它们作为RA治疗靶点的潜力。

最后,通过GO富集分析发现,小柴胡汤的抗RA作用靶点主要涉及参与细胞黏附的蛋白质复合体、核染色体分离、成纤维细胞生长因子结合等过程,体现了中医药治疗疾病的多途径的特点。WikiPathway的结果表明,小柴胡汤治疗RA的靶点集中在PI3K-Akt信号通路、NF-κB信号通路、PDGF通路等。相关研究表明,PI3K/AKT/mTOR信号通路参与RA患者滑膜液中信号素5A的改变,进而影响FLS活化的铁死亡进程[20]。NF-κB信号通路几乎存在于所有动物细胞中,在炎症、免疫反应、抗凋亡和促进细胞增殖等方面发挥重要作用。该信号通路的异常激活可导致M1表型巨噬细胞极化[21],炎性细胞因子和基质金属蛋白酶的水平上调及FLS的异常活化等[22-23],从而加重RA的发展进程。PDGF通路的失调可导致滑膜炎症、血管新生和RA的发生[24]。以上结果均提示,小柴胡汤可能作用于多种生物学过程和信号转导途径进而对RA有治疗作用。

综上,小柴胡汤通过多成分、多靶点、多途径来治疗RA,其中核心成分、核心靶点和核心途径分别为香豆雌酚、CKLF和CSTA及PI3K-Akt信号通路。本研究为今后中药方剂小柴胡汤治疗RA的机制研究提供了新方向,但已鉴定基因的生物学活性以及这些基因与RA的关系仍需进一步实验验证。