类风湿关节炎外周血单核细胞circRNA-miRNA-mRNA网络构建及关键基因预测

郭锦晨，刘健，王茎，周巧，黄旦

（1.安徽中医药大学新安医学教育部重点实验室，安徽 合肥 230038；2.安徽中医药大学第一附属医院风湿病科，安徽 合肥 230031）

类风湿关节炎（rheumatoid arthritis，RA）是一种病因尚未完全阐明以滑膜炎和血管炎为基本病理改变的慢性全身性炎症性自身免疫性疾病[1]，其病理机制尚未完全阐明，暂无根治之药，因此探索可靠、准确的生物学标志物对其早期诊断及疗效评估有重大意义。研究证实[2]，单核细胞/巨噬细胞可通过分泌肿瘤坏死因子-α（TNF-α）、白细胞介素1β（IL-1β）、白细胞介素6（IL-6）等炎性细胞因子来维持RA 患者的炎性环境，还可引起大量的免疫细胞迁移至关节炎，最终导致关节发生不可逆的损伤。研究发现[3]，环状RNA（circular RNA，circRNAs）在RA 患者外周血单核细胞中存在差异表达，部分学者认为可能成为诊断RA 的新型生物标志物，微小核糖核酸（microRNA，miRNA）在免疫细胞分化、成熟、稳态调节方面发挥着重要的作用，部分miRNAs 已证实与T、B淋巴细胞的激活、炎症因子和趋化因子的产生等自身免疫过程相关。circRNA 富含miRNA 的结合位点，大量研究表明可以起到海绵作用作为竞争内源ceRNAs，通过其miRNA 反应元件与miRNA 竞争性结合，然后间接调节miRNA 靶向信使RNA（mRNA）的表达，与RA 的发生发展密切相关[4]。本研究通过对基因表达综合数据库（Gene Expression Omnibus，GEO）中RA 患者和正常对照者的外周血单核细胞测序数据进行分析，研究circRNA、miRNA 及mRNA 差异表达谱，并构建内源性竞争RNA（ceRNA）网络，为探索RA 诊断和治疗的新靶点提供理论依据。

1 资料与方法

1.1 资料来源 在GEO 据库（http://www.ncbi.nlm.nih.gov/geo）下载RA 患者外周血单核细胞相关数据集：circRNA数据集（GSE189338）、miRNA数据集（GSE124373）和mRNA 数据集（GSE55457）。采用“Limma”软件包校正微阵列数据表达水平，采用“edgeR”软件包对RA 组与正常对照组样本进行差异表达分析，采用“pheatmap”软件包进行热图绘制，LogFC 过滤阈值为1，校正后P（false discoveryrate，FDR）＜0.05 设定为筛选标准，如果表达数值很大，需要对数据取Log2。

1.2 circRNA-miRNA-mRNA 网络的构建 根据组织样本信息，将circRNA 数据集GSE189338、miRNA数据集GSE124373 和mRNA 数据集GSE55457 分别分为RA 组和正常对照组，在3 个数据集中分别提取差异表达的circRNA、miRNA 和mRNA。将RA患者差异表达的circRNA 通过CSCD 数据库（https://gb.whu.edu.cn/cscd）数据库预测相应的靶标miRNA，结合GSE124373 中筛选得到的差异表达miRNA 进一步筛选目标miRNA。然后，利用miRDB（http://www.mirdb.org）和TargetScan（http://www.tar getscan.org/vert_72）数据库预测RA 患者目标miRNA对应的mRNA，再与GSE55457 中RA 患者差异表达的mRNA 进一步取交集以识别目标mRNA。采用“venn”软件包对各差异表达的miRNA、mRNA 与靶标取交集，基于ceRNA 理论通过Cytoscape 3.7.1 软件构建ceRNA 调控网络。

