陈雪芹，李灵兰，付国涛，韩莉欣，毋瑞朋，熊咏民，史传道，刘启玲，张荣强
1陕西中医药大学公共卫生学院流行病与卫生统计学教研室，陕西咸阳 712046 2陕西省卫生健康委员会疾病控制处，西安 710003 3西安交通大学公共卫生学院地方病研究所，西安 710061 4西藏民族大学医学公共卫生与预防医学教研室，陕西咸阳 712082
类风湿性关节炎(rheumatoid arthritis，RA)是一种以持续滑膜炎症、系统性炎症和关节不可逆损伤为主要临床表现的关节炎性疾病，是遗传及环境改变等多因素综合的结果[1]，可造成后期关节畸形和运动障碍，严重影响患者的生活质量[2]。我国RA患者约500万，其导致的致残率呈明显升高趋势[3]。然而，RA的发病机制仍不明确，积极探索其病因机制对RA防治具有重要的现实意义。自噬是一种以细胞“自食”的凋亡方式维持细胞稳态的分解代谢过程[4]。骨骼生长板中的软骨细胞位于相对缺乏氧气和营养的内部环境中，自噬与凋亡呈现负相关调控关系，自噬保护肥大的软骨细胞免于过早凋亡和骨异常替代，增加了软骨细胞的增殖和活性[5]。截至投稿前，已经发现了许多与RA相关的自噬信号通路，如磷脂酰肌醇3-激酶/蛋白激酶B/哺乳动物雷帕霉素靶蛋白、苄氯素1、自噬相关基因5(autophagy-related gene 5，ATG5)的调节、丝裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)/细胞外信号调节激酶、TP53与免疫信号调控和非编码RNA与自噬之间的关系[6]。RA引起的破骨功能异常激活促进肿瘤坏死因子-α(tumor necrosis factor alpha，TNF-α)升高，介导ATG7和自噬蛋白苄氯素1表达，促使破骨细胞直接分泌溶酶体至细胞膜，诱导破骨细胞相关基因的mRNA表达，有利于骨基质体外吸收和骨融合[1]。本研究纳入GEO公共数据库(https：//www.ncbi.nlm.nih.gov/geo)的数据集GSE55457和GSE55235，筛选出RA滑膜组织的差异表达基因(differentially expressed gene，DEG)，旨在通过鉴定差异表达自噬相关基因(differentially expressed autophagy-related gene，DE-ATG)的表达及转录因子(transcription factor，TF)-基因-miRNA模型的建立，并分析DE-ATG在RA的表达模式和诊断鉴别能力，RA自噬可能参与的信号通路和靶向miRNA的凋亡调控机制，RA与环境化学物、药物小分子之间的关系，为RA临床治疗提供诊疗依据与治疗靶点。
数据来源及预处理在美国国家生物技术信息中心的GEO“DataSets”版块以“Rheumatoid arthritis”为关键词，“Human sapiens”为检索类型，检索出两组RA滑膜组织的基因表达数据集GSE55235和GSE55457，共包含23个RA患者滑膜组织和20个健康者滑膜组织。本研究纳入的数据集属于RNA测序数据，多个探针定位到同一基因时，计算其平均数，删除无用或空白数据；
采用“background Correct”R包进行背景校正，使用“normalize”R包对数据归一化，使用“remove Batch Effect”包去除批次效应，确保每个样本的表达分布是同质的。人类自噬数据库(Human Autophagy Database，HADb)(http：//autophagy.lu/clustering/index.html)获得232个ATG。
差异表达基因的鉴定与筛选“limma”R包筛选RA组相比于健康组的差异表达基因，筛选标准为|logFC|>1，错误发现率<0.05。热图、火山图分别由“pheatmap”和“volcano”R计算和描绘，在Rstudio(版本4.0.3)完成分析。两组间的比较使用Wilcoxon秩和检验，多组比较采用Kruskal-WallisH检验，P<0.05为差异有统计学意义。
差异表达基因的富集分析上调和下调基因的基因本体功能(gene ontology，GO)富集分析和京都基因与基因组百科全书通路(Kyoto Encyclopedia of Genes and Genomes，KEGG)分析由R语言“ggplot2”、“ClusterProfiler”和“GOplot”软件包分析和绘制。GO功能注释由细胞成分、生物过程和分子功能组成，KEGG分析侧重于代谢途径和分子机制。
DE-ATG的筛选及表达DEG和ATG的交集DE-ATG。首先提取GSE55457和GSE55235的DE-ATG的mRNA表达数据，采用F检验进行方差齐性检验，独立样本t检验进行组间比较，P<0.05为差异有统计学意义。
蛋白互作网络分析STRING 11.5数据库(http：//www.genemania.org)用来分析DE-ATG之间的潜在交互作用[7]，信度设置为0.7，其余参数为默认设置，将结果导入Cytoscape软件3.8.0中进行可视化，建立蛋白交互作用网络。GeneMANIA数据库(http：//www.genemania.org)可以通过大量基因组学和蛋白质学数据查找与DE-ATG功能相似的基因[8]，设置阈值错误发现率<0.01，P<0.05，根据查询的预测值对每个基因组数据集进行加权，探索DE-ATG与其功能相似基因之间的关系。
TF-基因-miRNA共表达网络建立采用Network Analyst在线分析工具(https：//www.networkanalyst.ca)的miRTarBase、TarBase和miRecords预测DE-ATG和miRNA之间的调控关系。NetworkAnalyst中的ENCODE、JASPAR和ChEA数据库用于预测DE-ATG的TF。TF-基因-miRNA共表达网络的结果采用Cytoscape 3.8.0实现图形可视化。
DE-ATG与环境化学物的相互作用比较毒理学数据库(http：//ctdbase.org)是一个关于环境影响疾病的病因和分子机制的可检验假设的资源数据库[9]。将DE-ATG上传至比较毒理学数据库，探索DE-ATG与环境化学小分子之间的相互作用关系。
DE-ATG对RA诊断能力的受试者工作曲线采用受试者工作曲线分析评估DE-ATG鉴别诊断RA的敏感性和特异性，曲线下面积(area under the curve，AUC)>0.75提示DE-ATG对RA具有良好的预测准确性；
本研究按照样本量的7∶3比例将原始数据分为训练集和验证集，采用留一法进行交叉验证。
DE-ATG的靶向药物预测DrugBank(www.drugbank.ca)是一个由美国食物药品监督管理局批准程序的具有实验药物数据和综合药物靶点的数据库[10]。在DrugBank数据库中分析DE-ATG，以确定潜在治疗自噬相关的RA药物小分子。
数据预处理结果在背景校正、数据归一化和去除批次间差异后，数据稳定性良好，可保证后续分析结果的准确可靠(图1A)。
差异表达基因的筛选与正常组织相比，在RA滑膜组织共鉴定出485个DEG。其中上调277个、下调208个(图1B、1C)。
差异表达基因的富集分析结果对DEG进行GO和KEGG富集分析，GO富集分析结果显示，上调DEG主要参与T细胞激活、细胞间黏附的正向调控、淋巴细胞和白细胞增殖过程，下调DEG主要参与平滑肌细胞增殖和激素水平的调节等。KEGG富集分析结果显示，上调DEG主要涉及原发性免疫缺陷、RA和趋化因子等信号通路，下调DEG主要涉及酪氨酸代谢、破骨细胞分化和受体-因子信号传导等信号通路(图2)。
KEGG：京都基因与基因组百科全书通路；
GO：基因本体功能
DE-ATG的验证及表达在差异表达基因中，共有11个DE-ATG。和健康对照组相比，APOL1、CCR2、CXCR4、TNFSF10的mRNA水平在RA滑膜组织中明显升高，而CDKN1A、FOS、GABARAPL1、MAPK8、MYC、PPP1R15A和VEGFA的mRNA水平在RA中明显降低(图3)。
RA：类风湿性关节炎；
DE-ATG：差异表达自噬相关基因
蛋白互作网络分析结果STRING数据库揭示了DE-ATG预测的密集蛋白调控网络，其中包含31个节点和179条边(P<0.001)(图4A)。GeneMANIA揭示了DE-ATG与功能相似的20个基因建立了12 479条密切的关系网，其中包括93.65%的共表达作用，4.21%的物理相互作用，1.27%的共定位作用，0.75%的预测作用，0.07%的蛋白质区域共享和0.05%的通路分析。
这11个DE-ATG和前20个相似基因主要与细胞应激和RNA聚合酶Ⅱ对转录的作用有关(图4B)。
A.STRING建立的蛋白互作网络(红色椭圆是上调的DE-ATG，绿色椭圆是下调的DE-ATG，紫色椭圆是DE-ATG靶向调控的核心蛋白)；
B.GeneMANIA数据库预测的蛋白(结点代表蛋白的生物学功能，网络是蛋白间相互作用关系)
TF-基因-miRNA共表达网络构建结果通过NetworkAnalyst的3个数据库(ENCODE、JASPAR和ChEA)获得380个TF(图5A)，在预测miRNA的数据库中发现876个节点和1341条边(图5B)，在TF-基因-miRNA共表达网络中筛选到1245个节点和2264条边，并由Cytoscape的CytoHubba插件筛选根据度值排序的前200个TF和miRNA网络(图6)。其中，前10个关键TF(PPATG、GATA2、TP53、SOX2、RUNX1、HNF4A、FOXC1、TP63、E2F1、CREB1)和前10个关键mi-RNA(hsa-mir-124-3p、hsa-mir-1-3p、hsa-mir-17-5p、hsa-mir-16-5p、hsa-mir-155-5p、hsa-mir-335-5p、hsa-mir-34a-5p、hsa-mir-20a-5p、hsa-mir-107、hsa-mir-497-5p)的鉴定有助于阐明miRNA和TF在RA发展中的影响。
TF：转录因子；
红色菱形代表上调的DE-ATG，绿色菱形代表下调的DE-ATG，蓝色椭圆代表靶向预测的TF，黄色椭圆代表靶向预测的miRNA
红色菱形代表上调的DE-ATG，绿色菱形代表下调的DE-ATG，蓝色椭圆代表靶向预测的TF，黄色椭圆代表靶向预测的miRNA
DE-ATG与环境化学物的相互作用结果通过比较毒理学数据库发现1387个基因-环境化学物产物和2657条联系。根据度值，排名前10的环境化学物分别是四氯二苯二氧芑、二氧化硅、氧气、硫酸镍、(+)-JQ1化合物环、环孢菌素、苯并(a)芘、锌、丙戊酸、曲格列酮(图7)。
红色菱形代表上调的DE-ATG，绿色菱形代表下调的DE-ATG，紫色三角形代表DE-ATG相互作用的环境化学物
DE-ATG对RA诊断能力的ROC曲线ROC曲线显示，11个DE-ATG在验证集和训练集中AUC值均大于0.85，且均具有统计学意义(P<0.05)，提示本部分结果稳定可靠，表明11个DE-ATG对RA均具有良好的预测和指示能力(图8)。
图8 11个DE-ATG的受试者工作特征曲线
DE-ATG的靶向药物预测结果为了预测并开发用于治疗由ATG调控的RA，DrugBank数据库过滤到食物药品监督管理局认证的2860个相互作用的药物小分子，其中有氟米龙、轮状病毒疫苗、睾酮内酯、阿霉素、聚二甲基硅氧烷和青蒿琥酯等。
全球每年10万RA病例中有5～50例新发病例，其中50%的RA发生风险与遗传因素有关[11]，影响全球约2000万人的生活质量[12]。免疫和遗传因素都会导致免疫细胞激活自身抗原，促使自身抗体水平的持续升高，例如抗环瓜氨酸抗体和类风湿因子。在RA中发现瓜氨酸化的自蛋白的抗体反应，瓜氨酸肽的产生与自噬有关[13]，促使T细胞、B细胞、巨噬细胞、滑膜成纤维细胞、软骨细胞和破骨细胞异常的信号传导引起不同炎症介质的释放[14]，全身持续慢性炎症反应。
RA滑膜成纤维细胞的侵袭性和侵袭性表型的信号通路表达发生改变，加剧细胞凋亡和自噬之间的负相关程度，引起进行性骨和软骨破坏[15]。本研究深入挖掘GEO数据库的RA患者滑膜组织的RNA-seq数据，鉴定出差异表达上调基因277个、下调基因208个。这485个差异表达基因与T细胞活化、淋巴细胞增生分化、激素反应相关，参与肾脏血管内血脂调节、平滑肌增殖和代谢的调控反应，富集在某些特定信号通路，如核因子κB(nuclear factor κB，NF-κB)、原始免疫缺陷、TNF、破骨细胞分化和白细胞介素(interleukin，IL)-17等信号通路。Th17细胞产生的细胞因子，如TNF、IL-22、IL-21和粒细胞-巨噬细胞集落刺激因子，与疾病活动评分、C反应蛋白和抗硝酸蛋白抗体呈正相关[16]，刺激滑膜成纤维细胞和巨噬细胞产生大量的促炎因子，如IL-1、IL-6、TNF和前列腺素E2，引起软骨破坏和细胞外基质蛋白的降解，加剧滑膜炎症。激活的滑膜成纤维细胞通过IL-8、IL-17(RA早期高度表达)等炎症因子使中性粒细胞、单核细胞和巨噬细胞聚集在关节炎部位，Th17诱导滑膜基质和淋巴细胞分泌粒细胞-巨噬细胞集落刺激因子，加速骨基质的破坏[17]。IL-17、IL-1和TNF升高促使IL-6、基质金属蛋白酶(matrix metallo-proteinase，MMP)和Ⅰ型胶原C端肽的表达在关节滑膜上调、IL-17诱导成骨细胞分化、滑膜中核因子-κB 受体激活剂配体的表达和激活蛋白-1信号，通过调节MMP3、MMP13和血小板反应蛋白解整合素金属肽酶4的表达对软骨细胞造成损害[12，18]，促进破骨细胞的生成[17]。此外，IL-17介导的线粒体功能障碍还可以通过激活自噬加速RA滑膜成纤维细胞的凋亡[19]。在HADb数据库映射到RA的11种DE-ATG中，APOL1、CCR2、CXCR4和TNFSF10在RA中的表达显著升高，CDKN1A、FOS、GABARAPL1、MAPK8、MYC、PPP1R15A和VEGFA在RA中的表达显著降低。STRING预测到的ATG7在基因敲减后出现TNF介导的骨破坏受损，部分原因是通过抑制自噬减少了IL-6和IL-1的产生[20]。同时发现自噬相关蛋白视神经素通过调节NF-κB和干扰素-β信号传导负向调节破骨细胞生成[21]。
GeneMANIA表明DE-ATG与RNA聚合酶Ⅱ对转录的调控关系、蛋白激酶和MAPK信号级联反应具有关联性。此外，TF-基因-miRNA共表达网络显示DE-ATG与hsa-mir-155-5p有较强的调控关系。Jiang等[22]在TNF-α刺激的人滑膜肉瘤细胞和滑膜成纤维细胞中发现hsa-miR-147b-3p、hsa-miR-155-5p功能的丧失会下调TNF-α、IL-6、MMP3和MMP13的表达，通过NF-κB和JAK/STAT通路在类风湿性关节炎中发挥主要病理作用[12]。受DE-ATG调控的TF与RNA聚合酶的结合，抑制或促进DE-ATG在RA中的表达。例如在小鼠体内注射Runx1 mRNA，能够减少骨赘形成和软骨破坏，而Runx2作为成骨细胞增殖的TF[23]在促进软骨细胞肥大上发挥关键作用，直接调节MMP13的表达[12，24]。SoxC TF家族对软骨内骨化和关节软骨细胞降解有独特的方式，Sox4和Sox11经过视黄酸和IL-1处理后在关节软骨细胞中上调，直接与启动子结合上调Adamts4和Adamts5，使MMP13的表达上调，加速骨基质降解，破坏关节软骨[12，24-25]。
比较毒理学数据库对DE-ATG与环境因子相互作用的分析表明，CDKN1A相关的环境因素数量最高，其次是VEGFA、MYC和FOS。芳基烃受体(Aryl hydrocarbon receptor，Ahr)作为环境污染物的受体，介导二氧芑的毒性和生化作用[26]，在滑膜组织中高度表达，Ahr激活使Wnt/β-连环蛋白信号传导下调[27]，抑制胶原诱导性关节炎中成骨细胞分化和增殖[28]。四氯二苯二氧芑通过激活Ahr途径，抑制β-连环蛋白表达导致骨侵蚀，与RA早期成骨细胞的功能有密切关系[29]。Ilar等[30]基于人群的病例对照研究表明，易感个体暴露于二氧化硅等空气传播因子，瓜氨酸化触发抗原的翻译后修饰产生瓜氨酸化肽抗体，诱导类风湿因子(抗IgG)的产生，促进RA的发展。ROC分析显示AUC值均大于0.75，具有较好的预测效能，验证了11种DE-ATG参与到RA的炎症机制调控，其中CDKN1A、GABARAPL1、MAPK8、MYC、PPP1R15A和TNFSF10的AUC>0.9，提示可能是新的RA生物标志物候选物。DrugBank提示氟米龙、睾酮内酯和聚二甲基硅氧烷与DE-ATG靶点具有良好的结合活性。
综上，本研究鉴定出的11个DE-ATG创造了一个自噬调控RA滑膜表达的环境，这可能促进软骨细胞凋亡和分解，加深RA和自噬关系的理解。DE-ATG靶向TF和miRNA，加速RA骨基质降解，奠定了RA关节软骨破坏及炎症浸润基础，在一定程度上影响疾病的治疗和预后，对鉴别诊断RA具有一定的价值。11种新型生物标志物和药物的初筛为RA的诊断和治疗提供了新的思路和方法。
猜你喜欢滑膜软骨靶向新型抗肿瘤药物：靶向药物中老年保健(2022年1期)2022-08-17基于滑膜控制的船舶永磁同步推进电机直接转矩控制研究舰船科学技术(2022年10期)2022-06-17如何判断靶向治疗耐药保健医苑(2022年5期)2022-06-10SOX9在SD大鼠胚胎发育髁突软骨与胫骨生长板软骨中的时间表达研究昆明医科大学学报(2022年2期)2022-03-29高层建筑施工中的滑膜施工技术要点探讨建材发展导向(2021年14期)2021-08-23鞍区软骨黏液纤维瘤1例中国临床医学影像杂志(2019年5期)2019-08-27靶向治疗 精准施治 为“危改”开“良方”中华建设(2019年7期)2019-08-27髓外硬膜内软骨母细胞瘤1例中国临床医学影像杂志(2019年1期)2019-04-25靶向超声造影剂在冠心病中的应用医学研究杂志(2015年7期)2015-06-22关节镜在膝关节滑膜软骨瘤病诊治中的应用中国民族民间医药·下半月(2014年2期)2014-09-26