基于新发现的乳糜泻免疫相关生物标记物开发神经网络诊断模型和潜在药物

背景：

乳糜泻（Celiac disease，CD）是最常见的肠道炎症性疾病之一，其典型特征是由于摄入谷蛋白引起的自身免疫性疾病。尽管随着时间的推移，CD的发病率和患病率有所增加，但诊断方法和治疗方案仍然有限。因此，迫切需要研究CD潜在的生物标志物和靶向药物。

研究方法：

基因表达数据从GEO数据集下载。差异表达分析用于鉴定失调的免疫相关基因。采用随机森林（randomForest）、SVM-RFE和LASSO等多种机器算法选择核心免疫相关基因（Hub immune-related genes，HIGs）。基于HIGs构建了免疫相关基因评分（Immune-related genes score，IG score）和人工神经网络（Artificial neural network，ANN）。利用Enrichr平台和分子对接方法鉴定了靶向HIGs的潜在药物。

结果：

在全基因组水平鉴定了失调的免疫相关基因，并证明了它们在CD相关免疫途径中的作用。通过对几种机器算法的结果取交集进一步筛选Hub基因（MR1, CCL25和TNFSF13B）。同时，使用单样本基因集合富集分析（ssGSEA）和一致性聚类将CD患者分为免疫活性高或低的不同亚型。通过构建基于HIGs的IG评分，发现IG评分高的患者主要属于于高免疫活性亚组，这表明HIGs和CD患者的免疫活性之间有很强的联系。此外，新构建的人工神经网络模型显示了HIGs良好的诊断能力。从机制上讲，本研究证实了HIGs在调节CD的免疫和炎症状态中起着关键作用。通过靶向HIGs，还利用Enrichr平台和分子对接方法发现了潜在的抗CD治疗药物。

结论：

这项研究揭示了HIGs，并阐明了这些基因在 CD 中调控的网络，突出了HIGs在准确预测患者是否患有CD中的重要意义。因此，这项研究为CD的诊断生物标志物和治疗靶点的开发提供了广阔的前景。

1.数据获取

1.1  从GEO数据库（https://www.ncbi.nlm.nih.gov/geo/）中筛选并下载了两个独立的CD数据集GSE11501和GSE164883。GSE11501数据集包含来自英国的110名CD患者和22名健康对照。利用GPL6104平台标注信息将探针转化为相应的基因符号。GSE164883数据集包含来自德国的26例CD患者和22例健康对照。利用GPL10558平台标注信息将探针转化为相应的基因符号。GSE11501数据集作为训练集，GSE164883数据集作为验证集。

1.2  然后，从Immport数据库（https://www.immport.org/resources）下载了免疫相关基因列表，共计2483个免疫相关基因。我们从TISIDB数据库（http://cis.hku.hk/TISIDB/）下载了28个免疫基因集。

2.结果

2.1 研究流程图：

Step1: 鉴定差异表达的免疫相关基因以及对这些基因进行功能富集分析

Step2: 基于多种机器学习算法鉴定hub基因

Step3: 使用ssGSEA和一致性聚类将CD患者分为免疫活性高和免疫活性低的不同亚型

Step4: 基于HIGs构建了免疫相关基因评分和人工神经网络

Step5: 利用Enrichr平台和分子对接方法鉴定了靶向HIGs的潜在药物

2.2 CD中差异表达免疫相关基因的鉴定

  通过差异表达分析鉴定了了896个差异表达基因，其中369个上调基因，527个下调基因（图2A）。为了在差异表达基因中进一步筛选免疫相关基因，对896个差异表达基因与ImmPort数据库（https://www.immport.org/resources）获得的2483个免疫相关基因取交集，最终获得了58个差异表达的免疫相关基因（图 2B）。差异表达的免疫相关基因的表达谱如图2C所示。

  为了研究这些基因的作用，我们进行了GO和KEGG分析。因此，GO富集分析揭示了多种免疫相关的生物学过程，包括“T细胞和B细胞活化”、“T辅助细胞分化”、“白细胞介素-6/-8产生”、“白细胞介素-15介导的信号通路”、“MHC II类蛋白复合物结合”、“抗原加工和通过MHC II类呈递外源肽抗原”、“MyD88依赖性toll样受体信号传导途径”、“G蛋白偶联受体结合”、“肽基-酪氨酸自磷酸化”、“调节炎症反应中涉及的细胞因子产生”、“调节对抗原刺激的急性炎症反应”、“趋化因子（C-X-C基序）配体2产生”和“调节 NLRP3 炎症小体复合物的组装”。KEGG富集分析表明差异表达的免疫相关基因也富集到了多种免疫相关通络中，包括“T细胞受体信号传导途径”、“B细胞受体信号传导途径”、“Th1和Th2细胞分化”、“Th17细胞分化”、“自然杀伤细胞介导的细胞毒性”、“EGFR酪氨酸激酶抑制剂抗性”、“IgA产生的肠道免疫网络”和“癌症中的PD-L1表达和PD-1检查点途径”（图S1）。