1.3 ceRNA 网络节点热图和箱线图分析 采用“pheatmap”软件包对RA 患者外周血单核细胞ceRNA 网络节点中circRNA、miRNA 和mRNA 进行热图绘制。“reshape2”包是由Hadley Wickham 开发的一个R 包，利用reshape2 包可以融化数据进行重新整合，并对长宽数据格式进行相互转换，最后使用“ggpubr”包绘制circRNA、miRNA 和mRNA 箱线图。

1.4 网络中mRNA 的GO 和KEGG 功能富集分析为进一步分析mRNA 的功能和主要作用通路，将上述筛选得到的mRNA 输入DAVID 数据库（https://david.ncifcrf.gov/），选定物种为“homo sapiens”，设定阈值P＜0.05，将RA 关键基因从生物过程（biological process，BP）、细胞组成（cellular component，CC）、分子功能（molecular function，MF）3 方面进行基因本体功能（GO）富集分析。利用R 语言“clusterProfiler”包对RA 患者关键靶点进行KEGG（https://www.kegg.jp/）富集分析，并使用“ggplot2”包绘制气泡图。

2 结果

2.1 差异表达的circRNA、miRNA 和mRNA 在GSE189338 数据集发现13 个上调circRNA 和45个下调circRNA（图1A）；在GSE124373 数据集发现24 个上调miRNA 和11 个下调miRNA（图1B）；在GSE55457 数据集发现112 个上调mRNA 和38 个下调mRNA（图1C）。从CSCD 数据库中检索相关circRNA 数据，24 个circRNA 的基本结构模式见图2，24 个circRNA 的基本信息，包括位置、基因组长度、靶向基因名等见表1。

图1 RA 相关差异表达circRNA、miRNA、mRNA 热图

图1 RA 相关差异表达circRNA、miRNA、mRNA 热图（续）

表1 24 个差异表达circRNA 的基本特征

图2 筛选出的24 个circRNA 基本模式图

2.2 circRNA-miRNA-mRNA 网络的构建 通过CSCD 数据库对circRNA 进行靶向预测，筛选出1378 个靶向miRNA。将这1378 个靶标miRNA 与miRNA 数据集GSE124373 得到的35 个差异miRNA取交集，筛选出15 个目标miRNA（图3A）。再利用TargetScan 和miRDB 数据库对这15 个目标miRNA进行靶向预测，筛选出7972 个靶标mRNA，将这7972 个靶标mRNA 与mRNA 数据集GSE55457 得到的150 个差异mRNA 取交集，筛选出60 个目标mRNA（图3B）。根据ceRNA 网络关联数据并运用Cytoscape 3.7.1 软件，最终筛选出10 个circRNA、8个miRNA 和38 个mRNA 构建circRNA-miRNA-mRNA可视化网络，见图4。

图3 差异表达miRNA 及mRNA 筛选

图4 可视化ceRNA 网络

2.3 ceRNA 网络节点热图和箱线图 通过对ceRNA 网络节点进行热图和箱线图分析，发现6 个上调circRNA，4 个下调circRNA，其中，hsa_circ_0009012、hsa_circ_0003695、hsa_circ_0001589、hsa_circ_0092-125、hsa_circ_0074598 在RA 患者中高表达，hsa_circ_0058230、hsa_circ_0009127、hsa_circ_0067-492、hsa_circ_0001543 在RA 患者中低表达；发现4个上调miRNA，4 个下调miRNA，其中hsa-miR-6165、hsa-miR-98-5p、hsa-miR-619-5p、hsa-miR-186-5p 在RA 患者中高表达，hsa-miR-4516、hsamiR-3613-3p、hsa-miR-766-3p、hsa-miR-6737-3p 在RA 患者中低表达；发现27 个上调mRNA，11 个下调mRNA，其中，ADAM28、TREML2、PDCD1LG2、CD38、CD79A、RASGRP1、DAZL、ST8SIA1、E2F8、KIF23、TMEM156、BARD1、MS4A1、PNOC、LAMP3、PIM2、SEMA4D 在RA 患者中高表达，LGALSL、CYTH3、ADM2、KCNK3、TRHDE、MXD1、NR4A3 在RA 患者中低表达，见图5、图6。

