硬皮病(SSc)是一种可以引起皮肤和其他器官纤维化的自身免疫性疾病,是一种最具破坏性的风湿病,皮肤增厚或硬化是该病的一个显著特征
[1, 2]。研究表明,硬皮病的发病体现在基因遗传、环境因素和病毒感染的复杂且相互作用上,但硬皮病的发病机制尚未完全阐明
[3]。因此,探索硬皮病的潜在分子机制,发现更有效的诊断技术和更可靠的分子标志物对于改善硬皮病的诊断和治疗是至关重要的
[4]。
线粒体主要功能包括氧化磷酸化、能量转换及储存钙离子等。其中,线粒体最重要的功能是通过能量转化为细胞各项活动提供能力,在成纤维细胞的多个过程中起着至关重要的作用
[5, 6]。病理条件下,线粒体损伤促进固有免疫系统的异常激活,能导致自身炎症或自身免疫性疾病
[7]。研究表明,与硬皮病发病机制有关的几个因素,如缺氧或慢性炎症反应中的细胞应激,均可诱导线粒体损伤
[8]。而线粒体损伤可激活硬皮病成纤维细胞中Smad3信号增强的纤维化转录程序,且硬皮病中不受控制的TGF-β信号与线粒体损伤和成纤维细胞激活可能存在因果关系
[9]。
尽管线粒体和硬皮病的联系不断被证实,仍然缺乏系统分析线粒体相关基因与硬皮病关系的研究。如哪些线粒体相关基因在硬皮病中有关键作用等,仍有待进一步考证。因此,为探索线粒体相关基因对硬皮病的影响,我们基于基因表达综合数据库(GEO),使用3种机器学习算法鉴定出硬皮病关键基因,通过人工神经网络构建了硬皮病的诊断模型,以期该模型为探索硬皮病的发病机制提供一个新的视角。
1 材料和方法
1.1 GEO芯片数据收集及差异表达基因获取
从GEO数据库(
https://www.ncbi.nlm.nih.gov/geo/)下载硬皮病患者的基因表达数据。纳入GSE95065、GSE76807及GSE59785数据集。从MitoCarta3.0数据库获取2030个线粒体相关基因的数据。通过R(4.2.1版)软件“sva”包将GSE59785与GSE95065合并作为实验组,然后从实验组数据集中提取2030个线粒体相关基因的表达量。使用R中“limma”包对实验组进行差异分析,设置过滤条件为
|logFC
|>0.5,adj.
P<0.05,筛选得到上调和下调的差异表达基因(DEGs)。GSE76807数据集作为验证组(
图1,
表1)。
1.2 Metascape分析
使用在线工具Metascape(
http://metascape.org/gp/index.html)对DEGs进行分析,过滤条件为:
P<0.01;且每个功能或通路中至少有3个差异基因;富集因子>1.5,默认显示前20个富集最显著的功能或通路并绘制网络图和柱状图。
1.3 GO和KEGG富集分析
采用R软件“clusterProfiler”包进行GO和KEGG分析。设置参数P<0.05代表富集结果显著,用R“ggplot2”包将GO及KEGG富集分析结果可视化。
1.4 机器学习筛选硬皮病关键基因
用最小绝对收缩和选择算子(LASSO)、支持向量机(SVM)及随机森林算法筛选硬皮病关键基因。首先用R包“randomForest”进行DEGs筛选,基尼系数法计算重要性评分>0.5的基因被视为纳入标准,绘制随机森林树图形和基因重要性图形
[10]。再对DEGs使用LASSO回归和SVM算法筛选特征基因,用R中“glmnet”和“caret”包进行筛选。最后将3种机器学习算法所得结果基因用R软件中“venn”包进行韦恩分析,取交集后得到SSc的关键基因。
1.5 人工神经网络的构建与验证
根据特征基因的表达水平,对特征基因进行基因评分,若基因的表达水平高于基因表达中位值,基因得分记为1,否则记为0,之后输出所有基因得分的结果
[11]。使用R中“neuralnet”包和“NeuralNetTools”包构建特征基因人工神经网络(ANN),以疾病特征基因的表达量与基因权重作为疾病分类的依据。设置5个隐藏层作为模型参数,根据5个隐藏层节点和权重得到输出层从而构建ANN。用“Caret”包计算ANN模型的10折交叉验证,以优化模型去除过拟合,混淆矩阵函数计算结果的准确性,用“pROC”绘制受试者操作特征(ROC)曲线,计算曲线下面积(AUC)值,评估人工神经网络对疾病预测的准确性。再对实验集整体模型进行ROC验证。最后用外部数据集构建ROC曲线,进一步验证模型准确性。
1.6 RT-qPCR验证关键基因
动物实验设计:6~8周龄SPF级雌性小鼠6只[河南斯克贝斯生物科技股份有限公司,许可证号:SCXK(豫)2020-0005]。饲养于南阳理工学院动物实验中心,昼夜光照交替,自由饮食,温度22±2 ℃,湿度40%~60%,动物实验经南阳理工学院伦理委员会批准(批准号:南理工动伦审[2023]014)。6只小鼠随机分为模型组和对照组,3只/组。剃去中央区被毛后,对照组小鼠连续4周皮下注射0.9%氯化钠注射液0.1 mL/d。其余小鼠采用博来霉素,使用pH7.4的磷酸缓冲盐溶液(PBS)稀释,200 μg/mL,0.1 mL/d,连续皮下注射4周,建立SSc小鼠模型。取下背部皮肤以10%福尔马林溶液固定、常规脱水、石蜡包埋、制作切片进行HE染色及Masson胶原纤维面积对比,观察皮肤外观及组织病理学样变,确定造模成功。造模成功后,用trizol法提取总RNA,根据体系逆转录cDNA并扩增关键基因,采用2
-ΔΔct值检测关键基因mRNA相对表达量。trizol试剂盒、逆转录试剂盒、RT-qPCR试剂盒(赛默飞世尔)。引物由武汉合测基因科技有限公司设计合成,引物序列见
表2。
1.7 免疫细胞含量分析
采用CIBERSORT算法对浸润免疫细胞进行相关性分析。设置参数perm=1000,P<0.05为过滤条件,获取硬皮病组和正常组中22类免疫细胞表达的相对含量。分析2组免疫细胞的差异和各免疫细胞的相关性。使用R软件中的“corrplot”可视化。
1.8 已鉴定基因与免疫细胞浸润相关性分析
利用R中的Spearman等级相关分析,研究已鉴定的基因潜在生物标志物与免疫细胞浸润水平的相关性,用R中的“ggplot2“包可视化。
1.9 统计学方法
使用SPSS 21.0软件进行统计学分析,计量资料以均数±标准差表示,采用独立样本t检验分析SSc和对照组之间的基因表达差异,以P<0.05为差异有统计学意义。
2 结果
2.1 DEGs获取结果
共获取硬皮病DEGs 24个,其中上调基因11个,下调基因13个。绘制热图和火山图(
图2)。
2.2 Metascape分析
在Metascape分析中DEGs富集最显著的功能或通路为长时程抑制、正向调节蛋白磷酸化、AGE-RAGE信号通路在糖尿病并发症中的作用、炎症反应、CTLA4抑制信号传导、腺发育、细胞对肽的反应、酒精性肝病、线粒体组织的调节、ERBB2信号传导等(
图3)。
2.3 GO和KEGG通路富集分析结果
GO富集分析结果显示,DEGs的生物学功能主要与肽基-酪氨酸磷酸化、肽基-酪氨酸修饰、粘着斑、细胞-衬底结、蛋白质-高分子适配器活性、支架蛋白结合等有关。KEGG通路富集结果表明DEGs主要富集在长期抑制、EGFR酪氨酸激酶抑制剂耐药性、AGE-RAGE信号通路在糖尿病并发症中的作用、FOXO信号通路等(
图4)。
2.4 机器学习筛选硬皮病特征基因的验证
在随机森林树算法中,将参数“ntree”设置为500,计算交叉验证误差,选择30棵树作为交叉验证误差最小的最佳树数。用基尼系数法计算每个基因的重要性得分,得出重要性评分大于1的特征基因(
图5A、B)。通过Lasso回归和SVM算法,得到硬皮病的的特征基因(
图5C、D)。随后将3种机器学习算法得到的特征基因取交集,最终得到7个具有代表性的硬皮病关键基因(POLB、GSR、KRAS、NT5DC2、NOX4、IGF1、TGM2)(
图5E)。
2.5 ANN的构建与ROC验证结果
将3种机器学习算法筛选出的7个关键基因通过ANN分析优化各基因权重。ANN模型包括7个输入层、5个隐藏层和2个输出层。以基因权重之和与显著基因表达水平乘积为分类依据,ANN每条连线上的数字代表节点之间的权重。模型在建模组中识别正常和硬皮病样本的准确率分别为0.882和0.949。10折交叉验证的每个结果都用ROC曲线表示,平均AUC值超过0.980,证明了该模型的可靠性(
图6)。最后,建立一个ANN模型(
图7A),根据上述信息对硬皮病样本和对照样本之间的基因表达数据进行分类,该模型的整体AUC为0.984(
图7B)。在验证组数据集中观察7个模型基因的表达差异,AUC值为0.740(
图7C)。
2.6 RT-qPCR验证结果
使用RT-qPCR对硬皮病预测基因的mRNA表达情况进行检测,经过实验验证发现,与对照组相比,硬皮病中POLB(
P=0.004)、GSR(
P=0.029)、KRAS(
P=0.007)、NOX4(
P=0.019)、IGF1(
P=0.008)、TGM2(
P<0.0001)表达量明显上调,而NT5DC2(
P=0.001)表达量在硬皮病组中明显下调(
图8)。
2.7 免疫细胞浸润结果
使用CIBERSORT筛选出22种免疫细胞表达的相对含量。各免疫细胞正相关较高的为记忆静息CD4
+T细胞与滤泡辅助T细胞、树突状细胞活化与中性粒细胞;负相关较高的为活化NK细胞与静息肥大细胞及巨噬细胞M0(
图9)。
2.8 基因潜在生物标志物与浸润性免疫细胞的相关性分析
相关性分析结果显示,GSR与滤泡辅助T细胞呈负相关;IGF1与幼稚B细胞呈正相关、与静息树突状细胞呈负相关;KRAS与记忆激活CD4
+T细胞呈正相关、巨噬细胞M0呈负相关;NOX4与T细胞调节(Tregs)呈正相关、与单核细胞,肥大细胞激活呈负相关;NT5DC2 与CD4
+T细胞记忆静息呈负相关(
图10)。
3 讨论
在21世纪,硬皮病的研究和发展较前已经有了较大进步。但是在硬皮病线粒体相关基因方面的研究尚少且不成系统。有研究指出线粒体损伤和机体自身免疫相关, 但免疫介导的硬皮病与线粒体功能异常之间的关系尚未明确
[12]。因此,建立的基于线粒体相关基因的硬皮病诊断模型对后续临床诊疗的参考发展有十分重要的意义。
本研究比较了正常和硬皮病样本之间线粒体相关基因的表达差异。共筛选了7个线粒体相关基因(POLB、GSR、KRAS、NT5DC2、NOX4、IGF1、TGM2),构建了用于诊断的ANN模型。有研究表明,在小鼠脑和肾的高纯度线粒体中显示出POLB免疫反应性。另一方面,来自HEK293T细胞的线粒体制剂中POLB具有免疫反应性,当POLB过表达时,信号增加
[13]。有研究已经证明,POLB既存在于线粒体中,也是细胞器中大部分碱基切除修复活动所必需的
[14]。但尚未有系统研究表明POLB在硬皮病发病中的分子机制。GSR基因在氧化还原动态平衡和细胞对氧化应激的防御中起关键作用
[15]。GSR在人宫颈癌组织中的表达增加,GSR的RNA干扰可诱导宫颈癌细胞死亡
[16]。KRAS基因位于染色体12p12.1上,由6个外显子组成,KRAS的种系突变已在心脏-面部-皮肤(CFC)综合征、Noonan综合征和Costello综合征患者中被发现
[17, 18]。研究表明
[19],皮脂腺痣和线状皮脂腺痣综合征也可由受精卵后HRAS和KRAS突变引起。在皮脂腺痣的病例中,细胞增殖增加和肿瘤转化发生,可能是由于HRAS和KRAS基因变异导致Raf-MEK-ERK和磷酸肌醇-3激酶信号通路的异常激活
[20]。NT5DC2是一种50-核苷酸酶,可以催化核苷酸的水解,与精神疾病和癌症有关
[21]。近年来,NT5DC2被确定为一种癌基因,通过稳定表皮生长因子受体促进肝癌细胞增殖
[22]。曾有研究在肺纤维化的研究中检测到NT5DC2表达量明显改变
[23]。参与病理性纤维化过程发展的多种生长因子,最突出的是TGF-β,已被证明可以调节NOX的表达,特别是NOX4的表达
[24]。曾有研究证实了硬皮病真皮成纤维细胞中NOX4表达和ROS产生的增加,并表明NOX4刺激ROS产生介导了强有力的促纤维化作用,在这些细胞中观察到由于NOX4 siRNA抑制导致I型胶原生成的显著减少
[25]。TGM2参与许多病理反应,包括炎症、炎症因子趋化、肿瘤进展、伤口愈合、组织纤维化和免疫疾病
[26]。TGM2转酶活性通过与细胞外基质蛋白交联,促进胶原和细胞外基质的成熟和稳定,并参与疤痕组织、肝脏和心脏等多种组织的伤口愈合和纤维化
[27]。有研究发现,IGF1、RANTES和VEGF可视为早期或轻度硬皮病相关的血管生成细胞因子
[28]。研究发现
[29]与系统性红斑狼疮患者或健康对照组相比,硬皮病患者尤其是更严重的间质性肺疾病患者的血清IGF1水平显著升高。
免疫细胞浸润相关性分析结果显示特征基因与滤泡辅助T细胞、幼稚B细胞、静息树突状细胞、记忆激活CD4
+T细胞、巨噬细胞M0、T细胞调节、单核细胞、记忆静息CD4
+T细胞和肥大细胞激活等相关。免疫系统,尤其是细胞介导的免疫应答在硬皮病发病中起主导作用,而T细胞,特别是CD4
+T细胞的过度活化是硬皮病血管疾病、纤维化、体液免疫发生发展的关键因素
[30]。滤泡辅助T细胞是CD4
+T细胞的一个亚群,其特征是高水平的趋化因子受体CXCR5表达
[31]。结合目前已知的硬皮病发病机制和滤泡辅助T细胞的研究,示该细胞也可能通过包括IFNγ在内的机制参与硬皮病的发病
[30]。巨噬细胞可通过分泌促纤维化介质转化生长因子-b、CCL18、血小板衍生生长因子和成纤维细胞生长因子家族介质参与纤维化,诱导成纤维细胞增殖或胶原生成
[32]。肥大细胞可通过淋巴器官内的同源相互作用微调T细胞和B细胞活性,从而调节局部及全身炎症反应
[33]。有研究发现硬皮病患者的血液单核细胞和MDM不仅具有促纤维化特性,而且与促炎性巨噬细胞有一些共同的特征
[34]。而CD4
+T和CD8
+T细胞在未经治疗的早期弥漫性系统性硬化症患者的皮肤活检中积聚,并且数量超过浸润性B细胞
[35]。
综上所述,本研究发现硬皮病发病可能是通过调节线粒体相关基因POLB、GSR、KRAS、NT5DC2、NOX4、IGF1、TGM2实现的。我们还发现硬皮病可能通过调节滤泡辅助T细胞、幼稚B细胞、静息树突状细胞、记忆激活CD4+T细胞、巨噬细胞M0、T细胞调节、单核细胞、CD4+T细胞记忆静息和肥大细胞激活等功能来影响疾病的发展。本研究基于机器学习算法构建了ANN,为临床研究和基础实验提供了理论依据和数据支撑,并为探索硬皮病的发病机制提供一个新的视角。