基于CIBERSORT反卷积算法研究强直性脊柱炎差异基因的免疫浸润机制及相关中药预测的生物信息学分析*

丁 灿，沈 浮，刘 鑫，杨 雷

（1.湖南中医药大学第一附属医院，湖南 长沙 410007；2.永州市中医医院，湖南 永州 425000）

强直性脊柱炎（ankylosing spondylitis，AS）是一种常见的血清阴性型脊柱炎，以慢性炎症和新骨形成为主要特征，常累及脊柱及骶髂关节[1]。流行病学提示其全球发病率为0.1%～1.4%，男性发病率高于女性，晚期AS可导致活动受限和关节僵硬，甚者致残而严重影响患者生活质量[2-3]。AS的确切原因尚不清楚，目前一般认为其与遗传及环境因素导致的自身免疫机制紊乱所致的慢性炎症相关[4-5]。近年来，越来越多学者认为，免疫细胞浸润在AS的发生、发展和治疗中起着至关重要的作用。如树突状细胞、巨噬细胞、自然杀伤细胞、辅助性T1细胞和辅助性T17细胞都与AS的发展有重要关系[6]。因此，评估免疫细胞浸润和确定免疫细胞成分的差异对阐明AS进展的分子机制和开发新的免疫治疗靶点具有重要价值。在临床实践中，AS的治疗主要依赖于非甾体抗炎药及生物制剂[7]，而此类药物长期使用具有心血管、胃肠道及肾脏不良反应的潜在风险[8]，而中药的较好疗效、低副作用和低成本，以及对于免疫失衡的调节干预，使得其成为治疗AS较为优越的替代疗法。

随着高通量遗传分析的出现，基因表达谱分析成为了鉴定各种疾病差异表达基因的有效方法，有助于探索发病机制和开发生物标志物。迄今为止，已有多项研究使用微阵列表达来鉴定AS发病机制中的关键基因[9-10]。CIBERSORT是一种分析工具，它使用微阵列数据或RNA测序数据来评估样本中免疫细胞的表达，并获得各种免疫细胞比率[11]。本研究应用基因芯片数据库（gene expression omnibus，GEO）中AS患者与正常人群的全血差异枢纽基因进行相关生物信息学分析，采用广泛用于评估鉴定22种免疫细胞相对含量的CIBERSORT反卷积算法进行差异基因的免疫浸润分析研究，并筛选与免疫浸润相关的关键生物标志物，从而预测免疫相关靶点及中药复方，旨在为AS的诊断和治疗提供依据。

1 资料与方法

1.1 数据收集 以“Ankylosing Spondylitis”为关键词检索GEO基因表达综合数据库，筛选得到“GSE25101”芯片为本次研究分析的对象。该芯片主要包含16例AS患者及16名年龄性别对应的健康人全血基因表达谱，通过GPL6947 Illumina HumanHT-12 V3.0表达珠芯片平台将探针转换成基因名称，利用R语言程序中“Limmar”安装包对每个样本的微阵列数据进行预处理，剔除缺失值，当基因对应于多个探针时，取最大值[12]。

1.2 关键差异枢纽基因的蛋白网络互作（PPI）及功能注释 在STRING数据库（https://string-db.org/）获取蛋白互作网络模型（protein-protein Interaction，PPI）网络，在Cytoscape软件中利用CentiScape插件计算每个靶基因的DC（度）值，以平均度值的2倍进行核心靶基因的筛选并进行可视化[13]。运用R语言程序中的“clusterProfiler”插件包进行差异基因的基因本体论（Gene Ontology，GO）与京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。

1.3 免疫浸润矩阵整理 本研究利用CIBERSORT反卷积算法进行免疫细胞浸润分析，包括T细胞、B细胞、NK细胞等22种免疫细胞在内的转录特征矩阵的模拟计算，设定次数为100次，对P＜0.05的数据进行后续分析。