图5 ceRNA 网络节点circRNA、miRNA、mRNA 热图分析

图5 ceRNA 网络节点circRNA、miRNA、mRNA 热图分析（续）

图6 ceRNA 网络节点circRNA、miRNA、mRNA 箱线图分析

2.4 靶点通路分析 将上述筛选出的38 个差异mRNA通过DAVID 数据库进行GO 功能富集分析，共得到GO 条目289 个，在BP 中，基因主要富集于B cell activation、B cell differentiation、lymphocyte differentiation、mononuclear cell differentiation、B cell prolifera tion 等252 个条目中；在CC 中，基因主要富集在external side of plasma membrane、RNA polymerase II transcription repressor complex、lamellar body、pericentriolar material、ruffle 等5 个条目中；在MF 中，基因主要富集在fibronectin binding、cytokine binding、receptor ligand activity、signaling receptor activator activity 等32 个条目中；各类别排名前10 的GO 条目绘制气泡图，见图7。KEGG 通路富集筛选得到11 条信号通路（P＜0.05），见图8。

图7 差异表达mRNA GO 富集分析

图8 差异表达mRNA KEGG 富集分析

3 讨论

RA 是一种病因和发病机制尚未明确的慢性炎症性自身免疫性疾病，免疫炎症反应可破坏全身多关节软骨和骨骼，氧化应激反应可进一步加重炎症组织损伤，最终导致患者丧失劳动力、甚至致残[5]。CircRNA 是一类有较高的稳定性和进化保守性的内源性表达的具有闭合环状结构非编码RNA，越来越多的证据表明[6]，包括circRNAs、miRNA 和mRNA在内的调控网络在RA 的发病过程中起着重要的调节作用，参与自身免疫和炎症反应的调节。如Li B等[7]研究发现，hsa_circ_0001859 可以抑制miR-204/211 活性而上调激活转录因子2（activating transcription factor 2，ATF2），从而促进RA 的炎症发展。有专家对RA 患者PBMCs 进行高通量测序发现，在RA中很多差异表达的circRNA 可能是RA 的潜在诊断生物标志物，并且在RA 的发病过程中可能起调控作用[8]。因此，本研究基于ceRNA 理论构建circRNA-miRNA-mRNA 网络筛选参与调控RA 的潜在mRNA，为RA 临床治疗及诊断提供新靶标。

本研究中首先从GEO 数据库获得RA 相关circRNA、miRNA 和mRNA 测序数据集进行差异分析，筛选出差异表达circRNA、miRNA、mRNA 分别有58 个、35 个、150 个。基于ceRNA 理论进一步构建circRNA-miRNA-mRNA 网络，筛选出10 个circRNA、8 个miRNA 和38 个mRNA 组成的ceRNA 调控网络。通过对ceRNA 网络节点进行箱线图分析，发现hsa_circ_0009012、hsa_circ_0003695 等circRNA在RA患者中高表达，hsa_circ_0058230、hsa_circ_00091273 等circRNA 在RA 患者中低表达，hsa-miR-6165、hsa-miR-98-5p 等miRNA 在RA患者中高表达，hsa-miR-4516、hsa-miR-3613-3p 等miRNA 在RA 患者中低表达，ADAM28、TREML2、PDCD1LG2、CD38、CD79A、RASGRP1、DAZL 等mRNA在RA 患者中高表达，LGALSL、CYTH3、ADM2、KC NK3、TRHDE 等mRNA 在RA 患者中低表达。通过对差异基因集进行GO 和KEGG 富集分析，发现分子功能主要富集在纤维蛋白结合、细胞因子结合、受体配体活性、信号受体激活剂活性等，生物过程主要富集在B 细胞活化、分化及增殖、淋巴细胞分化、单核细胞分化等。KEGG 通路主要富集在造血细胞谱系、原发性免疫缺陷、B 细胞受体信号通路、Toll 样受体信号通路、T 细胞受体信号通路等。