2.3 CD患者被分成具有不同免疫活性的亚组

  乳糜泻是一种自身免疫性疾病，其中免疫系统在其发病机制中起着核心作用，这表明CD患者的免疫活性与CD的发展密切相关。因此，为了评估CD患者的免疫活性和CD风险之间的关系，进行了ssGSEA和一致性聚类。我们对TISIDB数据库（http://cis.hku.hk/TISIDB/）鉴定的28个免疫基因集进行了ssGSEA。计算了每个CD患者的28个免疫基因集的ssGSEA评分。基于28个免疫基因集的ssGSEA评分，使用一致性聚类的k-均值对GSE11501数据集检索到的CD患者进行分类。然后确定了最佳聚类数为K = 2。在k = 2时，CD患者被分为两个亚组，C1和C2，具有明确的界限，表明聚类稳定可靠。随后，进行UMAP和t-SNE分析以验证亚型分配，两种方法的结果表明一个亚组中的样本比另一个亚组中的样本更相似（图3D，E）。

  通过比较C2亚组和C1亚组之间的28个免疫基因组的ssGSEA评分时，发现15个免疫细胞存在显著差异。值得注意的是，对于免疫细胞，如“中枢记忆CD8 T细胞”、“自然杀伤细胞”、“自然杀伤T细胞”、“活化树突细胞”、“浆细胞样树突细胞”、“巨噬细胞”、“嗜酸性粒细胞”、“肥大细胞”和“嗜中性粒细胞”，C1亚组表现出明显更高的ssGSEA评分（图S4）。此外，我们通过比较C1和C2亚组各自的ssGSEA评分，检测了它们之间的免疫活性关系。我们的研究结果表明，与C2亚组相比，C1亚组的ssGSEA评分明显更高（图3F），这表明C1亚组的整体免疫活性超过了C2亚组。因此，基于他们的免疫特征，我们能够将CD患者分为两个不同的亚组:高免疫活性组和低免疫活性组。

2.4 基于HIGs的IG评分的构建与验证

  基于上面提到的HIGs（MR1, CCL25和TNFSF13B），计算IG评分，并进一步使用主成分分析来评估从GSE11501中检索的每个样本发生CD的风险。根据IG评分<或>0，将从GSE11501中检索的所有样本分为低IG和高IG亚组。同时，我们进行了相关分析，以评估IG评分与免疫活性的相关性。通过冲积图中每个样本的基本概况，其中包括高和低IG分数组以及由上述ssGSEA划分的高和低免疫活性患者。我们的研究结果显示，大多数IG评分高的CD患者属于C1亚组，而大多数IG评分低的CD患者属于C2亚组（图4A）。随后，我们还比较了三个亚组（包括健康对照组、C1亚组和C2亚组）的IG评分。结果显示，C1亚组患者IG评分高于C2亚组，CD患者（C1和C2） IG评分高于健康对照组（图4B）。此外，在进行相关分析时，我们发现IG评分与ssGSEA评分呈正相关（图4C）。我们进行ROC分析以进一步评估IG评分的预测能力。由于IG评分的ROC曲线下面积（AUC）为0.758，结果表明IG评分具有良好的可预测性能（图4D）。综上所述，这些数据不仅表明基于HIGs的IG评分可以预测个体发生CD的风险，而且还建立了HIGs与CD患者免疫活性之间的紧密联系。

  此外，我们评估了基于HIGs的IG评分在另一个独立的CD患者队列中是否具有相似的预测值。将GSE164883的CD患者纳入IG评分分析。同样，我们发现CD患者的IG评分高于对照组。AUC值为0.886的分数也显示了IG评分预测的性能（图S5）。这些结果再次表明，基于HIGs的IG评分具有较好的敏感性和预测性。

2.5 HIGs在CD的炎症状态中起关键作用

  为了进一步研究HIGs的调节作用，我们进行了GSEA分析以确定HIGs调节的信号通路。结果表明CCL25在“B细胞受体信号通路”、“结肠直肠癌”、“胃癌”、“IL-17信号通路”和“人T细胞白血病病毒1感染”中显著富集。TNFSF13B与“TNF信号通路”、“IL-17信号通路”、“炎症性肠病”、“T细胞受体信号通路”和“抗原加工和呈递”显著相关。MR1与“B细胞受体信号传导途径”、“原发性免疫缺陷”、“病毒性致癌作用”、“结肠直肠癌”和“Th1和Th2细胞分化”显著相关（图5A）。总之，我们发现HIGs可以调节免疫细胞受体信号通路，如B细胞或T细胞受体信号通路。其中，CCL25和TNFSF13B可以调节IL-17信号通路，特别是作为促炎细胞因子的IL-17，其在CD发展期间促进慢性炎症状态特征。TNFSF13B可调节TNF信号通路，并可能引发炎症性肠病。MR1能调节Th1和Th2细胞分化，产生丰富的促炎细胞因子。

  为了揭示HIGs的调节作用，我们进一步整合了KEGG途径分析和文献报道，并在Figdraw平台中绘制了HIGs在CD发展中的主要参与关系（图5B）。我们发现CCL25和TNFSF13B在HSA 04672-肠道免疫网络中被富集以产生IgA。CD发病机制中活化的T辅助细胞将活化B细胞，以诱导它们产生针对组织转谷氨酰胺酶、醇溶蛋白和内膜的IgA抗体，这加剧了CD的慢性炎症状态特征。还表明MR1在树突状细胞（DCs）上表达，树突状细胞向粘膜相关不变T （MAIT）细胞呈递抗原，导致MAIT活化和产生促炎细胞因子，如肿瘤坏死因子（TNF）-α、干扰素（INF）-γ、白细胞介素（IL）-23、IL-17。随后，这些介质将募集并激活其他免疫细胞，导致肠道组织炎症的诱导。