1.4 免疫细胞间相关性分析及组间差异分析 为了分析各类免疫细胞之间的相关性，本研究在CIBERSORT反卷积算法以P＜0.05筛选可信样本数据，然后计算可信样本数据中不同免疫细胞间的Pearson相关系数，并运用秩和检验分析AS患者与健康人群之间的差异。

1.5 关键差异枢纽基因与免疫细胞浸润水平关系的相关性分析 为进一步研究关键差异枢纽基因与疾病相关性程度，本研究将筛选所得的差异枢纽基因与差异性较为显著的免疫细胞进行相关性分析，分析靶基因与各免疫细胞表达是否存在相关性并计算相关性系数r值。|r|＜0.4为弱相关，0.4＜|r|＜0.7为中等强度相关，|r|＞0.7为高度相关。

1.6 AS免疫浸润过程中潜在中药及复方的筛选预测 将AS差异核心靶基因及GO富集的与免疫浸润相关生物学过程（biological process，BP）导入Coremine Medical数据库（https://coremine.com/）及本草组鉴数据库（http://herb.ac.cn/）中，筛选具有潜在效应机制的中药。Coremine Medical数据库中的药物筛选标准为P＜0.05；本草组鉴数据库中筛选标准为P-value＜0.05，且FDR_BH＜0.05。通过阅读文献及古籍，基于上述所预测的药物进行复方预测，使得复方尽可能的包含较多的预测药物。

2 结果

2.1 AS差异表达基因获取结果 GEO数据库检索及筛选结果显示，AS组与正常组的差异基因为135个，其中上调基因为99个，下调基因为36个。图1A热图展示了前20个上调及下调基因的情况，图1B为差异基因的分布情况。

图1 AS 相关基因探针差异表达的热图及火山图

2.2 AS关键枢纽基因获取 利用STRING数据库对差异基因进行蛋白网络互作分析（PPI），将分析结果的TSV文件导入Cytoscape中进行可视化展示。图2A中节点颜色越深代表该节点的度值越大，在PPI模型中的重要程度越高。图2B可视化展示了度值大于2倍平均度值（15.8）的关键枢纽基因，其中度值大于15.8的有14个，分别是SNRPG、LSM3、PSMA6、MRPL22、RPL26、COX7C、PSMA4、RPS17、RPL31、SLIRP、NDUFS5、RPS24、NDUFB3、SRP14。

图2 AS 差异基因的PPI 网络图及核心靶基因

2.3 GO及KEGG富集分析结果 运用R语言程序中的“ClusterProfiler”插件包分析AS差异基因进行GO及KEGG富集分析。差异基因的GO分析中共获取生物过程（biological process，BP）77条，细胞组成（cellular component，CC）28条，分子功能（molecular fuction，MF）17条，对其中基因富集数目最多的10条进行条形图展示。根据基因的富集数目，表1展示了与免疫相关的生物学过程，主要涉及中性粒细胞脱粒（Neutrophil degranulation）、免疫反应中的中性粒细胞激活（Neutrophil activation involved in immune response）、先天免疫反应的调节（Regulation of innate immune response）、核糖体生物发生（Ribosome biogenesis）、氧化磷酸化（Oxidative phosphorylation）、NF-κB转录因子活性的正调控（Positive regulation of NF-κB transcription factor activity）、中性粒细胞趋化性（Neutrophil chemotaxis）、细胞防御反应（Cellular defense response）、白细胞介素-1介导的信号通路（Interleukin-1-mediated signaling pathway）、T细胞介导细胞毒性（T cell mediated cytotoxicity）。（见图3）在KEGG富集分析中，差异基因主要涉及核糖体合成、氧化磷酸化等相关通路，与代谢、炎性反应相关通路联系密切。

表1 AS 与对照组差异基因免疫浸润相关GO 富集分析列表

图3 AS 差异基因的KEGG 及相关基因可视化

