通过生物信息学方法构建系统性红斑狼疮患者的基因表达亚群

马江磊，陈华秋，王光明，2*

（1.大理大学临床医学院，云南 大理 671000；2.大理大学第一附属医院基因检测中心，云南 大理 671000）

系统性红斑狼疮（systemic lupus erythematosus，SLE）是一种以免疫调节异常、自身免疫复合物沉积为特征的慢性自身免疫性疾病，常累及多脏器损害，甚至危及患者的生命。女性患病率高于男性，比例约为10∶1〔1〕。在过去的半个多世纪，随着人们对SLE认识逐渐增多，SLE患者的生存率和治疗效果有所提高〔2〕，但SLE临床症状表现多样，仍有部分患者因诊断延误或误诊，未得到及时治疗，从而产生严重后果。因此，建立准确、科学的分类和有针对性的治疗，对改善SLE患者的治疗效果有十分重要的意义。

SLE发病机制复杂，目前已知除了与紫外线照射、气候和地理分布等外界因素有关外，还发现多种与SLE有关的易感基因〔3〕。Wang等〔4〕通过对17 389名SLE患者和29 273名健康人群数据进行荟萃分析后发现，在欧洲及亚洲的总体人群中，信号转导和转录激活因子（signal transducer and activator of transcription，STAT）4基因rs7574865、rs10168266多态性与SLE发病风险密切相关。Sandling等〔5〕采用以信号通路为中心的方法，对瑞典958名SLE患者和1 026名健康个体的1 832个候选基因进行了免疫通路基因测序，确定了T淋巴细胞分化和先天免疫2个涉及SLE易感性的主要独立途径，通路风险评分有可能将SLE患者临床表现以外的分子进行分类，有助于临床随访和治疗的选择。

目前，针对SLE的基因治疗还不完善，为促进SLE治疗新策略的开发，本研究对1 254例SLE患者的转录组数据进行汇总分析，建立共识聚类表达分析，分成相应的亚群，并比较各亚群之间在临床特点、分子构成及细胞学功能方面的差异，为SLE的基因分型提供相关依据。

1 材料与方法

1.1 数据集的富集分析和主成分分析（principal componentanalysis，PCA）在NCBI的GeneExpression Omnibus（GEO）数据库（https：//www.ncbi.nlm.nih.gov/geo/）中搜索“systemic lupus erythematosus”，下载GSE121239、GSE65391和GSE154851的微阵列数据集，各芯片数据集中只包含SLE患者及健康对照者的组织。GSE121239由GPL13158平台文件提供，GSE65391和GSE154851分别来自GPL10558、GPL16699平台文件。使用R 4.1.3软件中的“limma”和“sva”数据包消除批次效应，进行归一化〔6-7〕。将数据合并后得到1 254例SLE患者（以下称“转录组”）和124例健康对照者（以下称“对照组”）的基因表达数据集。其中临床数据包括年龄、性别、SLE活动指数（systemic lupus erythematosus disease activity index，SLEDAI）等。

1.2 亚群的构建和亚群间临床特点的比较1 254例SLE患者的基因表达矩阵利用R软件中的“ggplot2”“rstatix”“reshape2”和“Consensus Cluster Plus”等数据包对矫正后的样本数据集进行分析，得到具有相似差异表达基因（differentially expressed genes，DEGs）的SLE患者分型。每个亚群的聚类得分越高，说明各亚群基因在基因表达上的相似性越高，然后把每个亚群视为一个聚类进行分析。

使用R软件中“ggplot2”和“ggpurb”数据包对各亚群的临床特点进行分析。以患者性别作为离散变量，分析男性和女性在各亚群中所占的比例；以年龄和SLEDAI作为连续变量，进行差异富集分析得到相关箱式图。

1.3 筛选亚群特异性基因及构建蛋白质相互作用（protein-protein interaction，PPI）网络将各亚群与对照组以及各亚群之间进行比较后，从中筛选出均值差异>0.2、调整后P<0.05的基因，得到各亚群对应的DEGs，将所得到的DEGs取交集，进一步筛选出只在某个亚群中表达上调的基因，即亚群的特异性基因。取各群前100个特异性基因上传到STRING（https：//stringdb.org/）网站，设置置信系数为0.4，生成特异性基因表达蛋白的PPI网络图，再将PPI网络图传输到Cytoscape软件中，以Degree算法计算得到节点数最多的前10个蛋白并对其进行分析。

1.4 各亚群基因集的富集分析对基因集富集分析（gene set enrichment analysis，GSEA）。GSEA前对基因列表文件和基因集文件进行对应转换，将每个组得到的列表和数据集文件传输到GSEA软件（version 4.1.0）中，设置最大运行值为5 000，得到各亚群中特异性基因在对照组中的表达分析结果。