2.6 人工神经网络模型的构建与验证

  为了验证HIGs的诊断作用，我们在GSE11501训练集和GSE164883验证集的诊断测试评估中检测了每种HIGs的receiver operating characteristic ROC曲线。在AUC值MR1=0.696, CCL25=0.860, TNFSF13B=0.839的情况下，这些HIGs都能在GSE11501训练集中以更高的准确率区分CD和健康对照。除了在GSE11501训练集，还在GSE164883验证集中计算了AUC值。当AUC值为MR1=0.741, CCL25=0.906, TNFSF13B=0.988时，这些HIGs在GSE164883验证集中能够以更高的准确率区分CD和健康对照（图S6）。此外，我们还构建了基于HIGs的神经网络来诊断CD的发病。

    ANN作为人工智能的一种突出形式，广泛应用于临床医学的各个专业领域。值得注意的是，许多研究强调了人工神经网络在诊断和治疗胃肠疾病方面的显著潜力。根据这些发现，我们将HIGs整合到一个人工神经网络框架中，以开发一个能够辨别样本是属于健康对照组还是CD组的预测模型。人工神经网络模型包括三个基本组成部分:输入层，隐藏层，输出层（图6A）。随后，我们将人工神经网络模型生成的预测与样本的实际分组信息进行了比较。ANN对训练集和验证集的预测准确率如表2所示，分别为0.9146和0.92。此外，我们采用 ROC分析来评估人工神经网络模型在训练集和验证集上的预测能力。训练集的曲线下面积（AUC）值为0.793（图6B），验证集的AUC值为0.821（图6C）。总之，人工神经网络模型显示出巨大的潜力，可以作为CD的独立诊断预测器。

2.7 靶向HIGs的潜在药物的识别和对接

  为了找到靶向HIGs的药物，我们使用了enrichment平台（https://maayanlab.cloud/Enrichr/）进行在线分析和筛选。基于DSigDB数据库，鉴定出7种靶向HIGs的药物，p值< 0.05。接下来，我们利用分子对接方法（MDM）研究药物与其靶向HIGs的结合亲和力，其结合能如表4所示。结果表明，Tetradioxin能够靶向所有的HIGs蛋白，且三种HIGs蛋白与Tetradioxin的结合能绝对值均高于与其他药物分子的结合能绝对值（图7A-C），依次为CCL25-Tetradioxin （-5.6 kcal/mol）、MR1-Tetradioxin （- 6.8 kcal/mol）、TNFSF13B-Tetradioxin（-6.62 kcal/mol）。其中，CCL25的氨基酸结合位点LEU-43和TRP-47与Tetradioxin的距离最近，分别为3.46 Å和3.28 Å，而PRO-68、LYS-69、LEU-89和LEU - 90与Tetradioxin的距离相对较远。TNFSF13B的氨基酸结合位点TYR-206、ASN-235和LEU-240与Tetrtradioxin的距离均较远，只有PRO-237与Tetrtradioxin的距离较近，为3.33 Å。相比之下，在MR1的7个氨基酸结合位点中，PHE-30与Tetradioxin有3个氢键形成相互作用，LEU-32与Tetradioxin的距离为3.47 Å， ILE-45与Tetradioxin的距离为3.55 Å， GLN-115与Tetradioxin的距离为3.00 Å， ALA135与Tetradioxin的距离为3.74 Å。总体而言，MR1与四曲霉毒素之间最近的氨基酸结合位点（GLN-115 3.00 Å）的距离小于CCL25（TRP47 3.28 Å）和TNFSF13B （PRO-237 3.33 Å）的最近距离。因此，从HIGs氨基酸结合位点与四曲霉毒素的距离来看，MR1与四曲霉毒素之间存在多个且距离较近的氨基酸结合位点，这些氨基酸可以形成更多的相互作用，即MR1相对于CCL25和TNFSF13B表现出更强的结合能，为6.8 （kcal/ mol）。此外，CCL25-CROTONALDEHYDE （- 3.03 kcal/mol）、MR1-DMBA（- 6.8 kcal/mol）和TNFSF13B-FENRETINIDE （- 6.56 kcal/mol）这三个配合物的结合能绝对值相对较高（图7D-F），而 MR1- demecoline （- 5.94 kcal/mol）和MR1-cyclophosphamide （- 4.64 kcal/mol）以及TNFSF13B- demecoline （- 6.01 kcal/mol）和TNFSF13B- diuron （- 5.4 kcal/mol）这四个配合物的结合能绝对值相对较低，提示这些药物可能对HIGs有调节作用。