2.4 AS患者与正常组的免疫细胞分布特征 通过CIBERSORT反卷积算法以P＜0.05样本进行筛选，最终获取32个可信样本，如图4A所示。热图左侧16个样本（GSM616684～GSM616699）为正常组，右侧16个样本（GSM616668-GSM616683）为AS组，结果显示，单核细胞（Monocytes）、CD8 T细胞（T cells CD8）、中性粒细胞（Neutrophils）、初始CD4+T细胞（T cells CD4 naive）在两组中含量均较高，且AS组中单核细胞（Monocytes）含量高于正常组。图4B显示，免疫细胞比例箱图中进一步展示了22种免疫细胞的浸润比例。

图4 AS 组与正常组中免疫相关细胞的浸润组成

2.5 免疫细胞间相关性分析 未活化的CD4+记忆性T细胞（T cells CD4 memory resting）与γδT细胞（T cell gamma delta）呈中度正相关（r=0.52），CD8+T细胞（T cells CD8）与活化的自然杀伤细胞（NK cell activated）呈中度正相关（r=0.51），而CD8+T细胞（T cells CD8）与中性粒细胞（Neutrophils）呈中度负相关（r=-0.65），初始CD4+T细胞（T cells CD4 naive）与未活化的CD4+记忆性T细胞（T cells CD4 memory resting）呈中度负相关（r=-0.69）。（见图5）

图5 AS 组中免疫细胞间的相关性分析

2.6 AS组与正常组的免疫细胞浸润差异分析 通过小提琴图对不同组织样本免疫细胞浸润差异分析进行可视化，以P＜0.05为显著差异性。结果显示，AS组血清中单核细胞（Monocytes）明显多于正常组（P=0.002），而正常组血清中调节性T细胞（T cells regulatoy，Tregs）明显少于AS组（P=0.03）。（见图6）

图6 AS 组与正常组间的免疫细胞浸润差异分析小提琴图

2.7 关键差异枢纽基因与免疫细胞浸润水平关系的相关性分析 将AS关键枢纽基因与上述所得差异性显著的免疫细胞（单核细胞、调节性T细胞）进行相关性分析，对相关性较为显著的基因进行展示。结果显示，SNRPG（r=-0.26）、LSM3（r=-0.14）、PSMA6（r=-0.27）、MRPL22（r=-0.18）与单核细胞呈负相关（P＜0.05），但相关性较弱；SNRPG（r=-0.30）、LSM3（r=-0.16）、PSMA6（r=-0.20）、MRPL22（r=-0.16）与调节性T细胞呈负相关（P＜0.05），但相关性较弱。（见图7）

图7 AS 关键枢纽基因与免疫细胞浸润水平相关性分析

2.8 干预AS免疫细胞浸润相关生物途径的中药及复方预测通过Coremine Medical及本草组鉴预测具有潜在治疗作用的中药。图8显示了AS核心差异靶基因及相关免疫反应生物学过程相关的中药，其中雷公藤、蜈蚣、冬虫夏草、红花、人参、银杏叶等度值较高，与AS相关免疫细胞浸润的关系最密切。药物归经多归为肝、肾、肺经，如图9A所示；功效多与清热补虚及祛风湿相关，如图9B所示。表2为AS免疫细胞浸润相关生物途径的中药复方预测结果，涉及补肾活血汤、化痰逐瘀散结汤、寄生汤等多种方剂，功效多为补肾健脾、活血通络，兼祛风散寒化痰。

表2 AS 免疫细胞浸润相关的中药复方预测

图8 AS 免疫细胞浸润相关的中药预测网络图

图9 AS 免疫细胞浸润相关的中药归经功效统计

3 讨论

