本研究采用整合的计算毒理学方法，系统性地探究了苯并[a]芘暴露与类风湿关节炎之间的关联。分析鉴定了五个核心靶基因，揭示了这些基因在关键免疫通路中的富集，并验证了BaP-蛋白质结合的稳定性，阐明了环境污染物诱导类风湿关节炎（RA）的潜在分子机制。
本文内容
摘要
摘要
多环芳烃（PAHs）是一类普遍存在的环境污染物，被认为是类风湿关节炎（RA）发病过程中重要的环境影响因素。苯并[a]芘（BaP）作为PAHs的关键组分，可能与RA的发病相关；然而，其潜在的毒理学机制尚未完全阐明。本研究采用网络毒理学、机器学习与分子对接相结合的综合计算方法，系统地填补了这一知识空白。首先，基于BaP的分子结构开展网络毒理学分析，通过整合并筛选多个数据库中的靶点信息，最终鉴定出15个BaP潜在的RA相关靶基因，并构建了其相互作用网络。GO和KEGG富集分析显示，这些基因显著富集于白细胞迁移、免疫细胞信号转导等生物学过程，并与NF-κB信号通路、T细胞受体信号通路等相关。随后利用STRING数据库和Cytoscape软件进行拓扑学分析，筛选出5个核心基因（LCK、ZAP70、ITK、GZMA和ITGAL），并通过机器学习进一步验证了这些基因的重要性。分子对接与分子动力学模拟结果表明，BaP与这些靶基因编码的蛋白产物具有较强的结合亲和力，能够形成构象稳定的复合物。综上所述，本研究通过整合计算方法揭示了BaP可能参与RA发生发展的潜在机制，为今后针对环境污染物相关RA的预防与治疗研究提供了理论依据。
引言
类风湿关节炎（RA）是一种常见的自身免疫性疾病，其特征为慢性滑膜炎和血管翳形成，进而导致软骨进行性破坏和骨侵蚀。这些病理变化可引起关节功能障碍、病理性骨折风险增加，并最终导致残疾，严重影响患者的生活质量1。尽管RA的病因尚未完全明确，但其发病机制通常归因于遗传因素（如HLA-DR4和HLA-DR1基因亚型）与环境因素（如吸烟、饮酒、EB病毒（Epstein-Barr virus）感染、空气污染暴露）的共同作用2。流行病学研究已证实，空气污染中的多环芳烃（PAHs）与RA发病风险升高之间存在显著关联3。
多环芳烃（PAHs）是一类常见的空气污染物，来源于煤炭、石油、天然气和烟草等物质的不完全燃烧，被认为是类风湿关节炎（RA）发病过程中重要的环境介导因子。PAHs可与免疫细胞中的芳香烃受体（AHR）复合物结合，导致核定位信号暴露，进而促使配体-AHR复合物转位至细胞核内。在RA发病机制中，与PAH结合的AHR作为关键的环境感应器，通过多个相互关联的通路驱动免疫系统失调。AHR激活后，会与ARNT形成异源二聚体，并调控CYP基因的表达，启动炎症级联反应4。同时，AHR的激活会破坏促炎性T细胞与调节性T细胞之间的平衡：一方面，其触发AHR/Jag1/Notch信号轴，增强Th17细胞的细胞因子释放5；另一方面，AHR直接结合GOT1启动子，上调GOT1表达，从而诱导FOXP3位点的高甲基化，抑制Treg细胞的分化6。由此导致的Th17/Treg失衡有利于促炎环境的形成。此外，AHR的激活还可通过上调CCR8表达以及IL-4、IL-13等细胞因子，进一步放大Th2型免疫反应，共同加剧滑膜炎症和组织损伤7。因此，AHR的激活构成了连接环境PAH暴露与Th17/Treg/Th2免疫失衡的关键枢纽，建立了基因型、环境因素与RA病理之间的机制性桥梁。
在大气中，多环芳烃（PAHs）以复杂的混合物形式存在，其中苯并[a]芘（BaP）是关键组分之一。在人巨噬细胞中，BaP可通过促进芳香烃受体（AHR）与CXCL8启动子的结合，诱导CXCL8（IL-8）的产生，进而上调中性粒细胞趋化因子的表达8。此外，BaP能够以剂量依赖的方式上调类风湿关节炎（RA）患者成纤维样滑膜细胞（FLS）中Slug的表达，从而加剧关节炎的进展9。在野生型小鼠中，BaP通过诱导CYP1A1酶活性，促进核因子κB受体活化因子配体（RANKL）介导的破骨细胞（OC）活化，最终导致骨量丢失10。然而，BaP毒性在RA发病机制中所起作用的具体机制仍不明确。我们假设，BaP通过直接作用于关键的免疫相关靶蛋白，干扰T细胞活化、Th17/Treg平衡以及炎症性细胞因子生成等多个信号通路，从而促进RA的发生发展，将环境中的BaP暴露与滑膜炎症及关节破坏联系起来。相较于传统的单一通路实验方法，网络毒理学更适用于本研究，因为BaP可能作用于多个免疫相关靶点及交叉通路。与通常一次仅研究一个通路或少数几个靶点的传统实验方法相比，网络毒理学能够系统性地揭示多靶点相互作用和整体效应，但其预测结果依赖于数据库信息，仍需通过实验加以验证。
目前关于BaP在类风湿关节炎（RA）中作用的研究大多局限于单一通路或线性机制描述，缺乏对多靶点、多层次网络调控特征的整合分析。因此，需要采用网络毒理学等整体性方法，以揭示BaP暴露与RA发病机制之间的复杂关联11,12,13。然而，针对环境污染物诱导疾病的网络毒理学研究仍十分匮乏。本研究的创新之处在于整合网络毒理学、机器学习与分子对接技术，系统性地探究BaP介导的RA发病机制，而非局限于单一通路或孤立靶点。研究通过拓扑学分析结合机器学习识别关键枢纽基因，并首次在分子水平上验证了BaP与各核心基因表达产物之间的结合模式及热力学稳定性。通过系统性地鉴定BaP可能促进RA发生与发展的潜在分子机制，本项计算研究旨在为理解RA的环境触发因素以及开发靶向治疗策略提供理论依据。与传统的单通路分析相比，这种整合方法能够系统评估多靶点相互作用，在研究复杂环境疾病机制方面具有更广泛的适用性。但需指出的是，本方法通过交集分析和拓扑分析优先筛选高置信度的枢纽基因，可能无意中排除了未同时满足筛选阈值但具有生物学意义的候选基因。未来的研究可探索互补策略，例如将机器学习应用于预测靶点的并集、整合其他组学数据或开展靶向实验验证，以进一步确认并拓展本研究的发现。
方案
伦理声明
本研究未直接涉及任何人类参与者或动物受试对象。
苯并[a]芘靶标获取
通过整合多个数据库的数据对苯并[a]芘（BaP）进行了表征。在PubChem数据库（https://pubchem.ncbi.nlm.nih.gov/）中使用“Benzo[a]pyrene”作为关键词检索，获取其化学结构和标准二维结构（SMILES字符串：C1=CC=C2C3=C4C(=CC2=C1)C=CC5=C4C(=CC=C5)C=C3）14。潜在的BaP靶标从ChEMBL（https://www.ebi.ac.uk/chembl/）、SEA（https://sea.bkslab.org/）和PharmMapper（http://lilab-ecust.cn/pharmmapper）数据库中获取15,16,17。所有预测的靶标均限定于智人（Homo sapiens）蛋白质组。预测的BaP靶标完整列表（n = 474）见补充表S1。完整的分析流程在图1中以示意图形式展示。
图1．本文数据集分析的流程图，展示了包括数据获取、预处理、差异表达分析、网络构建及验证步骤在内的整体工作流程。请点击此处查看该图的放大版本。
获取类风湿关节炎相关靶点
本研究通过在NCBI基因表达综合数据库（GEO，https://www.ncbi.nlm.nih.gov/gds/）中使用关键词“类风湿关节炎”和“智人”检索，获取了五个类风湿关节炎（RA）数据集18。根据数据集规模和实验设计，选择GSE77298（RA：16个样本；对照：7个样本）、GSE1919（RA：5个样本；对照：5个样本）和GSE55235（RA：10个样本；对照：10个样本）作为训练集，用于鉴定差异表达基因（DEGs）；而GSE12021（RA：24个样本；对照：13个样本）和GSE55457（RA：13个样本；对照：10个样本）则作为验证集。这些数据集的更多详细信息，如平台、样本及GSE编号，见表1。
使用 GEO2R 在线工具对数据进行标准化处理，生成 log2 转换后的表达矩阵用于后续分析。为消除不同实验批次带来的干扰，基于参数化经验贝叶斯框架，采用 SVA 软件包中的 ComBat 函数对数据集之间的系统性偏差进行校正。随后通过主成分分析（PCA）验证校正效果，结果显示批次间样本聚类显著改善，从而证实批次效应得到有效去除。合并并校正后的数据矩阵用于后续差异分析。
| GSE 系列 | 样本 | 平台 | 分组 |
| GSE77298 | 16 例类风湿关节炎和 7 例对照 | GPL570 | 训练队列 |
| GSE1919 | 5 例类风湿关节炎和 5 例对照 | GPL91 | 训练队列 |
| GSE55235 | 10 例类风湿关节炎和 10 例对照 | GPL96 | 训练队列 |
| GSE12021 | 24 例类风湿关节炎和 13 例对照 | GPL96 | 验证队列 |
| GSE55457 | 13 例类风湿关节炎和 10 例对照 | GPL9 | 验证队列 |
表1：本研究使用的五个GEO数据集的汇总信息。
该表格列出了每个数据集的GEO登录号（GSE系列）、样本组成（类风湿关节炎患者和健康对照的数量）、平台标识符（GPL），以及所属的训练队列或验证队列。
加权基因共表达网络分析（WGCNA）
采用加权基因共表达网络分析（WGCNA）评估与类风湿关节炎（RA）相关的差异表达基因（DEGs）的共表达网络特征19。基于校正批次效应后的表达矩阵，首先进行数据预处理：剔除标准差小于0.5的低变异基因，并使用评估优质样本和基因的函数对样本和基因质量进行评价。随后，采用层次聚类方法识别并剔除离群样本。为构建加权共表达网络，使用系统评估软阈值幂次的函数，对1至20范围内的软阈值幂次进行系统评估。最终选择幂次=12作为最优软阈值（无尺度拓扑拟合指数R²=0.90），以确保网络拓扑结构符合无尺度准则。基于该幂次值构建邻接矩阵，并计算拓扑重叠矩阵（TOM）。对基因进行层次聚类，并利用动态树切割算法识别初始基因模块。随后，通过对模块特征基因（module eigengenes）进行聚类，合并相似模块，从而构建稳健的基因模块网络。所有分析均使用专用于加权共表达网络分析的R软件包完成，以确保网络构建的可靠性与可重复性。进一步对差异表达基因/WGCNA枢纽基因与预测的BaP靶基因的交集进行分析，以鉴定与RA发病机制相关的BaP核心靶点，并利用Venn图软件进行可视化展示。
鉴定与类风湿关节炎发病机制相关的 BaP 靶点
使用 R 语言的 Venn 图分析包进行交集分析，以确定 BaP 靶点中与类风湿关节炎（RA）发病机制重叠的部分。将这些靶点导入 STRING 数据库，构建蛋白质-蛋白质相互作用（PPI）网络，物种设定为“Homo sapiens”，相互作用置信度评分设定为 > 0.7，以确保网络的高可靠性20。选择该阈值是因为其对应于 STRING 数据库中的“高置信度”水平，能够在保留生物学相关相互作用的同时，最大限度地减少通常与较低置信度评分相关联的假阳性结果。在毒理学网络研究中，采用 > 0.7 的截断值已被广泛采纳，用于优先筛选出稳健且可重复的蛋白质关联。从蛋白质-蛋白质相互作用数据库（STRING）下载生成的 TSV 文件，并导入网络可视化软件 Cytoscape 进行网络可视化分析。通过 CytoHubba 插件中的 Degree 算法对网络中的核心蛋白进行排序，并根据排序结果确定核心蛋白，用于后续分析。
KEGG 和 GO 富集分析
与 BaP 调控和 RA 发病机制相关的基因缩写通过 R 中的 "org.Hs.eg.db" 注释包转换为 Entrez ID。随后使用 clusterProfiler 工具进行 KEGG 通路富集分析，显著性阈值设定为 0.05。同时，GO 功能注释覆盖了三个主要的 GO 分类：生物过程（BP）、细胞组分（CC）和分子功能（MF），并使用 enrichGO 函数进行分析，P 值和 q 值的截断值均设定为 0.05。需要注意的是，本研究未进行多重检验校正，因为此次探索性分析的主要目的是最大限度地发现潜在相关的重要生物学通路和功能术语，从而为后续实验验证生成更广泛的可检验假设。最后，使用 enrichplot 包中的 barplot 和 dotplot 函数对富集分析结果进行图形化展示。
基于机器学习的核心基因验证
为评估与BaP和RA相关的核心基因的预测能力，并保持模型的透明性，我们实施了一套系统的机器学习工作流程。利用所选核心基因的表达谱，采用11种不同的机器学习算法构建了预测模型：Lasso回归（LR）、支持向量机（SVM）、随机森林（RF）、glmBoost、逐步广义线性模型（GLM）、岭回归、弹性网络（Enet）、梯度提升机（GBM）、线性判别分析（LDA）、极端梯度提升（XGBoost）以及朴素贝叶斯。通过五折交叉验证对超参数进行优化，并采用分层抽样将数据划分为训练集和内部验证集。在整个机器学习流程中使用固定的随机种子（set.seed(123)），以确保数据划分、交叉验证折次和模型训练的可重复性。各算法的关键超参数详见补充表S2。模型性能通过多个指标进行评估，包括曲线下面积（AUC）、准确率和F1分数。为克服单一模型方法固有的局限性，我们采用堆叠集成策略，整合表现最优的基础模型的预测结果。鉴于许多机器学习模型具有“黑箱”特性，我们应用SHapley加性解释（SHAP）算法量化每个基因对预测结果的贡献。通过SHAP值的大小和方向来解释基因在分类决策中的重要性，从而增强模型输出的可解释性。
BaP与核心靶点的分子对接
为探究BaP与核心基因产物之间的结合特性，进行了分子对接模拟。BaP（配体）的三维结构从PubChem数据库获取，格式为SDF。对应核心靶点的蛋白质结构从RCSB蛋白质数据库（https://www.rcsb.org/）以PDB格式下载，并根据其UniProt编号进行选择，优先选用包含共结晶配体或高分辨率坐标的结构。在对接前，使用PyMol进行蛋白预处理，移除水分子、共结晶配体以及离子等非蛋白组分，以避免干扰21。对于原始PDB结构中包含共结晶配体的蛋白，其活性位点中心根据结合配体的原子坐标确定；对于无共结晶配体的蛋白，则依据文献报道的关键残基坐标确定活性位点中心，这些残基通常对催化活性或抑制剂结合至关重要。对接网格以确定的活性位点坐标为中心，对每个靶点采用25 × 25 × 25 Å的立方体盒子。该标准的25 Å盒子尺寸可充分覆盖每个活性位点，并为配体采样提供足够空间，同时避免过高的计算成本。所有对接计算均使用AutoDock Vina（版本1.2.5）完成。选择Vina评分最优的构象作为代表性结合模式，并记录相应的结合能。使用PyMol（版本2.5.7）生成三维结合构象，使用Discovery Studio（版本2021）生成二维相互作用图，以可视化关键相互作用，包括氢键和疏水接触。
分子动力学模拟
使用 Gromacs 2025.3 进行分子动力学模拟，以对接获得的复合物作为初始结构。蛋白质原子采用 AMBER14SB 力场进行建模，水分子使用 TIP3P 模型表示。每个蛋白质-配体复合物被置于立方水盒子中溶剂化，蛋白质表面与盒子边界之间的最小距离为 1 nm。根据需要添加钠离子或氯离子以实现体系的电中性。首先结合最陡下降法和共轭梯度算法进行初始能量最小化，每种算法最多运行 10,000 步。长程静电相互作用通过粒子-网格埃瓦尔德（Particle-Mesh Ewald, PME）方法计算，范德华相互作用和短程静电相互作用均采用 1.0 nm 的截断距离。在能量最小化之后，体系在 NVT（恒定体积和温度）和 NPT（恒定压力和温度）条件下逐步进行平衡。随后在恒定温度和压力下进行 100 ns 的生产模拟，时间步长为 0.002 ps（2 fs），总共 50,000,000 步。每项模拟仅执行一次（无重复），主要目的是评估结合复合物在标准条件下的稳定性。温度通过 V-rescale 温控器维持，压力由 Parrinello–Rahman 压力耦合器控制。在整个模拟过程中，非键相互作用始终采用 1.0 nm 的截断距离。为了评估结构的稳定性和柔性，我们计算了原子位置的均方根偏差（RMSD）、每个残基的均方根涨落（RMSF）、作为结构紧密性度量的回转半径（Rg），以及溶剂可及表面积（SASA）。所有图表均使用 QtGrace 生成。
结果
BaP 靶点获取
BaP 的分子结构数据来自 PubChem 数据库（图 2A）。通过整合来自三个互补数据库——ChEMBL、PharmMapper 和 SEA 的信息，系统性预测了 BaP 的潜在生物学靶点，共鉴定出 474 个潜在靶点（图 2B）。
![多环芳烃结构，维恩图；苯并[a]芘化合物，基因重叠分析。](/files/ftp_upload/70636/70636fig2.jpg "figure-results-1")
图 2．与类风湿关节炎（rheumatoid arthritis，RA）相关的苯并[a]芘（BaP）靶蛋白的鉴定。（A）BaP 的化学结构，显示其分子构型。（B）通过 CHEMBL、PharmMapper 和 SEA 数据库进行靶点预测，展示用于鉴定潜在 BaP 相关靶点的工作流程。请点击此处查看该图的放大版本。
获取类风湿关节炎相关靶点
为了减少批次效应，使用 limma R 软件包中的 normalizeBetweenArrays 函数对三个数据集 GSE77298、GSE1919 和 GSE55235 进行合并和标准化处理。主成分分析（PCA）显示，标准化后的数据分布更加均匀，聚类模式也更为清晰（图 3A、B）。为定量验证标准化效果，我们计算了校正前后由批次因素解释的方差以及批次聚类的平均轮廓系数；两项指标在校正后均显著降低，表明成功消除了批次间的偏差。差异表达分析采用 |log₂FC| > 0.585（对应 1.5 倍变化）和校正后 p 值 < 0.05（Benjamini-Hochberg FDR 校正）作为阈值，共鉴定出 1379 个在类风湿关节炎患者中表达发生显著变化的基因。这些基因通过火山图和热图进行可视化展示（图 3C、D）。在加权基因共表达网络分析（WGCNA）中，共获得 13 个不同颜色标记的基因模块（图 3E）。模块-性状关联分析发现某些特定模块与类风湿关节炎之间存在显著相关性（p < 0.05）（图 3F）。通过将差异表达基因与 WGCNA 获得的基因取交集，最终鉴定出 500 个与类风湿关节炎相关的基因（图 3G）。
图3. (A主成分分析散点图显示，在批次校正前，GSE77298、GSE1919 和 GSE55235 数据集之间存在明显分离，表明存在批次效应。B) 经批次校正后的主成分分析散点图，显示三个数据集的整合显著降低了批次效应。C火山图显示基于logFC值和统计显著性筛选的差异表达基因。红色圆点表示上调基因，蓝色圆点表示下调基因，灰色圆点表示表达水平无统计学显著差异的基因。D热图显示不同样本中差异表达基因（DEGs）的表达模式；红色表示上调，蓝色表示下调。E) WGCNA 基因树状图，显示基于共表达关系的层次聚类；底部的颜色条代表不同的基因模块。F模块-性状关系热图，显示通过WGCNA鉴定的模块与样本性状（正常）之间的相关性 与类风湿性关节炎）；方框内的数值表示相关系数 P-值。 (G维恩图显示差异表达基因（紫色）与WGCNA分析所得基因（黄色）的分布情况；重叠的棕色区域表示两组数据共有的基因。 请点击此处以查看此图的放大版本。
与类风湿关节炎发病机制相关的 BaP 靶点鉴定
BaP 靶蛋白与类风湿关节炎相关基因的交集分析确定了 15 个可能在 BaP 诱导类风湿关节炎中起关键作用的靶点（图 4A、B，表 2）。如表 2所示，这些靶点包括 T 细胞信号分子（LCK、ZAP70、ITK）、免疫效应分子（GZMA、ITGAL），以及参与细胞骨架调控（CORO1A）、B 细胞信号传导（MS4A1）、细胞因子受体信号传导（ERBB4、IL1R1）、细胞周期（AURKA、KIF11）、DNA 修复（RAD51）、核受体信号传导（THRB、NR4A2）和蛋白水解（CTRC）的蛋白质。这 15 个基因随后被用于 GO 和 KEGG 富集分析（图 4C、D）。GO 分析显示，这些与 BaP 相关的靶点在白细胞迁移、免疫应答激活的细胞表面受体信号通路（生物过程）、免疫突触、膜筏、膜微区（细胞组分）、组蛋白激酶活性和蛋白酪氨酸激酶活性（分子功能）等生物学过程中显著富集。KEGG 通路分析表明，这些关键靶点主要与免疫信号转导、免疫效应过程及免疫缺陷相关通路有关，包括 NF-κB 信号通路、Th17 细胞分化/T 细胞受体信号通路以及自然杀伤细胞介导的细胞毒性。
图4．与类风湿关节炎（RA）相关的BaP靶基因的鉴定。（A）维恩图比较与BaP暴露（紫色）和RA（黄色）相关的基因，显示有15个重叠基因。（B）蛋白-蛋白相互作用网络，展示重叠基因之间的关系；节点代表基因，边代表预测的相互作用。（C）重叠基因在三个基因本体（Gene Ontology）类别中的功能注释：生物过程（BP）、细胞组分（CC）和分子功能（MF），表示基于显著性阈值富集的功能术语，展示最显著富集的通路。x轴表示基因数量；颜色梯度对应校正后的P值，红色越深表示显著性越高。（D）基于KEGG通路富集分析得到的与重叠基因相关的生物学通路。x轴表示基因比例，点的大小表示基因数量。颜色梯度对应校正后的P值，红色越深表示显著性越高。请点击此处查看该图的高清版本。
| 编号 | 蛋白质名称 | 基因名称 |
| 1 | Tyrosine-protein kinase Lck | LCK |
| 2 | Granzyme A | GZMA |
| 3 | Tyrosine-protein kinase ZAP-70 | ZAP70 |
| 4 | Tyrosine-protein kinase ITK/TSK | ITK |
| 5 | cDNA FLJ57691，与Integrin alpha-L高度相似 | ITGAL |
| 6 | Coronin | CORO1A |
| 7 | 跨膜4区A亚家族成员14 | MS4A1 |
| 8 | 表皮生长因子受体4 | ERBB4 |
| 9 | 白细胞介素1受体I型 | IL1R1 |
| 10 | Aurora kinase A | AURKA |
| 11 | DNA修复蛋白RAD51 | RAD51 |
| 12 | 驱动蛋白家族成员11 | KIF11 |
| 13 | 甲状腺激素受体β | THRB |
| 14 | 核受体亚家族4组A成员2 | NR4A2 |
| 15 | 糜蛋白酶C | CTRC |
表2：从蛋白质-蛋白质相互作用（PPI）网络中筛选出的15个核心靶基因列表。
针对每个基因，提供了相应的蛋白质名称。这些基因代表了苯并[a]芘（BaP）相关靶点与类风湿关节炎（RA）相关基因的交集。
基于机器学习的核心基因验证
根据蛋白质相互作用网络度值排序，针对前五个基因（LCK、ZAP70、ITK、GZMA 和 ITGAL）共构建了104个预测模型，以进行全面的机器学习分析。
这些模型通过对11种不同的机器学习算法在多个训练-测试划分、参数设置和特征预处理组合下进行评估而生成，共产生了104种独特的模型配置。为确保模型的稳健性，采用独立验证数据集和五折交叉验证程序对模型性能进行评估，作为预测准确性的内部对照。通过五折交叉验证评估了模型在不同运行间的性能变异性，曲线下面积（AUC）值以五折结果的均值±标准差（SD）表示。线性判别分析（LDA）模型在训练和验证阶段均表现出优异性能（图5A），在平均AUC和稳定性方面优于所有其他算法。受试者工作特征（ROC）曲线分析证实了这些核心基因具有良好的诊断价值（AUC > 0.8，图5B）。这些基因在类风湿关节炎（RA）中的差异表达模式通过火山图进行可视化展示（图5C）。采用SHAP可解释性分析对核心基因的贡献程度进行量化。根据所有样本的平均绝对SHAP值（图5D、E），LCK对模型预测的总体影响最强（0.158），其次是GZMA（0.095）、ITGAL（0.073）、ZAP70（0.034）和ITK（0.032）。对一个代表性样本的分析（图5G）进一步说明了各基因对单个预测结果的贡献：GZMA的正向贡献最大（+0.142），其次是ZAP70（+0.107）、ITGAL（+0.0684）、ITK（+0.0642）和LCK（+0.0537）。模型的基线值E[f(x)]为0.587，这些基因的联合贡献使预测值上升至f(x) = 0.958，这一变化在图中直观反映了所分析基因对预测性能的累积效应（图5D–G）。
图5. 与BaP诱导的RA相关核心基因的机器学习验证。(A) 模型性能（曲线下面积（AUC）热图，显示各模型间的性能比较），(B) ROC曲线（诊断效能），(C) 差异表达基因（DEG）火山图，(D) SHAP特征重要性，(E) SHAP汇总图，(F) SHAP依赖性图。(G) SHAP力图。请点击此处查看此图的高清版本。
苯并[a]芘与核心靶点的分子对接
为了确定苯并[a]芘（BaP）是否能够与五个核心基因编码的蛋白质结合，进行了分子对接分析。对接结果预测所有五种BaP-蛋白质复合物均具有有利的结合能，其值 consistently 低于 −7 kcal/mol：LCK = −7.9 kcal/mol，ZAP70 = −8.4 kcal/mol，ITK = −10.6 kcal/mol，GZMA = −7.5 kcal/mol，ITGAL = −10.2 kcal/mol（表3）。对结合构象的分析显示，每种BaP-蛋白质复合物均形成了稳定的对接构象（图6）。
| 配体 | 靶标名称 | PDB ID | 结合能 (kcal/mol) |
| Bap | LCK | P06239 | -7.9 |
| Bap | ZAP70 | P43403 | -8.4 |
| Bap | ITK | Q08881 | -10.6 |
| Bap | GZMA | P12544 | -7.5 |
| Bap | ITGAL | P20701 | -10.2 |
表3：分子对接分析结果，显示苯并[a]芘（BaP）与五个核心类风湿关节炎相关靶基因的蛋白产物之间预测的结合亲和力。
图6. BaP与类风湿关节炎（RA）相关核心基因之间相互作用的分子对接验证。BaP与GZMA（A）、ITGAL（B）、ITK（C）、LCK（D）和ZAP70（E）的对接结果，展示了预测的结合构象及相互作用模式。请点击此处查看该图的高清版本。
分子动力学模拟
为评估 BaP 与核心靶标蛋白形成的结合复合物的构象稳定性和动态行为，我们进行了分子动力学模拟。在初始构象调整阶段（0–20 ns）之后，RMSD 曲线趋于平稳，并在随后的 80 ns 模拟过程中仅出现微小波动（图 7A）。这表明复合物在模拟时间范围内已达到相对稳定的状态。值得注意的是，BaP–ZAP70 和 BaP–ITK 复合物的波动幅度最小，提示其在模拟中具有更优的结构稳定性。RMSF 分析揭示了蛋白残基水平上的灵活性差异（图 7B）。LCK、ZAP70 和 GZMA 激酶结构域中的关键残基表现出较低的 RMSF 值，表明 BaP 的结合在模拟中并未显著破坏其催化功能区域的结构刚性。相比之下，ITGAL 和 ITK 出现明显的 RMSF 峰值，提示非催化区域（如 N 端或 C 端）的残基在模拟环境中具有更高的灵活性，可能与其在蛋白质-蛋白质相互作用中的功能有关。对反映蛋白整体结构紧密性的 Rg 分析（图 7C）显示，所有复合物的 Rg 值在整个模拟过程中保持稳定，未出现持续上升或下降的趋势。SASA（图 7D）用于量化配体-蛋白复合物表面的溶剂可及性。尽管不同靶标蛋白因其天然结构差异导致绝对 SASA 值存在固有不同，但所有复合物的 SASA 值在模拟过程中均未发生剧烈变化，表明结合界面区域的溶剂可及性保持相对稳定。自由能景观分析（图 7E、F）从能量角度进一步支持了复合物在模拟条件下的稳定行为。所有体系均呈现出明显且集中的低能盆地，表明在模拟时间范围内，BaP 与各靶标蛋白的结合预计可形成热力学稳定的复合物。
图 7. BaP 与靶标蛋白的分子动力学模拟。（A）五种靶标蛋白与 BaP 结合过程中的均方根偏差（RMSD）曲线，反映各蛋白结构在模拟期间的稳定性变化。（B）五种靶标蛋白与 BaP 结合过程中的均方根波动（RMSF）曲线，反映蛋白残基水平的波动程度。（C）五种靶标蛋白与 BaP 结合过程中的回转半径（Rg）曲线，显示各蛋白整体构象紧密性的动态变化。（D）五种靶标蛋白与 BaP 复合物的溶剂可及表面积（Area）随时间演变图，横轴为模拟时间（ns），纵轴为表面积（nm2），反映构象稳定性和溶剂暴露特性。（E）五种靶标蛋白与 BaP 结合的自由能景观图，显示蛋白–BaP 结合的能量分布特征。（F）五种靶标蛋白与 BaP 结合的三维自由能景观图（FEL），展示不同构象状态之间的能量差异，反映蛋白–配体结合的能量景观特征。请点击此处查看该图的放大版本。
综上所述，这些结果表明，苯并[a]芘（BaP）与五个核心免疫相关靶点（LCK、ZAP70、ITK、GZMA 和 ITGAL）以及类风湿关节炎（RA）相关的关键通路（包括T细胞受体信号通路和NF-κB信号通路）存在关联。分子对接结果显示，所有BaP-靶点复合物均具有良好的结合能；分子动力学模拟进一步证实，这些复合物在100 ns时间尺度内具有构象稳定性。这些发现支持了以下假说：苯并[a]芘等环境污染物可能通过多靶点调控机制参与类风湿关节炎的发病过程。
讨论
类风湿关节炎（RA）是一种复杂的自身免疫性疾病，由遗传易感性与环境因素相互作用所致。在众多环境风险因素中，多环芳烃（PAHs）作为最常见的空气污染物之一，被认为是连接环境暴露与RA发病的重要环节。早期研究初步揭示，PAHs可通过芳香烃受体（AHR）通路影响免疫细胞分化的平衡，并诱导氧化应激。苯并[a]芘（BaP）是PAHs中具有高度代表性的组分，因其毒性强且在环境中广泛存在而受到广泛关注22。然而，BaP在特定细胞和信号网络中的直接作用靶点及核心作用机制尚不明确，这限制了对RA风险的预警能力以及针对环境因素的早期干预设计。为弥补这一空白，本研究提出假设：BaP可能通过直接结合五个核心免疫相关蛋白（LCK、ZAP70、ITK、GZMA和ITGAL），进而干扰T细胞受体、NF-κB和Th17信号通路，促进RA的发生发展。本研究通过整合网络毒理学、机器学习与分子对接技术，尝试构建一个系统的计算框架，突破单一通路分析的局限，为理解BaP通过多靶点调控机制参与RA发病提供理论依据。
通过这一整合分析，从最初的15个基因中鉴定出5个与类风湿关节炎（RA）相关的核心BaP靶基因：LCK、ZAP70、ITK、GZMA和ITGAL。LCK、ZAP70和ITK在T细胞受体信号起始过程中发挥关键作用；GZMA介导细胞毒性颗粒诱导的细胞凋亡；ITGAL则促进淋巴细胞黏附和跨内皮迁移。这些基因的共同失调可能加剧滑膜炎症。机器学习进一步验证了该基因组合的诊断潜力。SHAP分析显示，LCK对模型预测的总体贡献最大，其次是GZMA和ITGAL（图5D、E）。该排序突显LCK为最具影响力的预测因子，而GZMA和ITGAL也具有重要作用。分子对接模拟预测BaP与这些核心基因编码的蛋白质产物之间具有良好的结合亲和力。随后的分子动力学模拟表明，BaP-蛋白质复合物在100 ns的模拟过程中保持了构象稳定性，从动力学角度支持了预测的结合模式。
KEGG 和 GO 富集分析表明，鉴定出的 BaP 相关靶点富集于免疫相关通路，包括 T 细胞受体信号通路、NF-κB 信号通路、Th17 细胞分化以及自然杀伤细胞介导的细胞毒性（图 4C、D）。这些富集模式提示，BaP 暴露可能与类风湿关节炎（RA）发病机制相关的免疫过程有关。为进一步验证这一假设，可考虑采用其他方法，例如体外结合实验以直接测定 BaP 与蛋白质的相互作用、在 BaP 处理的免疫细胞中进行细胞功能实验，或建立暴露于 BaP 的 RA 动物模型。具体而言，LCK、ZAP70 和 ITK 是 TCR 信号级联中的关键组分，而 GZMA 和 ITGAL 参与免疫效应功能。预测的 BaP 与这些蛋白的结合，结合其在 RA 相关通路中的富集，提示了一种可能的作用机制：BaP 的结合可能改变这些激酶的构象或活性，从而导致异常的 TCR 信号传导、持续的 NF-κB 激活以及 Th17/Treg 分化失衡。类似地，BaP 与 GZMA 和 ITGAL 的相互作用可能影响细胞毒性颗粒的释放和淋巴细胞黏附，进而促进滑膜炎症。综上所述，这些发现将这五个靶点定位为 BaP 可能影响免疫过程的候选介导分子。
在免疫效应阶段，GZMA 和 ITGAL 通过不同的机制参与炎症反应的放大与维持。GZMA 是细胞毒性颗粒的主要成分，在活化的细胞毒性 T 淋巴细胞和自然杀伤细胞中发挥重要作用23。这与关于自身免疫中致病性免疫细胞亚群的广泛研究相一致24。BaP 可能通过影响 GZMA 的表达或功能，增强其对关节滑膜细胞和软骨细胞的细胞毒性作用，从而直接导致关节组织损伤。ITGAL 则是淋巴细胞趋化、黏附和跨内皮迁移的关键介质25。ITGAL 与 BaP 的相互作用可能促进炎症细胞向关节滑膜的募集和浸润，加剧局部炎症反应。上述过程可能共同汇聚于 NF-κB 信号通路。例如，T 细胞受体（TCR）信号的异常激活可能通过蛋白激酶 C（PKC）和 CARD11 等接头分子直接激活 NF-κB 通路26。此外，由细胞毒性作用和炎症细胞浸润引起的组织损伤可能通过模式识别受体间接激活该通路。NF-κB 通路激活后，可促进多种促炎细胞因子（如 TNF-α、IL-6 和 IL-1β）的转录与表达，从而形成慢性炎症环境。
综上所述，本研究通过整合网络毒理学、机器学习和分子动力学模拟，鉴定出 LCK、ZAP70、ITK、GZMA 和 ITGAL 为 BaP 可能影响与类风湿关节炎（RA）发病机制相关的免疫过程的候选靶点。本文提出的计算预测结果为进一步提出可验证的假设奠定了基础，用以探讨 BaP 等环境污染物在分子层面影响自身免疫性疾病的机制。对这些预测相互作用及其功能后果的实验验证仍是下一步的关键工作。
本研究存在一定的局限性。研究主要基于对公共数据库的生物信息学分析，其预测结果需要实验验证。尽管分子对接和分子动力学模拟提供了预测的结合能和构象稳定性，但这些方法无法获得严格的定量结合自由能估计值。因此，所报告的结合亲和力应被解释为基于对接的评分和模拟稳定性指标，而非绝对的自由能数值。未来采用自由能计算方法（如 MM/GBSA）的研究可为 BaP-蛋白相互作用提供更准确的估计值27。其次，对接分析未包含与已知不参与 BaP 相互作用的蛋白质进行基准比较。缺乏此类阴性对照限制了对结合特异性的评估。未来研究若能引入针对诱饵蛋白集或经实验验证不发生相互作用的靶标进行对接分析，将有助于明确所观察到的 BaP-蛋白结合的特异性。为进一步探究所提出的机制，未来研究可采取以下具体步骤：对 BaP-靶标复合物进行 MM/GBSA 计算以获得定量的结合自由能；纳入预期不与 BaP 结合的阴性对照蛋白以评估对接特异性；开展体外激酶活性实验，检测 BaP 暴露后 LCK 和 ZAP70 的磷酸化状态；建立 BaP 暴露的类风湿关节炎动物模型，以验证所鉴定靶标的体内相关性。
需要注意的是，我们选择通过交集分析差异表达基因（DEGs）和WGCNA衍生基因的方法，优先考虑特异性而非敏感性。该策略要求一个基因必须同时满足在类风湿关节炎（RA）样本与对照样本之间差异表达，并且属于与RA相关的稳健共表达模块，从而提高了所鉴定基因的可信度。相比之下，基于并集的方法（合并两种方法所得的所有基因）虽然可提高敏感性，但会引入更多假阳性结果，可能削弱后续分析的生物学相关性。鉴于本研究的目标是识别高可信度的候选靶点以进行机制研究和实验验证，因此采用交集策略更为合适。尽管我们也考虑了其他替代策略，例如对并集数据应用机器学习，或在特征筛选后进行交集分析，但这些方法未在本研究中实施，因其需要额外的参数优化和交叉验证步骤，超出了当前研究的范围。未来的研究可探索这些互补性方法，以进一步验证并拓展我们的研究发现。
结论：
本研究表明，苯并[a]芘（BaP）可能通过直接作用于五个核心免疫相关靶点（LCK、ZAP70、ITK、GZMA 和 ITGAL），并干扰 T 细胞受体、NF-κB 及 Th17 信号通路，从而参与类风湿关节炎（RA）的发病机制。结合网络毒理学、机器学习与分子动力学模拟的综合计算框架，为优先筛选环境污染物的作用靶点以供实验验证提供了系统性方法。这些发现为开发针对污染物诱发的自身免疫性疾病的早期预警生物标志物和靶向预防策略提供了理论依据。
数据可用性：
本研究中分析的所有原始数据集均可在 NCBI 基因表达综合数据库（GEO）中公开获取，登录编号分别为 GSE77298(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77298)、GSE1919(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1919)、GSE55235(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55235)、GSE12021(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12021) 和 GSE55457(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55457)。
原始数据已作为单独的补充文件随本修订稿提供。进一步的询问可直接联系通讯作者。
披露
作者在本研究中无利益冲突。
致谢
本工作得到国家自然科学基金 [资助号 82274435、82074223]；中央级重点专项：贵重中药资源可持续利用能力建设 [资助号 2060302]；以及 2022 年第五批全国中医药临床优秀人才研修项目 [资助号 2022178] 的支持。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| AutoDock Vina | https://vina.scripps.edu | 1.2.5 (SCR_011958) | 分子对接软件 |
| clusterProfiler (R 软件包) | https://bioconductor.org/packages/clusterProfiler | 4.10.0 (SCR_016884) | 用于富集分析的 R 软件包 |
| CytoHubba (Cytoscape 插件) | https://apps.cytoscape.org/apps/cytohubba | 0.1 | 用于枢纽基因识别的插件（Degree 算法） |
| Cytoscape | https://cytoscape.org | 3.10.1 (SCR_003032) | 网络可视化软件 |
| Discovery Studio | Dassault Systèmes BIOVIA | 2021 | 用于生成二维相互作用图的软件 |
| enrichplot (R 软件包) | https://bioconductor.org/packages/enrichplot | 1.22.0 (SCR_021165) | 用于富集结果可视化的 R 软件包 |
| GROMACS | https://www.gromacs.org | 2025.3 (SCR_014565) | 分子动力学模拟软件 |
| limma (R 软件包) | https://bioconductor.org/packages/limma | 3.58.1 (SCR_010943) | 用于差异表达分析的 R 软件包 |
| org.Hs.eg.db (R 软件包) | https://bioconductor.org/packages/org.Hs.eg.db | 3.18.0 (SCR_006442) | 用于人类基因标识符注释的 R 软件包 |
| PyMol | Schrödinger, Inc | 2.5.7 (SCR_000305) | 分子可视化软件 |
| QtGrace | https://sourceforge.net/projects/grace/ | 0.2.6 | 用于轨迹分析的绘图工具 |
| R (编程环境) | https://www.r-project.org | 4.3.1 (SCR_001905) | 统计计算软件 |
| STRING 数据库 | https://string-db.org | 12 (SCR_005223) | 蛋白质-蛋白质相互作用数据库 |
| venn (R 软件包) | https://cran.r-project.org/package=venn | 1.11 | 用于生成维恩图的 R 软件包 |
| WGCNA (R 软件包) | https://cran.r-project.org/package=WGCNA | 1.72 (SCR_003302) | 专用于加权共表达网络分析的 R 软件包 |
参考文献
- Zhao, Y., Chen, G. Y., Fang, M. Research trends of rheumatoid arthritis and depression from 2019 to 2023: a bibliometric analysis. J Multidiscip Healthc. 17, 4465-4474 (2024).
- Venetsanopoulou, A. I., Alamanos, Y., Voulgari, P. V., Drosos, A. A. Epidemiology of rheumatoid arthritis: genetic and environmental influences. Expert Rev Clin Immunol. 18 (9), 923-931 (2022).
- Beidelschies, M., et al. Polycyclic aromatic hydrocarbons and risk of rheumatoid arthritis: a cross-sectional analysis of the National Health and Nutrition Examination Survey, 2007 - 2016. BMJ Open. 13 (5), e071514(2023).
- Xi, X., et al. Polycyclic aromatic hydrocarbons affect rheumatoid arthritis pathogenesis via aryl hydrocarbon receptor. Front Immunol. 13, 797815(2022).
- Xia, M., et al. Vehicular exhaust particles promote allergic airway inflammation through an aryl hydrocarbon receptor-notch signaling cascade. J Allergy Clin Immunol. 136 (2), 441-453 (2015).
- Sun, L., et al. Particulate matter of 2.5 µm or less in diameter disturbs the balance of TH17/regulatory T cells by targeting glutamate oxaloacetate transaminase 1 and hypoxia-inducible factor 1α in an asthma model. J Allergy Clin Immunol. 145 (1), 402-414 (2020).
- Hew, K. M., et al. Childhood exposure to ambient polycyclic aromatic hydrocarbons is linked to epigenetic modifications and impaired systemic immunity in T cells. Clin Exp Allergy. 45 (1), 238-248 (2015).
- Podechard, N., et al. Interleukin-8 induction by the environmental contaminant benzo(a)pyrene is aryl hydrocarbon receptor-dependent and leads to lung inflammation. Toxicol Lett. 177 (2), 130-137 (2008).
- Lee, J., et al. A role for benzo[a]pyrene and Slug in invasive properties of fibroblast-like synoviocytes in rheumatoid arthritis: a potential molecular link between smoking and radiographic progression. Joint Bone Spine. 80 (6), 621-625 (2013).
- Iqbal, J., et al. Smoke carcinogens cause bone loss through the aryl hydrocarbon receptor and induction of Cyp1 enzymes. Proc Natl Acad Sci USA. 110 (27), 11115-11120 (2013).
- Chen, G. Y., et al. Mechanisms of total glucosides of paeony in alleviating methotrexate-induced liver injury. Drug Des Devel Ther. 19, 3407-3423 (2025).
- Chen, G. Y., et al. Prediction of Rhizoma Drynariae targets in the treatment of osteoarthritis based on network pharmacology and experimental verification. Evid Based Complement Alternat Med. 2021, 5233462(2021).
- Chen, G., Yan, Z., Wang, Y., Tao, Q. Rheumatoid Arthritis and Fibromyalgia Syndrome: A Bibliometric and Bioinformatics Perspective on Comorbidity Research. J Multidiscip Healthc. 18, 6811-6827 (2025).
- Kim, S., et al. PubChem 2023 update. Nucleic Acids Res. 51 (D1), D1373-D1380 (2023).
- Zdrazil, B., et al. The ChEMBL database in 2023: a drug discovery platform spanning multiple bioactivity data types and time periods. Nucleic Acids Res. 52 (D1), D1180-D1192 (2024).
- Keiser, M. J., et al. Relating protein pharmacology by ligand chemistry. Nat Biotechnol. 25 (2), 197-206 (2007).
- Liu, X., et al. PharmMapper server: a web server for potential drug target identification using pharmacophore mapping approach. Nucleic Acids Res. 38 (2), W609-W614 (2010).
- Barrett, T., et al. NCBI GEO: archive for functional genomics data sets—update. Nucleic Acids Res. 41 (D1), D991-D995 (2013).
- Langfelder, P., Horvath, S. Fast R functions for robust correlations and hierarchical clustering. J Stat Softw. 46 (11), 1-17 (2012).
- Szklarczyk, D., et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 51 (D1), D638-D646 (2023).
- wwPDB consortium. Protein Data Bank: the single global archive for 3D macromolecular structure data. Nucleic Acids Res. 47 (D1), D520-D528 (2019).
- Gentner, N. J., Weber, L. P. Intranasal benzo[a]pyrene alters circadian blood pressure patterns and causes lung inflammation in rats. Arch Toxicol. 85 (4), 337-346 (2011).
- Ordóñez, D., et al. Cell-mediated cytotoxicity in Lyme arthritis. Arthritis Rheumatol. 75 (5), 782-793 (2023).
- Zheng, Y., et al. Role of the granzyme family in rheumatoid arthritis: current insights and future perspectives. Front Immunol. 14, 1137918(2023).
- Wang, H., et al. Triple knockdown of CD11a, CD49d, and PSGL1 in T cells reduces CAR-T cell toxicity but preserves activity against solid tumors in mice. Sci Transl Med. 17 (782), eadl6432(2025).
- So, T., et al. Antigen-independent signalosome of CARMA1, PKCθ, and TNF receptor-associated factor 2 (TRAF2) determines NF-κB signaling in T cells. Proc Natl Acad Sci USA. 108 (7), 2903-2908 (2011).
- Genheden, S., Ryde, U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opi Drug Discov. 10 (5), 449-461 (2015).
重印与许可
申请许可以重复使用本 JoVE 文章的文本或图表
申请许可