1.5 GO功能富集分析和KEGG信号通路分析选择矫正后的基因数据集、各亚群中特异性上调基因文件以及含有连续变量的临床特征文件，使用R软件的“WGCNA”数据包进行分析。对样本临床特征与DEGs和模块之间的相关性分析后通过动态树切割法〔8〕，将基因数据集分成颜色不同的基因模块。应用R软件对分成的基因颜色模块进行GO功能富集分析和KEGG信号通路分析。GO功能富集分析主要对生物学过程、细胞成分和分子功能进行富集分析。KEGG信号通路分析主要对分子生物学途径进行分析。对KEGG信号通路分析的结果，筛选出每个模块中差异最显著的信号通路，分析各亚群与生物学信息之间更密切的联系。

1.6 统计分析用SPSS 26.0软件进行数据分析，数据以（±s）表示。组间比较采用单因素方差分析，不同组间的比较使用独立样本t检验，P＜0.05为差异有统计学意义。

2 结果

2.1 转录组数据的预处理及分组对“GSE121239”“GSE65391”和“GSE154851”3个数据集进行预处理及批次矫正后，得到总数据集的PCA图。3个数据集在批处理前，各组间无明显相关性。见图1A。消除批次效应后，总体基因分布均匀集中。见图1B。根据每个样本中的基因与SLE相关情况进行聚类分析，设置2~10个聚类，根据一致性评分结果，可将样本分为3个亚群，分别为亚群Ⅰ、亚群Ⅱ及亚群Ⅲ。表达差异富集模式图的分析结果也表明分为3个亚群时，组间表达模式显著不同。见图2。

2.2 亚群间临床特点的比较结果本研究纳入“GSE65391”和“GSE154851”数据集中（“GSE121239”数据集中未标注性别信息）患者1 066例，男124例，女942例，男女性别比为1.00∶7.60。把性别作为每个亚群的离散变量进行亚群间分析比较后发现，性别在各亚群间无显著差异，说明性别不是区分各亚群的要素。在3个亚群中，亚群Ⅰ的平均年龄为（15.24±5.10）岁，亚群Ⅱ的平均年龄为（16.09±5.76）岁，亚群Ⅲ的平均年龄为（15.09±6.00）岁，富集比较结果显示，亚群Ⅰ和亚群Ⅱ、亚群Ⅱ和亚群Ⅲ之间年龄差异有统计学意义（P<0.05）。在3个亚群中，亚群Ⅰ的SLEDAI为（4.03±3.11），亚群Ⅱ的SLEDAI为（1.88±2.80），亚群Ⅲ的SLEDAI为（2.68±3.22）。SLEDAI在各亚群中差异均有统计学意义（P<0.05）。见图3。

2.3 各亚群特异性基因与PPI网络图分析3个亚群的特异性上调基因数量分别为：亚群Ⅰ（2 134个）、亚群Ⅱ（105个）、亚群Ⅲ（2 289个），进一步筛选出各亚群中前10个表达差异最显著的DEGs。见表1。选择各亚群中上调前100个的特异性基因，对其表达的蛋白质进行PPI网络分析，使用Degree算法得到节点>10的蛋白枢纽图，其中，排名前10的蛋 白 分 别 是STAT3、TLR4、BRIX1、TLR2、PTEN、HSPA8、NHP2、NOP58、CAT、DDX18。见图4。

表1 各亚群中前10个上调的差异表达基因

2.4 各亚群的GSEA结果通过对每个亚群中的DEGs与对照组进行比较分析，发现各组中黑色线及灰色线均靠左聚集，GSEA结果显示各亚群的P值和FDR值均小于0.01，说明各亚群中的特异差异基因在对照组也存在明显的表达差异。见图5。

2.5 GO功能富集和KEGG信号通路分析结果在测试1~20阈值功率的结果中得到最佳功率值为9，在功率值为9的情况之下，基因之间的连通性符合无标度网络分布，可分为8个颜色模块。各基因颜色模块的分析结果见表2～5。

表2 生物学过程富集分析结果

GO功能富集分析结果显示，在生物学过程中，粉红色模块中的基因与血液凝固、止血途径以及血小板聚集等联系密切。在细胞成分中，黏着斑、细胞-基质连接以及血小板α颗粒等与粉红色模块密切相关。分子功能结果显示，粉红色模块主要与结合肌动蛋白、结合整合素、结合胶原等有关。从KEGG信号通路富集结果中，筛选出每个颜色模块中表达差异最显著的信号通路，分别为：自然杀伤细胞介导的细胞毒作用、核糖体、线粒体自噬、细胞凋亡、血小板活化、造血细胞谱系、破骨细胞分化、甲型流感，以进一步分析与亚群之间的关系。