AS是一种多种复杂因素诱发的疾病，随着现代医学的进步，AS的诊断与治疗已获得了飞速的发展，但由于其病因的不确切性，目前尚无治愈的疗法[14]。因此深入研究AS的发病机制，为AS的治疗提供潜在的靶点对于AS的诊治具有重要意义。本研究通过GEO数据库筛选合适的基因芯片，比较正常人群与AS患者血清样本的基因表达谱，得到差异基因并通过蛋白网络互作获得其核心靶基因。核心靶基因分别是SNRPG、LSM3、PSMA6、MRPL22、RPL26、COX7C、PSMA4、RPS17、RPL31、SLIRP、NDUFS5、RPS24、NDUFB3、SRP14。有研究表明SNRPG与肠道中巨噬细胞活化相关[15]，而产生巨噬细胞迁移抑制因子（MIF）的巨噬细胞在AS患者肠道中富集，MIF在AS患者的关节液中处于高水平，由此推测SNRPG干预AS的机制可能与激活肠道中MIF相关。LSM3基因在AS患者中表达下降，LSM3基因作为预催化的前体参与前体mRNA剪接，而前体mRNA在加工和剪接过程中可能会影响AS患者其他基因的表达[16]。PSMA4、PSMA6基因是蛋白酶体功能所需的20S蛋白酶体核心，蛋白酶体复合物是一种蛋白水解系统，通过主要组织相容性复合物在AS适应性免疫系统的功能中发挥关键作用[17]。MRPL22在细胞质中生成，属于基因编码线粒体核糖体蛋白的成员，对线粒体氧化磷酸化至关重要，且在AS中调节凋亡诱导因子方面发挥关键的作用[18]。RPL26通常构成参与核糖体亚基翻译的细胞过程，而敲低该基因的表达，能使成纤维细胞样滑膜细胞中的NF-κB信号通路去磷酸化，进一步降低AS相关风险基因的转录[19]。COX7C是细胞色素氧化酶（COX）的组成部分[20]。这对能量的生产至关重要，与本研究中的KEGG富集分析中的核糖体生物发生、氧化磷酸化密切相关，并且MRPL22、RPS17、RPL31、NDUFS5、RPS24等核糖体相关蛋白可能与AS发病的能量代谢有关。而其他基因尚缺乏与AS相关的研究，可能是AS的潜在治疗靶点，需进一步发掘验证。本研究中差异基因蛋白互作中的部分枢纽基因有AS的相关基础研究的结果支持，进一步说明了本研究预测结果的可靠性，但其他基因的表达或功能在AS的发病基础研究中鲜有提及。这些基因可能为AS诊断和治疗靶点的研究提供新的思路。

本研究进一步对AS相关差异基因进行GO及KEGG富集分析。结果显示，差异基因的生物学过程不仅有T细胞介导细胞毒性反应参与，还有更为丰富的中性粒细胞活化介导的相关免疫反应，包括中性粒细胞趋化性、中性粒细胞脱粒等生物过程。目前现代医学认为，AS是一种由T淋巴细胞介导的慢性炎症性疾病，尤其是CD4+T细胞的紊乱在其病程的发展中发挥关键性作用，最近的研究表明，先天免疫细胞在炎症和新骨形成过程中也很重要，而中性粒细胞越来越被认为是自身炎症和自身免疫性疾病的介质[21]，且中性粒细胞通过诸如TLR受体、吞噬、脱颗粒和活性酶的释放来实现上述过程，即细胞以网状物的形式释放其染色质[22]。NF-κB介导的免疫调节在AS中发挥着枢纽作用，Th1/Th2的平衡在免疫调节中发挥重要作用。I-κB激酶的激活使得NF-κB的抑制作用解除，从而激发不同Th因子的合成和释放，进而调节免疫紊乱[23]。白细胞介素-1（IL-1）介导的信号通路在AS差异基因的富集上也表现出良好的聚集性。IL-1由炎症滑膜中活化的巨噬细胞分泌，能启动免疫细胞的募集和炎症反应。有研究显示，IL-1多态性与AS易感性相关，IL-1基因簇中的IL-1A、IL-1B和IL-1RN可触发AS先天性免疫反应[24]。除此之外，AS的差异基因KEGG通路富集中涉及核糖体生物发生、氧化磷酸化等相关通路，与代谢、炎症反应相关通路联系密切。