CD38 一种广泛分布于不成熟的造血细胞和活化的T、B 淋巴细胞表面的Ⅱ型穿膜糖蛋白，参与T细胞的共刺激，调控B 细胞的成熟、分化与凋亡，在细胞活化、信号转导、介导细胞因子产生等过程中发挥重要的作用[9]。研究发现[10]，RA 患者外周血淋巴细胞中CD38+细胞表达显著升高，CD38 能够调节IL-1α 和IL-1β 的分泌水平介导RA 免疫紊乱。CD79A 是参与B 淋巴细胞信号转导的基因，富集于T 细胞受体复合物，可与B 细胞受体结合并在识别抗原后产生信号[11]。RAS 鸟嘌呤释放蛋白1（RAS guangreleasing protein 1，RASGRP1）是一种作用于MAPKs 和T 细胞受体等信号通路的RAS 鸟氨酸交换因子，可结合钙离子参与细胞内外信号传递[12]。E2F8 通过介导CyclinD1 在调节细胞周期、增殖、分化及凋亡和DNA 损伤修复过程中发挥重要作用，参与Notch 与NF－κB 信号通路的调控[13]。溶酶体相关膜蛋白3（LAMP3）是自噬和促进脂质生成过程中不可缺少的一个分子，在特定的细胞核分化阶段表达，在脂代谢、细胞降解途径起重要作用[14]。PIM2 一种对RA 有抑炎作用的丝/苏氨酸蛋白激酶，在RA 中发现外源性过表达PIM2 能够活化哺乳动物雷帕霉素靶蛋白1（mam malian target of rapamycin 1，mTORC1），从而阻断NF-κB 和COX-2 等信号通路介导的炎症反应[15]。信号素4D（Semaphorin 4D，Sema4D）是一个跨膜型同源二聚体蛋白质，在神经细胞导向、免疫调节、骨代谢及血管再生等过程中发挥重要作用，免疫系统所有T 细胞、活化B 细胞均表达Sema4D，血清Sema4D 过表达可能与RA 患者肺间质病变（interstitial lung disease，ILD）的发生存在一定联系[16]。NR4A3 属于核激素受体（NHRs）超家族中一种具有转录激活与抑制的双向转录调控因子，能够参与调节细胞生物学行为、炎症反应、应激反应、能量代谢、血内皮细胞损伤等[17]。TLRs 是一个对适应性免疫、固有免疫都具有重要作用的单跨膜受体家族，Toll 样受体信号通路由髓样分化因子（MyD88）依赖与非依赖性通路组成，均经NF-κB 诱导炎性细胞因子而产生[18]。B 细胞受体信号通路在RA 中起着重要作用，可分泌TNF-α、IL-6 等细胞因子、CXCL13 等趋化因子引发关节炎症，还可充当抗原递呈细胞活化抗原特异性T 细胞，导致关节组织炎性损伤。T 细胞受体信号通路属于免疫调节的重要通路之一，活化的CD4+T 细胞与人类白细胞抗原（HLA）、组织相容性复合体（MHC-Ⅱ）等分子相互作用可诱导PI3K-Akt 信号通路，导致CD4+T 细胞的成熟，从而激活CD8+T 细胞抗原促进炎症发生。此外，CD4+T 细胞的凋亡异常还会导致RA 患者大量分泌致炎因子，加重炎症反应，促进疾病活动。

综上所述，RA 是一种由多种因素导致的免疫系统疾病，发病机制复杂。本研究筛选出RA 患者外周血单核细胞中差异表达的lncRNA、circRNA和mRNA，构建了ceRNA 调控网络，分析了差异表达mRNA 的生物学功能和调控通路情况。