将上述信号通路与对照组及3个亚群中的富集表达情况进行分析，结果显示，自然杀伤细胞介导的细胞毒作用、核糖体和造血细胞谱系在亚群Ⅲ中高表达，而在亚群Ⅰ中低表达；亚群Ⅱ在各信号路径中表达无明显差别，说明从生物学功能层面可以区别不同的亚群。

表3 细胞成分富集分析结果

表4 分子功能富集分析结果

续表4

表5 KEGG信号通路富集分析结果

3 讨论

本研究通过对SLE患者的转录组基因进行分析，根据共识聚类表达，将所有基因分为3个亚群，对各亚群的临床特点进行对比分析，探讨了年龄、性别以及SLEDAI与各亚群间的联系。利用GO功能和KEGG信号通路富集，对各亚群的生物学过程、细胞成分、分子功能和相关信号通路的关联情况进行了分析，为SLE基因亚型的认识提供相应的依据。在所纳入的1 254例SLE患者中，性别在各亚群中差异无统计学意义，而年龄和SLEDAI在各亚群中差异有统计学意义，说明年龄和SLEDAI在SLE的病理生理过程及预后中起到重要作用。有研究〔9〕发现，年龄越小，其临床表现越严重，且发生治疗副作用的风险更高，病死率也会增加。SLEDAI是评估SLE临床状态和治疗效果的关键指标，多年来，在临床工作中被广泛应用〔10〕。SLEDAI-2K是对SLEDAI的修改，当SLEDAI-2K>4时，提示存在狼疮疾病活动〔11〕。

利用STRING网站分析，得到了网络节点最多的蛋白，通过表达蛋白研究其对应的基因可以看出，排名靠前的分别是STAT3、TLR4、BRIX1、TLR2。其中，STAT3是节点最多的基因，但它并不在每个亚群的前10个特异性上调基因中。STAT3是转录信号传导子与激活子通路的重要成员，在调节炎症和免疫应答反应中起到关键作用〔12〕。当T细胞功能缺陷时，导致STAT3信号通路中白细胞介素（interleukin，IL）-17和IL-21的异常表达，IL-17细胞因子可以介导中性粒细胞和其他促炎细胞因子的释放，从而引起炎症反应的发生〔13-14〕。Chen等〔15〕发现在SLE患者中，IL-17表达上调且与STAT3强度呈正相关，说明STAT3和SLE的发生有着密切联系，这与枢纽图分析的结果相同。

先天免疫系统是宿主抵御病原体侵袭的第一道防线，Toll样受体（Toll-like receptor，TLR）是固有免疫受体的重要组成部分，其在B细胞参与自身免疫性疾病的机制中起到产生抗体、T细胞的抗原呈递和促进细胞因子生成等重要作用〔16-17〕。目前已知存在于人类的TLR有11种（TLR1~11），TLR4是一种跨膜受体，属于富含亮氨酸重复序列的超家族蛋白，在哺乳动物中是第一个被定义的TLR成员，它既能在质膜上发挥作用，也能在核内体上发挥作用〔18-19〕。TLR2是另一种TLR受体，当TLR2被识别并与配体结合后，会激活细胞内信号转导通路，引起免疫应答，可介导多种生物学过程，包括多发性硬化症、天疱疮、特发性皮炎、白癜风和1型糖尿病等〔16，20〕。Ma等〔21〕发现TLR4在狼疮肾炎的发生发展过程中具有重要的致病作用，可能为狼疮的治疗提供新的治疗策略。然而，也有研究发现，TLR7和TLR9是SLE患者中B细胞的效应功能的主要作用因子，而并不是TLR4和TLR2〔16-17，22〕。综上所述，对于TLR4和TLR2在SLE的发病发展中的作用还需要进一步的证实。

在过去的十几年中，随着精准医疗的迅速发展，越来越多治疗手段被应用到疑难疾病的诊疗当中，特别是与基因相关的研究中〔23〕。2021年，我国科研工作者对SLE患者的转录组数据进行分析，建立了一个关于SLE患者动脉粥样硬化风险预测模型〔24〕。在本研究中，通过利用SLE患者在转录组数据聚集高表达分析得到的3个基因亚群，并对各亚群及对照组进行差异基因表达、信号通路以及临床特点等分析比较，阐述了SLE可能的亚型，为SLE的精准治疗提供依据。当然，本研究仍存在几点不足之处，从GEO数据库中获取的数据集样本缺乏更加详细的临床信息；研究中的样本年龄偏小，还需要结合更多数量的成年SLE样本一起分析；基因的亚型是多方面的，需要更多的样本量以及其他组学的相互结合分析作为依据。