本研究利用CIBERSORT反卷积算法探索免疫细胞浸润在AS发病机制的作用，结果显示，AS组患者血清中单核细胞较正常组明显增多，而调节性T细胞（Treg）则明显减少。Treg主要分为胸腺来源的Treg细胞（tTreg）及诱导调节性T细胞（iTreg），其是CD4+T细胞的一小部分，在维持免疫耐受和自身免疫方面起着关键作用[25]。Treg在AS中的含量变化是有争议的，但多数学者认为Treg比例下降是诱发AS重要因素，其机制可能为AS患者Treg细胞中STAT5磷酸化相对较少，FOXP3的平均荧光强度降低，导致Treg细胞不能有效抑制幼稚T细胞的增殖，进而导致疾病的进展[26]。单核细胞是专业抗原提呈细胞的前体。有研究显示，AS患者的血单核细胞的微阵列分化标记物上调，与正常单核细胞相比，AS单核细胞参与了多种炎症途径的蛋白质水平升高[27]。免疫浸润相关性分析结果中显示，CD8+T细胞与活化的自然杀伤细胞呈中度正相关性。在一项AS患者外周血淋巴细胞亚群变化规律的研究中，其NK细胞含量伴随T细胞显著升高，表明AS的免疫亢进伴有非特异性[28]。γδT细胞介于固有免疫与获得性免疫之间，在免疫反应中发挥重要作用。本研究中，CD4+记忆性T细胞与γδT细胞呈中度正相关，说明两者在AS的发病中发挥协同作用。而CD8+T细胞与中性粒细胞呈中度负相关性的基础研究尚缺乏大规模的验证，故本预测具有一定程度的参考意义。

AS在中医学属“大偻”“偻痹”等范畴。病机为五脏内伤、外感六淫，又以肝肾亏虚、督脉虚空为要，病性多为本虚标实[29]。中医药治疗AS疗效稳固且持久，对于调节患者免疫失衡具有独特的优势。本课题组通过Coremine Medical及本草组鉴预测具有潜在治疗作用的中药，其中雷公藤、蜈蚣、冬虫夏草等出现的频数较高且度值较高，提示与AS相关免疫浸润的关系最密切。药物归经多归肝、肾、肺经，功效多与清热补虚及祛风湿相关，此与AS的病机相符。其中雷公藤中的主要成分雷公藤甲素在体外可减少胶原形成、碱性磷酸酶活性和钙沉积，其对AS的免疫调节作用已经得到证实。雷公藤甲素可能通过上调外周血中的CD4+CD25+CD127+Treg和下调IL-17来实现免疫调节作用[30-31]。蜈蚣在骨痹的治疗中使用频率较高，因其走窜力强，脏腑经络及气血不通之处皆能通之，并且蜈蚣在关节炎的局部炎症细胞浸润中表现为显著的抑炎反应及软骨的保护作用[32]。冬虫夏草素目前尚缺乏在AS免疫调节方面的基础研究，但其在类风湿关节炎小鼠模型中有着显著的平衡免疫的功能，具体机制为调控MMP-3、TIMP-1而抑制炎症细胞因子IL-1β、IL-10，进而平衡Ⅱ型胶原免疫紊乱[33]。在方剂预测中，补肾活血汤、寄生汤中多数药物与预测中药相符合。有临床研究表明，补肾活血汤能显著改善AS患者的疼痛、晨僵不适及活动性炎性指标[34]，且有报道称补肾活血汤含药血清对CD4+CD25+Treg增殖具有明显促进作用[35]。而在本研究中，枢纽差异基因SNRPG、LSM3、PSMA6、MRPL22与Tregs呈患者负相关，但相关性较弱。故笔者推测，上述靶基因可能作为补肾活血汤干预AS免疫浸润相关机制的靶点，但需进一步实验验证。

本研究局限之处在于数据分析来源依赖于GEO数据库中所纳入的AS探针样本，需要进一步实验验证，而单一芯片纳入样本的数量局限可能造成偏倚。总之本课题组利用生物信息学方法对AS的免疫浸润机制进行分析并预测相关生物途径所涉及的中药及复方，对于后续研究及临床用药具有参考意义。