骨关节炎(Osteoarthritis,OA)的临床表现具有显著的异质性,其典型症状包括手、膝、髋及脊柱等关节的疼痛、功能障碍和结构畸形,这些症状不仅严重降低了患者的生活质量,还对社会经济造成了沉重的负担。近年来,研究揭示滑膜组织中Dkk3、OSM和ADAM12等基因的异常表达与OA的病理发展存在紧密联系
[1-3]。同时,众多研究者借助生物信息学的方法,构建了OA相关的分子互作网络,为揭示其发病机制提供了新的研究方向。
在表观遗传学领域,DNA甲基化作为调控基因表达的关键机制,参与多种生物学过程。研究者通过整合基因表达谱与DNA甲基化数据,发现特定基因的甲基化修饰模式与OA进展密切相关,这一发现为深入理解OA的分子机制研究提供了新的视角。
在分子机制的研究中,随着高通量测序技术的飞速发展,研究者发现circRNA凭借其独特的共价闭合环状结构,通过结合RNA结合蛋白、调控基因剪接及修饰亲本基因表达等方式参与多种病理过程
[4-6]。特别是,circRNA作为竞争性内源RNA,能够通过“分子海绵”机制调控miRNA功能,这一作用模式在动脉粥样硬化、朊病毒病和强直性肌营养不良等疾病中已得到验证
[4]。在OA研究中,circRNA芯片筛查显示关节软骨中存在71个差异表达的circRNA,其中功能实验证实circRNA-CER能够吸附miR-136促进MMP13表达上调,参与OA发病
[5]。此外,在肿瘤机制研究中,分子互作实验揭示高表达的CIRC-Amotl1通过与c-Myc相互作用促进肿瘤的恶性进展,为肿瘤治疗提供了新的潜在靶点
[6]。
大多数研究主要集中在骨关节炎中的单个基因事件及其相关机制上。事实上,骨关节炎的发病机制可能涉及众多分子和信号通路,分子和信号通路之间的异常串扰导致软骨合成代谢和分解代谢的失衡,最终导致骨性关节炎的发生和发展。因此,寻找可靠的生物标志物对OA的早期诊断和治疗具有重要意义。
本研究基于基因表达综合数据库(GEO,
https://www.ncbi.nlm.nih.gov/geo/)中的OA患者和正常对照组的滑膜组织微阵列数据进行了生物信息学分析,以确定OA患者滑膜组织的关键生物标志物。构建了circRNA-miRNA-mRNA调控网络,为骨关节炎的治疗提供了潜在的靶点和候选生物标志物。
1 材料和方法
从GEO数据库检索基因表达数据和DNA甲基化测序数据。DNA甲基化测序数据来自人的软骨组织数据集GSE63695(78个OA样本和19个正常对照组样本),表达数据来自人的滑膜组织数据集GSE55235、GSE55457中的OA和正常对照组样本各10例,并且这些表达数据来自同一个GPL平台。
1.1 差异表达基因的数据预处理与鉴定
首先,利用cbind函数将上述两个表达数据进行合并以便进行后续分析。使用removeBatchEffect函数去除批次效应,这有助于消除由不同实验批次引入的技术差异
[7]。通过制定适当的设计矩阵和批次信息,确保数据在进行差异分析之前是一致的和可比的。根据Affymetrix公司提供的GPL96平台注释文件将探针名称转换为基因符号,并删除缺失数据。使用R中的LIMMA包在合并的数据集中识别OA患者和对照组之间的DEGs。
1.2 甲基化数据处理与分析
CHAMP包适用于450 K和850 K的甲基化数据分析。首先,将数据导入R并进行质量控制,之后使用BMIQ对获得的原始数据进行归一化,并使用ComBat对批次效应进行校正,以确定差异甲基化位点(Differentially methylated position,DMP)。DMP被定义为带有|△β|>0.1和校正后的P<0.05的胞嘧啶-磷酸-鸟嘌呤(CpG)位点,含有差异甲基化位点的基因被鉴定为DMGs。
1.3 基因表达和DNA甲基化联合分析
将鉴定的DEGs与DMGs取交集,以获得受DNA甲基化状态改变影响的基因。疾病的发生和发展是极其复杂的,它是由多种调节因子共同作用的。所以,有必要探究DNA甲基化与基因表达之间的关联程度。相关系数是反映变量之间相关程度的统计指标,它的取值范围是[-1,1],当取值为0时表示不相关,当取值为[-1,0)时表示负相关,取值为(0,1]时表示正相关。
1.4 功能富集分析
富集分析使用clusterProfiler包来计算,通过KEGG和GO富集分析,研究差异基因最显著的富集途径和生物过程,统计显著性阈值设定为P<0.05。
1.5 PPI网络的构建
蛋白质-蛋白质相互作用网络(PPI)由蛋白质彼此之间相互作用构成,可以用来识别关键基因和重要模块。检索相互作用基因的搜索工具(STRING,
https://string-db.org)数据库是一个全球资源,提供了PPI的验证和预测信息。本文基于STRING数据库构建了受DNA甲基化影响的差异基因的PPI网络,使用Cytoscape中的CytoHubba来确定关键基因。
1.6 关键基因的验证
依据关键基因在训练集的表达数据,通过pROC包绘制ROC曲线来识别关键基因,曲线下面积AUC值越大表示诊断效果越好。
1.7 circRNA-miRNA-mRNA调控网络
最近的研究表明,大约80%的人类基因组被转录为非编码RNA(ncRNA),这些ncRNA已被证实能够调控多种重要的生物学功能,包括microRNA(miRNA)、长链非编码RNA(lncRNA)和环状RNA(circRNA)
[8]。作为一种新型的非编码RNA,circRNA作为其靶向miRNA的竞争性内源性RNA发挥作用,调控靶基因的表达
[9]。因此本研究也探索了潜在的circRNA-miRNA-DEGs的调控网络。circBase(
http://www.circbase.org/)预测了相交的差异基因的circRNA。针对交叉点基因的miRNAs从两个在线数据库miRDB和TargetScan获得,这两个数据库收集了实验验证的miRNA-genes相互作用。通过circBANK预测了circRNA和miRNA之间的海绵相互作用,利用Cytoscape可视化了circRNA-miRNA-mRNA调控网络。
2 结果与讨论
2.1 基因表达分析
在GSE55235和GSE55457中OA患者与对照组相比,
P<0.05且|log
2FC|>1被认为是DEGs筛查的临界值。最终得到的差异基因共576个,其中红色的上调基因309个,绿色的下调基因267个。
图1(a)和(b)显示了火山图和集群热图。
2.2 DNA 甲基化分析
在甲基化数据集GSE63695中,从OA和对照组之间筛选出15069个DMP,包括10296个高甲基化位点和4773个低甲基化位点,火山图如
图2(a)所示。CpG位点分布在TSS1500、TSS200、5'UTR、1stExon、Body和3'UTR等不同基因区域。如
图2(b)所示,在1stExon区域有151个CpG位点的DMGs,在3'UTR区域有505个CpG位点的DMGs,在5'UTR区域有695个CpG位点的DMGs,在Body区域有2779个CpG位点的DMGs,在TSS1500区域有795个CpG位点的DMGs,在TSS200区域有261个CpG位点的DMGs。高甲基化基因和低甲基化基因不同区域的CpG位点分布相似,约50%的甲基化改变发生在基因Body区域。将基因相关特征TSS1500、TSS200、5'UTR和第一外显子(1stExon)统称为启动子区。启动子区的甲基化状态可以影响基因的转录水平,从而影响基因的表达和功能。这种调控机制在生物体的发育、分化和疾病等过程中发挥着重要作用,所以甲基化与启动子在调控基因表达方面密切相关。在这里,我们仅讨论非基因Body区域的甲基化情况与基因表达的调控关系。
将高低表达的DEGs和高低甲基化的DMGs的其余5个区域分别取交集,发现57个基因是高甲基化低表达(Hyper-down),37个基因是低甲基化高表达(Hypo-up),32个基因是低甲基化低表达(Hypo-down),69个基因是高甲基化高表达(Hyper-up)。Hyper-down基因的CpG位点按基因区域分布为:1stExon(2个基因)、3'UTR(5个基因)、5'UTR(7个基因)、TSS1500(10个基因)、TSS200(2个基因),见
图3(a)。CpG位点分布在Hypo-up基因的不同区域,包括1stExon(1个基因)、3'UTR(4个基因)、5'UTR(4个基因)、TSS1500(8个基因),见
图3(b)。Hypo-down基因的CpG位点按基因区域分布为:1stExon(1个基因)、3'UTR(3个基因)、5'UTR(7个基因)、TSS1500(4个基因),见
图3(c)。CpG位点分布在Hyper-up基因的不同区域,包括1stExon(3个基因)、3'UTR(2个基因)、5'UTR(8个基因)、TSS1500(12个基因),见
图3(d)。由于有的基因在两个或两个以上的区域甲基化,所以上述结果存在基因名重复的问题,经过整理最终得到145个交集基因(Methylated and differentially expressed genes,mDEGs)。结果表明,异常甲基化表达基因的CpG位点主要分布在TSS1500区域,特别是在Hyper-up和Hyper-down基因中。Hypo-down和Hypo-up基因的CpG位点在TSS200区域没有且在1stExon区域数量最少。
2.3 相关性分析
为了探究基因表达水平和不同区域甲基化之间的潜在联系,用Spearman相关系数对基因表达数据与对应区域甲基化数据进行计算,根据相关系数|
r|>0.3,
P<0.05筛选出36个显著相关的重要基因,结果如
表1所示。
2.4 功能富集分析
利用clusterProfiler包对筛选出的145个重要基因进行GO功能分析和KEGG通路分析。结果显示,共富集到GO功能通路有242条,其中汇集到BP的有217条、CC有5条、MF有20条,可视化前10个富集项,结果如
图4(a)所示。主要富集到的BP有:细胞黏附的正向调节、白细胞活化的正向调节、细胞活化的正向调节、白细胞趋化性、脂肪细胞分化、上皮细胞增殖、细胞趋化作用、淋巴细胞活化的正向调节、白细胞游走、细胞数量的内稳态等;主要富集到的CC有:含胶原的细胞外基质、基底膜、阳离子通道络合物、膜筏、膜微域、电压门控钾通道复合物等;主要富集到的MF有:受体配体活性、硫化合物结合、蛋白酪氨酸/苏氨酸磷酸酶活性、MAP激酶酪氨酸磷酸酶活性、MAP激酶酪氨酸/丝氨酸/苏氨酸磷酸酶活性、核糖皮质激素受体结合、肝素结合、MAP激酶磷酸酶活性、糖胺聚糖绑定、核受体活性等。富集到的KEGG通路有41条,可视化前10条通路,结果如
图4(b)所示,主要包括甲状旁腺激素的合成、分泌及作用、AMPK信号通路、类风湿性关节炎、MAPK信号通路、肥厚型心肌病、PI3K-Akt信号通路、胰腺分泌、阿米巴病、蛋白质消化吸收、皮质醇的合成和分泌等。
2.5 PPI网络分析
为了解基因与所编码的各个蛋白质之间的相互关联功能,使用STRING数据库分析了上述145个重要基因的PPI网络,并用Cytoscape软件可视化,如
图5(a)所示,进一步选择插件cytoHubba选取前30个TOP基因作为关键基因,见5(b)。
2.6 hub基因的鉴定及验证
将基因表达数据和不同区域甲基化数据筛选出的36个重要基因与PPI网络筛选出的30个TOP基因联合分析,共得到5个hub基因,见
图6,它们分别是COL1A1、COL8A1、NR4A2、PTHLH、PTN。
2.7 关键基因的验证
为了确认以上5个基因的诊断价值,利用pROC包分析了训练集GSE55235、GSE55457的合并表达数据,绘制ROC曲线,如
图7所示。结果显示,在训练集中这5个基因的
AUC值均在0.7以上,见
图7(a);同样在甲基化数据集中,它们的
AUC值均在0.8以上,见
图7(b)。以上结果说明,这5个基因可以作为OA的关键基因,在OA中表现出较高的诊断价值。
2.8 ceRNA网络的构建
此外,从GEO数据库下载了2个数据集circRNA(GSE175959)和miRNA(GSE105027),并识别差异表达的circRNAs(Differentially expressed circRNAs,DECs)和差异表达的miRNAs(Differentially expressed miRNAs,DEMs)。将DEMs和DEGs靶向的miRNA取交集,获得了目标DEMs(Differentially expressed target miRNAs,DETMs),将DEGs与DEMs靶向基因交叉获得了靶DEGs(Differentially expressed target miRNAs,DETGs)。由24个circRNAs(3个下调,21个上调)、42个miRNA(11个上调,31个下调)和84个mRNA(30个下调,54个上调)构建了一个circRNA-miRNA-mRNA调控网络,如
图8所示。
图8中,深绿色的菱形代表低表达的circRNA,浅绿色的菱形代表高表达的circRNA;深红色的三角形代表低表达的miRNA,浅红色的三角形代表高表达的miRNA;深蓝色的圆形代表低表达的mRNA,浅蓝色的圆形代表高表达的mRNA。
在hub基因中,COL1A1的表达在circRNA高表达、miRNA低表达和低甲基化的调节下上调,COL1A1与3个miRNA相互作用,分别是hsa-miR-6740-5p、hsa-miR-6756-5p、hsa-miR-98-5p,与此同时,这3个miRNA也与特定的circRNA结合,它们是hsa_circ_0082680、hsa_circ_0000787与hsa_circ_0027914、hsa_circ_0088200。NR4A2与两种miRNA相互作用,分别是hsa-miR-5001-3p、hsa-miR-1224-3p。有趣的是,这两种miRNA同时与hsa_circ_0032363相互作用,该circRNA可能是调控的关键节点,有待进一步研究。
3 结论
骨关节炎(OA)是最常见的慢性骨关节病和退行性关节病之一。近年来的研究表明,基因表达动力学的表观遗传调控机制正逐渐成为重要的研究方向。本研究使用两个GEO基因表达谱数据集和一个甲基化谱数据集进行了生物信息学分析,以确定OA发病过程中涉及的生物学机制。共鉴定出145个常见DEGs,包括79个上调DEGs和66个下调DEGs。通过对这些基因进行GO和KEGG分析,发现差异基因最显著的富集功能和途径。
对得到的145个共同基因进一步分析,最终确定了5个hub基因,分别是COL1A1、COL8A1、NR4A2、PTHLH和PTN。对hub基因在训练集上绘制ROC曲线,发现其在骨关节炎中具有良好的诊断性能。
骨骼的主要有机成分是Ⅰ型胶原,它由1条Ⅰ型胶原α2链基因(COL1A2)编码的1条前α2(Ⅰ)链和2条胶原α1链基因(COL1A1)编码的2条前α1(Ⅰ)链组成。COL1A1是几种胶原蛋白之一,在RNA和蛋白质水平方面表现出丰富的差异性。胶原蛋白是软骨的主要结构成分,胶原失衡对OA有致病作用。Ⅰ型胶原蛋白是瘢痕组织中最丰富的胶原蛋白,是组织愈合的最终产物。有研究表明,COL1A1在晚期OA患者的滑膜以及OA诱导的小鼠和转化生长因子-
β(Transing growth factor-
β,TGF-
β)刺激的人成纤维细胞的滑膜中上调
[10]。既往研究表明,Ⅰ型胶原蛋白在正常软骨细胞中很少见,但在OA软骨中高表达,尤其是在OA晚期
[11]。
COL8A1是一种Ⅷ型胶原蛋白,主要在角膜、软骨和血管组织中表达。在正常的生理活动中,COL8A1与血管生成密切相关,并促进平滑肌细胞的增殖和迁移。COL8A1可能参与稳定软骨细胞的细胞外基质,沉默COL8A1会抑制COL2A1基因的表达,从而加速OA的发展进程。此外,COL8A1还与肿瘤发生密切相关,一些研究表明,COL8A1在肺癌、结肠癌和肝癌中的表达异常升高。研究发现,COL8A1可能通过调节上皮-间充质转化(EMT)和黏着斑相关通路促进人胃腺癌和结肠腺癌的发展
[12-13]。
NR4A2全称为细胞核受体第四亚科A组第二成员。NR4A2基因编码一个能结合DNA的锌指蛋白,在关节炎的发生和发展中起一定的作用。NR4A2(Nurr1)可能是TNF-α和核因子κB(NF-κB)信号通路下游的一个有希望的治疗靶点。该转录因子与NR4A1(Nur77)和NR4A3(NOR1)是NR4A受体家族的成员。在炎症反应中,NF-κB和环磷酸腺苷反应元件结合蛋白(CREB)直接与NR4A2启动子结合,并迅速诱导其在软骨细胞、滑膜细胞、内皮细胞和免疫细胞中的表达
[14-17]。NR4A2在类风湿关节炎和银屑病关节炎患者的炎症滑膜组织中以及OA患者的软骨中也高度表达
[18-21]。
甲状旁腺激素样激素(PTHLH)也称为PTHrP,在终末软骨细胞分化和软骨内骨形成中起抑制作用,miR-195下调可加速OA的发生和发展,而PTHrP是miR-195的新靶点
[22-23]。研究发现,PTHLH是miR-29b的靶标,PTHLH过表达会导致CDK4增强,从而促进CDK4诱导的RUNX2泛素化。与另一研究结果一致,PTHLH通过诱导CDK1介导的RUNX2泛素化和蛋白酶体降解来保护软骨细胞免于过早肥大
[24]。
促胰激素(PTN)是一种分泌的肝素结合肽,参与细胞生长、迁移和血管生成。除了直接促进血管生成外,PTN还通过促进内皮祖细胞的增殖和迁移来影响其他细胞的行为,并促进血管微环境的形成。血管生成在OA的发病机制中起关键作用,通过促进炎症细胞的迁移和局部疼痛受体的增加,导致组织损伤和疼痛。因此,PTN可以通过促进血管生成来加重骨关节炎的发展
[25]。然而,在骨关节炎的早期阶段,PTN在软骨细胞中重新表达,并可以在患者的滑液中原位检测到。PTN通过刺激细胞外基质合成,减少降解基质金属蛋白酶和诱导其抑制剂来增强软骨形成,PTN还可以减少促炎因子,例如一氧化氮和血管内皮生长因子
[26]。此外,PTN刺激软骨细胞聚集和增殖。因此,PTN似乎介导骨关节炎软骨的修复和保护过程,可能是治疗骨关节炎的一个有希望的因素。
通过circRNA-miRNA-mRNA调控网络,揭示了与DEGs结合的miRNA与DEGs来源的circRNA之间潜在的“海绵”关系,并随后构建了滑膜组织中circRNA-miRNA-mRNA调控网络,确定了一些可能为OA的诊断和治疗提供新靶点的miRNA和circRNA。
综上所述,本研究揭示了OA潜在的关键基因、通路和circRNA-miRNA-mRNA调控网络,结果有助于提高对OA潜在发病机制的认识。此外,这些hub基因和circRNA可能有助于OA的早期诊断和治疗。