生物信息联合免疫分析鉴定幼年皮肌炎相关基因

赵威 魏广友

1安徽理工大学医学院 安徽淮南 232000;2亳州市人民医院儿科

幼年皮肌炎(juvenile dermatomyositis,JDM)是一种罕见的慢性自身免疫性疾病,是特发性炎症性肌病之一。本病发病于儿童期,又称为儿童皮肌炎,发病率约为2.5/100万,患病率为2.5/10万［1］,平均发病年龄为81.97个月,以女童为主［2］。JDM特征是对称性近端肌肉无力,血清肌酶浓度升高,以及病理性皮疹［1］。JDM皮肤病表现出更大的活动性,比肌肉病更难治疗,并且儿童钙质沉着症和溃疡的风险较高［3］。JDM可导致严重的后遗症,包括营养不良性钙化、持续性皮肤炎症和功能受限的器官受累［4］。有研究表明,皮肌炎诊断后5年内发生肿瘤的风险很高,26%～70%的患儿在诊断皮肌炎后1年内发生肿瘤［5］。然而,JDM发病机制在很大程度上仍然是未解之谜,包括易感基因、环境应激源、免疫和非免疫介导机制在内的一系列因素,被认为可能与JDM的发病机制有关［6］。

近几年,随生物信息学发展,探索基因表达谱在疾病发病机制中的作用是一种有效的方法。JDM队列的基因表达谱和临床信息从NCBI基因表达综合数据库［7］下载(GEO;https://www.ncbi.nlm.nih.gov/geo/)。此外,为了确定与JDM相关的中枢基因,对患者和健康对照样本中的差异表达基因( differentially expressed gene,DEG)进行了比较。免疫系统参与疾病发展已为人所知多年,JDM由自身免疫结缔组织病变引起的,因此,本研究试图通过生物信息与免疫分析鉴定JDM的相关基因,把免疫系统分为免疫细胞和免疫功能两部分,并使DEG与其进行免疫评分及相关性分析,随后运用PPI网络选出hub基因,最后筛选出与免疫最为相关的基因。运用选出的基因做风险列线图、预测模型及验证,从而探索相关基因在诊断及治疗JDM的潜在价值［8］。

1 对象与方法

1.1研究对象 GEO是目前全球最大最全面的公开可用的基因表达数据库,保存着世界各国研究机构上传的高通量测序数据。在本研究中,通过替代变量分析方法将两个数据集(GSE11083和GSE11971)合并为一个集成数据集,来消除批次差异。GSE11083使用GPL570平台(基于:Affymetrix Human Genome U133 Plus 2.0 Array)生成,GSE11971使用GPL96平台(基于:Affymetrix Human Genome U133A Array)生成。另外一组数据集GSE100152使用GPL6884平台(基于:Illumina HumanWG-6 v3.0 expression beadchip)生成,用来进行预测模型的验证。GSE11083数据集包括28例健康对照儿童血液样本和27例JDM患者血液样本。GSE11971数据集包括4例健康对照儿童肌肉组织样本和19例JDM患者肌肉组织样本。GSE100152数据集包括9例健康对照儿童血液样本和40例JDM患者血液样本。所有数据集样本年龄均<18岁,且符合 1975年 Bohan / Peter提出的多发性肌炎、皮肌炎诊断标准［9］。

1.2研究方法

1.2.1差异基因的筛选 数据集GSE11083和GSE11971在R软件中进行预处理,包括背景校准、归一化和log2转换。如果预处理多个探针对应于共同基因,则使用平均值当做相应的表示值。此外,Bioconductor的替代变量分析(SVA)软件包用于消除批次效应。然后,运用Limma软件包选择筛选JDM和健康对照儿童之间的DEG,并设置DEGs的倍数变化| log2Fold change(FC)|>1和P值<0.05为筛选条件。最后,分别使用 R语言pheatmap软件包和gglot软件包绘制热图和火山图。

1.2.2ssGSEA计算免疫评分 单样本基因集富集分析(single sample gene set enrichment analysis, ssGSEA)。首先对上述所得DEG进行秩次标准化,然后利用经验累积分布函数计算富集分数(ES)。运用R语言GSVA包实现ssGSEA分析,即可得到每个样本基于ssGSEA算法的基因集。结合免疫细胞与免疫功能的数据分析,使用R语言pheatmap软件包绘制数据集与免疫细胞及免疫功能的热图。

1.2.3疾病组与正常组之间免疫评分的差异分析 首先对疾病组样本与对照组样本进行免疫细胞之间的差异性分析,然后再进行免疫功能之间的差异性分析。P<0.05表明差异有统计学意义,使用ggpubr软件包比较和可视化JDM和对照样本之间14个免疫细胞的水平,13个免疫功能的水平。

1.2.4PPI网络及筛选Hub基因 使用检索相互作用基因数据库(String)的搜索工具(版本11.5;www.String-db.org)［10］和最小相互作用分数构建DEG的PPI相互作用网络,≥ 0.4定义为截止值。通过Cytoscape可视化PPI网络,最后运用cytoHubba插件中的MCC算法来筛选10个Hub基因［11］。

1.2.5筛选免疫最相关基因 上述方法选出10个Hub基因,运用ggcorrplot软件包,使其与免疫细胞及免疫功能进行相关性热图分析,最后筛选出5个与免疫最为相关的基因。

1.2.6构建预测模型与模型验证 利用上述步骤选出的5个与免疫最相关基因构建预测模型,使用rms软件包建立疾病与5个基因的列线图来预测发病风险,并利用校准曲线评估列线图的可靠性。使用ROCR软件包分别绘制出整体和单独5个基因与JDM的ROC曲线图来预测诊断价值,并计算出ROC曲线下面积(Area under the ROC,AUC)及其95%可信区间。选取另外一组独立数据集GSE100152用于预测模型的验证。再次使用上述5个基因,通过ROCR软件包绘制整体和单独5个基因ROC曲线,并计算出AUC及其95%可信区间。AUC越大,说明预测性能越好。

2 结果

2.1DEGs的筛选鉴定 将GSE11083和GSE11971两个数据集结合起来,以消除批量效应,共获得29个基因差异表达,其中22个基因上调,7个基因下调。前5位上调基因为ISG15、IFIT5、IFI44L、OAS2和LGALS3BP,前5位下调基因是(EIF1AY)、DDX3Y, USP9Y, KDM5D, RPS4Y1。DEGs的热图和火山图见图1。

图1 DEGs的热图和火山图

2.2ssGSEA算法计算免疫评分 使用ssGSEA来评估免疫与基因之间的关联,计算出的29个DEG与14种免疫细胞以及13种免疫功能的相关性热图,见图2。

图2 DEG与免疫细胞及免疫功能分析热图

2.3免疫评分的差异分析 免疫细胞的差异分析箱型图显示,与对照组相比,JDM患者的激活树突状细胞(activate dendritic cells,aDCS)表达水平较高,见图3。免疫功能的差异分析箱型图显示,与对照组相比,JDM患者的树突状细胞的抗原提呈细胞(antigenpresentingcell,APC)、非经典炎症(Parainflammation)、I型干扰素反应(T1- IFN),这三种免疫功能表达水平较高,见图4。

图3 JDM与对照组的免疫细胞差异分析

图4 JDM与对照组的免疫功能差异分析

2.4PPI网络及筛选Hub基因 去除孤立节点后,构建了由25个节点和84条边组成的PPI网络,由Cytoscape可视化PPI网络,见图5。PPI网络中通过MCC算法得到的10个DEG被定义为hub基因,分别为IFI27、IFIT3、IFI44L、ISG15、IFIT1、MX1、MX2、OAS2、OAS1、GBP1,见图6。

图5 DEG的PPI网络

图6 PPI网络的Hub基因

2.5筛选免疫最相关基因 由上述选出的10个Hub基因,根据与免疫相关性,最终筛选出5个最为相关的基因,分别为GBP1、IFIT1、IFIT3、OAS2、MX2,见图7。

图7 10个Hub基因与免疫细胞及免疫功能的相关性

2.6构建预测模型与验证 队列中基因的列线图表明模型具有较高的发病率,见图8A;列线图校准曲线在该队列中显示出良好的一致性,见图8B。ROC曲线整体预测模型的AUC为0.825(95%CI:0.730～0.920)。各基因AUC分别为:GBP1为0.783(95%CI:0.677～0.889)、IFIT1为0.757(95%CI:0.648～0.866)、IFIT3为0.733(95%CI:0.621～0.845)、OAS2为0.790(95%CI:0.692～0.888)、MX2为0.773(95%CI:0.671～0.874),见图8C。通过另一组样本验证整体模型的AUC为0.889(95%CI:0.758～1.000),单各基因的AUC分别为:GBP1为0.622(95%CI:0.471～0.773)、IFIT1为0.467(95%CI:0.296～0.638)、IFIT3为0.486(95%CI:0.314～0.658)、OAS2为0.828(95%CI:0.659～0.996)、MX2为0.535(95%CI:0.347～0.722),见图8D。最终得出OAS2基因可能具有较高的预测诊断价值。

图8 预测模型与验证图

3 讨论

JDM由自身免疫结缔组织病变引起,通常包括自身抗体阳性和免疫异常等疾病,其特征是皮肤损伤、肌痛或及进行性肌肉无力,导致严重的皮疹,其并发症甚至可能危及生命。然而,JDM的确切发病机制尚未明确描述。在本研究中我们通过2个数据集进行整合后,得到JDM的差异基因。运用DEGs与免疫细胞及免疫功能做相关性分析,得出DCS、APC、T1-IFN及非典型炎症有较强的差异性。通过PPI网络选出10个hub基因,运用hub基因再次筛选出与免疫最相关的5个关键基因:GBP1、IFIT1、IFIT3、OAS2、MX2,所绘制列线图表明模型具有较高的发病率。且校准图显示具有良好的一致性。通过预测模型及验证分析,最后认为OAS2基因可能具有较高诊断JDM的潜在价值。

JDM是一种具有干扰素(IFN)特征的系统性自身免疫性疾病,发病机制及其干扰素特征的病因尚不清楚［12］。有研究表明该机制可能通过干扰素途径引起信号转导的上调［13］。Brian等［1］研究学者发现,肌肉和外周血中I型干扰素(α)诱导的DCS基因上调似乎是疾病发病的核心。抗-MDA5抗体阳性的JDM在血清和受损皮肤中表现出高T1-IFN特征。MDA5抗体阳性DM患者血清和血管中的高T1-IFN特征表明,T1-IFN可能在这些患者的血管病变中起重要作用［14］。Emily等［15］发现JDM的外周血单核细胞中I型干扰素信号增加;Roberson等［16］研究发现,临床上JDM患者的外周血单核细胞具有持续的免疫激活作用,并富含IL-1信号,JDM皮肤和肌肉都显示了I型干扰素激活的证据,以及与抗原提呈和细胞呼吸基因表达降低相关的基因。血液中的I型干扰素途径标志性生物标记物与JDM患者的皮肤皮肌炎病区和严重程度指数评分高度相关［17］。结合我们所得到的免疫功能分析显示,JDM样本与健康组样本中I型干扰素差异有统计学意义。综上表明T1-IFN可能是JDM的潜在特异性分子途径。

OAS2即寡腺苷酸合成酶2,其属于OAS基因家族。OAS2是一个核苷酸结合的寡聚结构域,包含蛋白质2(NOD2)结合伙伴,导致RNAse-L活性增强,提示OAS2和其他先天免疫信号通路之间可能存在必要相关性［18］。干扰素诱导的双链RNA激活酶是OAS蛋白。OAS基因家族的表达在JDM中受到高度调节,类似于对dsRNA病毒感染的免疫反应,有研究报道发现JDM肌肉活检增加了OAS家族基因的表达,特别是OAS1、OAS2、OAS3和OAS4［19］。较早有研究表明过度角化的细胞可能是皮肌炎患者皮损形成的原因,其中OAS基因可能激活一种凋亡细胞死亡机制［20］。同时,Ouyang等［21］研究发现,JDM的所有OAS基因都受到调节,并且99%的OAS基因家族网络显著上调。因此,筛选的OAS2基因可能为JDM的一个重要发病基因。

综上所述,本研究通过生物信息联合免疫分析鉴定了OAS2基因可能成为JDM诊断的生物标记物及潜在的治疗靶点,并提出T1-IFN可能是调节JDM的潜在分子途径。此次研究为JDM发病机制的探索提供新的思路,并且可能挖掘出相关基因在诊断及治疗JDM的潜在价值。然而,该研究存在一些局限性。首先,我们只关注了免疫分析相关的基因,可能忽略了其他与JDM相关的潜在基因。此外,本研究是基于生物信息分析与解读,还需体内外实验来进一步验证。