克罗恩病(CD)是一种胃肠道慢性炎症性疾病，主要累及末端回肠和邻近结肠，临床症状为腹痛、腹泻、瘘管、肛门病变、直肠出血、发热、体质量减轻和疲乏等^[1-2]。研究^[3]显示，中国CD的发病率约为0.15%。CD的发病涉及多因素，可能是由环境因素、免疫系统、易感基因和肠道微生物群变化等相互作用引起。研究^[4]显示，核苷酸寡聚结合域蛋白(NOD2)基因位点功能缺失突变是对CD影响最大的风险等位基因，对携带NOD2风险等位基因的CD患者进行直接离体分析，发现NOD2的缺失导致活化的成纤维细胞和巨噬细胞的稳态失调。研究^[5]表明，携带岩藻糖基转移酶2(FUT2)变异体的人会减少抗原的分泌，与细菌的相互作用会改变，更容易发生CD。黏蛋白2(MUC2)的缺乏可以改变黏液的组成，从而影响CD发生^[6]。目前尚无能够彻底治愈CD的特效药物或统一的治疗方案，其临床诊断具有一定的挑战性，并且已发现患者会出现严重的诊断延迟，这可能导致较差的临床结果^[7]。因此，探索CD可能的分子机制及相关基因，有助于提高CD患者的诊断准确性和治疗效率，为筛选出CD药物新靶点和免疫联合治疗提供新思路。

线粒体是关键的细胞器，主要功能是氧化磷酸化产生三磷酸腺苷(ATP), 以维持细胞的生物能量。线粒体功能异常可能会导致多种疾病的发生，其中包括CD。肠道通透性是CD的标志之一，与线粒体功能障碍相关^[8]。线粒体功能障碍在癌症、神经退行性疾病等多种慢性炎症性疾病中起着重要作用^[9]。线粒体活性氧影响稳态信号通路，以控制细胞增殖和分化，并有助于适应性应激信号通路^[10]。线粒体功能障碍被证明在CD发病机制中起着重要作用，线粒体功能障碍引发ATP合成受阻，进而导致能量供应短缺，影响肠道功能的正常运作，具体表现为肠细胞分化过程的紊乱、紧密连接稳定性的削弱，以及丁酸盐氧化代谢的受阻，共同对肠道健康构成了威胁。线粒体功能障碍还通过活性氧生成的增加和线粒体自噬的破坏，导致有缺陷的线粒体在体内累积，这个过程伴随着内部抗原成分的积累和炎症反应的激活^[11]。特别是在潘氏细胞中，线粒体功能对于维持细胞功能至关重要。线粒体功能障碍可导致功能失调的潘氏细胞与肠道干细胞分化，导致CD复发。本研究使用生物信息学和机器学习方法鉴定CD的线粒体相关基因，通过筛选的相关基因构建人工神经网络诊断模型，为CD临床诊断与免疫治疗提供新视角。

1.   材料与方法

1.1   数据获取

从高通量基因表达(Gene Expression Omnibus, GEO)数据库检索并获取CD数据集GSE186582(正常样本25例, CD样品196例), GSE102133(正常样本12例, CD样本65例), GSE95095(正常样本12例, CD样本24例)。其中GSE186582、GSE102133作为训练数据, GSE95095作为验证数据，见表 1。通过Mito Carta 3.0^[12]数据库获取人类线粒体相关基因1 136个。

表  1  选择GEO数据集的信息

GEO数据集	平台	CD样本/例	对照样本/例	组织来源	属性
GSE186582	GPL570	196	25	回肠	训练集
GSE102133	GPL6244	65	12	回肠	训练集
GSE95095	GPL14951	24	12	CD受累组织	验证集

下载: 导出CSV 

| 显示表格

1.2   差异表达基因(DEGs)筛选

将GSE186582、GSE102133数据合并，使用R语言中的“SVA”包(版本3.48.0)对其进行批次校正后作为训练集。使用R语言中的“limma”包(版本3.56.2)比较CD患者组和正常组中各个基因的表达水平，通过lmFit函数拟合线性模型model.matrix[~0+factor(Type)], 使用makeContrasts函数创建疾病组和对照组的对比矩阵定义并将其应用到之前拟合的线性模型上，使用经验贝叶斯方法调整基因表达差异的估计，最后通过错误发现率(FDR)方法对P值进行矫正，从而筛选出DEGs。筛选标准为矫正后P<0.05且|logFC|>1。

1.3   DEGs分析

通过在线工具Metascape分析筛选出的DEGs。在Metascape中选择基因所属的物种为人类。Metascape分析中选择Express Analysis方法，其中结合了多个数据库进行分析，包括KEGG Pathway、GO Biological Processes、Reactome Gene Sets、Canonical Pathways、CORUM、WikiPathways、PANTHER Pathway。Metascape中过滤条件为: P<0.01且富集因子﹥1.5;在功能或通路中至少存在3个差异基因。再使用R语言的“clusterProfiler”包(版本4.8.3)和“org.Hs.eg.db”包(版本3.17.0)进行基因本体论(GO)功能分析和京都基因与基因组百科全书(KEGG)通路富集分析^[13], 进一步研究DEGs的生物学作用。在GO、KEGG分析中使用了人类的基因注释库(org.Hs.eg.db)。筛选标准为矫正后P<0.05。

1.4   机器学习筛选CD特征基因

使用R语言“VennDiagram”包(版本1.7.3)对DEGs和线粒体基因取交集，得到CD线粒体相关的DEGs。通过最小绝对收缩和选择算子(LASSO)回归算法及随机森林算法(RF)进一步筛选CD的生物标志物。LASSO回归是一种统计技术，可用于研究临床变量对结果预测的影响。LASSO回归选择性地突出重要的预测变量，并通过将其系数缩小到0来有效删除不太重要的预测变量^[14]。RF是一种基于多种分类树构造的集成学习方法，通过构建多个决策树并集成他们的预测结果来提高预测精度和泛化能力^[15]。LASSO回归使用R语言中“glmnet”包(版本4.1.8)筛选。RF采用R语言中“randomForest”包(版本4.7.1.1)筛选，纳入基尼系数法计算重要性评分>2的基因。最后通过R语言“VennDiagram”包对LASSO和RF算法的结果基因取交集，以更好地筛选出CD的特征基因。

1.5   人工神经网络构建

先通过CD线粒体特征基因的表达水平对其进行评分。基因的|logFC|﹥0, 则其表达量大于中位值时得分为1, 小于中位值得分为0;基因的|logFC|<0, 表达量小于中位值时得分为1, 大于中位值得分为0^[16]。再通过R语言中“neuralnet”包(版本1.44.2)和“NeuralNetTools”包(版本1.5.3), 根据训练数据集的样本属性和所筛选出的特征基因的基因评分构建人工神经网络诊断模型，通过基因的表达量与基因权重对疾病组和对照组进行分类验证。模型中隐藏层参数设置为5, 输出层通过隐藏层节点与权重绘制，最终得到人工神经网络模型。采用R语言“pROC”包(版本1.18.4)绘制模型的受试者操作特征(ROC)曲线，评估模型对疾病预测的准确性。ROC的曲线下面积(AUC)﹥0.7表示结果较好。采用验证集数据构建ROC曲线，验证模型准确性。

1.6   免疫细胞含量分析

CIBERSORT作为目前被广泛使用的免疫细胞浸润算法，其卓越之处在于运用线性支持向量回归的方法来解析组织表达矩阵。该方法能够精确地量化每个样本中22种免疫细胞的丰度分数，从而提供了对免疫反应的深入理解^[17]。获取CD组和正常组中的免疫细胞表达含量，使用“corrplot”包(版本0.92)、“ggpubr”包(版本0.6.0)以热图的形式绘制免疫浸润细胞相关性。免疫细胞之间的相关性判断标准为免疫细胞之间相关系数的绝对值越接近1，表示2个细胞之间的相关性越强; 越接近0, 表示相关性越弱。对筛选出的生物标志物进行免疫细胞浸润分析。生物标志物与免疫细胞之间的相关系数绝对值越高，表示二者的相关性越强。

2.   结果

2.1   CD的DEGs筛选

对数据集GSE186582和GSE102133共261例CD样本以及37例正常对照样本，以矫正后P<0.05和|logFC|﹥1为过滤标准进行差异分析，获得551个DEGs。在疾病组中, 275个基因的表达明显上调, 276个基因的表达明显下调。前30个DEGs利用R语言“pheatmap”包(版本1.0.12)、“ggplot2”包(版本3.4.4)绘制DEGs的热图，并用火山图对DEGs可视化。见图 1。

图  1  CD的DEGs差异分析

A: DEGs热图，红色代表高表达，蓝色代表低表达; B: DEGs火山图，红色表示上调，绿色表示下调，黑色表示无差别基因。CD: 克罗恩病; DEGs: 差异表达基因。

下载: 全尺寸图片 幻灯片

2.2   DEGs分析

通过Metascape分析得到20个聚类(簇), 281条GO或通路，其中每个聚类富集的最为显著的GO或通路为对细菌的反应、编码细胞外基质(ECM)相关蛋白的基因集合、编码核心ECM的基因集合、核受体介导的信号转导通路、有机阴离子输运、对异生物质刺激的反应、激素水平调节、正调控细胞运动性、有机羟基化合物代谢过程、循环系统过程等，见图 2。GO富集分析富集到939个条目，其中富集到735个生物过程(BP), 54个细胞组分(CC)以及149个分子功能(MF)条目。BP、CC、MF前10位的条目显示(图 3A)，在BP中主要包括中性粒细胞迁移、中性粒细胞趋化性、白细胞趋化性等，在CC中主要包括细胞的顶端部分、顶端质膜、刷状缘等，在MF中主要包括ECM结构成分、金属肽酶活性、四吡咯结合等。KEGG富集分析富集到34条通路，表明DEGs主要富集在细胞因子-细胞因子受体相互作用、蛋白质消化吸收、脂质和动脉粥样硬化、趋化因子信号通路、白细胞介素-17(IL-17)信号通路等(图 3B)。

图  2  DEGs的Metascape分析

A: DEGs富集网络图，图中所示为聚类(簇)而非具体的功能或通路，其中的节点代表功能或者通路，颜色以每个聚类中P值最小的功能或通路名称命名; B: DEGs富集柱状图，纵坐标为功能或通路的名称。DEGs: 差异表达基因。

下载: 全尺寸图片 幻灯片

图  3  DEGs的GO和KEGG富集分析

A: GO柱状图; B: KEGG气泡图。GO: 基因本体论; KEGG: 京都基因与基因组百科全书; DEGs: 差异表达基因。

下载: 全尺寸图片 幻灯片

2.3   机器学习筛选CD特征基因

CD的DEGs和线粒体基因取交集得到20个基因(RMND1、AGXT2、MTHFD2、BPHL、PXMP2、GLS、SOD2、NAGS、MAOB、HMGCS2、MAOA、ACSL1、OAT、ARG2、EPHX2、GATM、BCL2A1、AKR1B10、OTC、FABP1)。通过LASSO回归得到CD线粒体特征基因15个(RMND1、AGXT2、MTHFD2、BPHL、PXMP2、SOD2、MAOB、HMGCS2、MAOA、ACSL1、OAT、ARG2、EPHX2、BCL2A1、OTC)。在RF算法中，首先构建500棵决策树，随后通过执行交叉验证来评估模型的误差，并基于评估过程选定了54棵树作为最优数量。利用基尼系数作为度量标准，计算了每个基因在模型中的重要性得分^[16]。最终识别出了11个特征基因(BPHL、PXMP2、RMND1、NAGS、SOD2、AGXT2、MAOA、MAOB、HMGCS2、MTHFD2、AKR1B10), 其重要性评分均超过了2。对LASSO回归算法和RF算法得到的基因取交集，最终得到9个CD特征基因(SOD2、MTHFD2、BPHL、PXMP2、RMND1、AGXT2、MAOA、HMGCS2、MAOB)。见图 4。

图  4  机器学习算法筛选出的关键基因结果

A: DEGs和线粒体基因的相交维恩图; B: 2种算法相交基因的维恩图; C: LASSO选择的特征基因; D: 随机森林树数与模型误差的相关性图; E: 随机森林分类器的基尼系数法的结果图。

下载: 全尺寸图片 幻灯片

2.4   人工神经网络模型构建

构建人工神经网络模型(图 5A), 该模型依据9个特征基因的评分作为输入。网络结构包括1个输入层，由9个神经元组成，对应筛选出的9个特征基因; 1个隐藏层，包含5个神经元，用于处理并学习输入数据之间的复杂关系; 1个输出层，设有2个神经元，用于输出模型的预测结果。在训练组中，该模型对于正常样品和CD样品的识别准确率分别达到了0.973和0.908。采用ROC曲线分析进一步验证人工神经网络模型的可行性和性能。在训练组中，模型的AUC为0.956(图 5B); 在验证组中，模型的AUC为0.736(图 5C)。

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

A: 人工神经网络模型可视化; B: 训练集的ROC曲线; C: 验证集的ROC曲线。

下载: 全尺寸图片 幻灯片

2.5   免疫细胞浸润

基于CIBERSORT算法评估免疫细胞浸润情况。其中，幼稚B细胞与幼稚CD4⁺T细胞呈正相关关系，激活肥大细胞与中性粒白细胞呈正相关关系，静息肥大细胞与激活肥大细胞呈负相关关系，静止记忆CD4⁺T细胞与中性粒白细胞呈负相关关系，见图 6。相关性分析结果显示, SOD2与中性粒白细胞呈正相关(R=0.52, P<0.001), MTHFD2与活性记忆CD4⁺T细胞呈正相关(R=0.60, P<0.001), BPHL与活性树突状细胞呈负相关(R=-0.50, P<0.001), PXMP2与活性记忆CD4⁺T细胞呈负相关(R=-0.56, P<0.001), RMND1与静止记忆CD4⁺T细胞呈正相关(R=0.34, P<0.001), AGXT2与静止记忆CD4⁺T细胞呈正相关(R=0.46, P<0.001), MAOA与静止记忆CD4⁺T细胞呈正相关(R=0.54, P<0.001), MAOB与CD8⁺T细胞呈正相关(R=0.43, P<0.001), HMGCS2与中性粒白细胞呈负相关(R=-0.37, P<0.001)。见图 7。

图  6  免疫细胞表达相对含量相关性分析

下载: 全尺寸图片 幻灯片

图  7  特征基因与免疫细胞浸润相关性

A: SOD2; B: MTHFD2; C: BPHL; D: PXMP2; E: RMND1; F: AGXT2; G: MAOA; H: MAOB; I: HMGCS2。

下载: 全尺寸图片 幻灯片

3.   讨论

CD在世界范围内发病率不断上升，具体病因尚未完全清楚。研究^[18]认为肠黏膜中正常免疫反应失调是CD发病机制的基础。研究^[16]指出线粒体损伤与免疫反应存在紧密联系, CD与免疫系统密切相关，但关于CD与线粒体损伤的关系尚未明确。因此，筛选CD线粒体相关基因并建立相关基因的诊断模型对CD的后续研究具有重要意义。

本研究通过机器学习方法对CD线粒体相关基因进行筛选，得到9个线粒体相关基因(SOD2、MTHFD2、BPHL、PXMP2、RMND1、AGXT2、MAOA、HMGCS2、MAOB), 构建了人工神经网络诊断模型。SOD2基因属于铁/锰超氧化物歧化酶家族，已发现SOD2的遗传变异与许多疾病有关，包括神经退行性变、线粒体功能障碍、过早衰老、血管生成和癌症^[19-20]。自身免疫性疾病的病理生理学潜在机制是自由基生成与自由基清除系统之间的不平衡，而SOD2在清除自由基方面同样起着重要作用。MTHFD2是一种线粒体单碳代谢酶，在肿瘤中高表达，靶向MTHFD2已被用作肿瘤治疗的靶点。研究^[21]表明MTHFD2协调体内巨噬细胞极化，在免疫细胞的命运和功能中起着至关重要的作用。在患有各种炎症和自身免疫性疾病的个体全血的RNA测序数据集中, MTHFD2在多种疾病中呈一致的过表达，包括溃疡性结肠炎、CD等^[22]。研究^[23]表明MTHFD2缺乏可以降低各种炎症状况模型中的疾病程度, MTHFD2抑制剂对扩展到B细胞功能的炎症具有保护作用。BPHL基因是一种分子量为30 kDa的丝氨酸水解酶，抑制BPHL表达可以降低肿瘤细胞的恶性程度，包括减少增殖、增加细胞凋亡、减少集落形成和减少转移和细胞周期阶段分布改变^[24]。

PXMP2位于过氧化物酶体膜中，参与多种恶性肿瘤的进展，其相关途径包括过氧化物酶体脂质代谢和4-羟脯氨酸降解。RMND1是一种核编码的线粒体蛋白，是线粒体正常功能所需的核基因，其基因产物形成线粒体翻译所必需的高分子量复合物，可能协调线粒体核糖体的组装或维持。RMND1致病变异可引起广泛的线粒体翻译缺陷^[25]。RMND1的基因突变也可能会导致严重的卵巢功能不全。RMND1致病基因会引起多种氧化磷酸化缺陷^[26]。AGXT2是一种多功能线粒体氨基转移酶, AGXT2单核苷酸多态性与糖尿病、轻度认知障碍和生物学参数(血压、血糖和肌酐)等有关^[27]。MAOA、MAOB是蛋白质编码基因。MAOA的相关途径包括神经递质清除和细胞色素P450氧化。MAOA在肿瘤生长、迁移和转移中起着至关重要的作用。抑制MAOA的表达可以降低癌细胞的侵袭性^[28]。MAOB是一种主要在星形胶质细胞中表达的线粒体酶，可导致阿尔兹海默症大脑中的神经元功能障碍和神经退行性变^[29]。HMGCS2是负责生酮作用限速反应的酶。研究^[30]表明缺乏生酮作用与一系列炎症性疾病存在相关性。HMGCS2的沉默可以促进细胞活力，抑制细胞凋亡、炎症和氧化应激。在对脓毒性心肌损伤研究^[31]中发现HMGCS2通过刺激巨噬细胞M2极化以减少脓毒症引起的炎症。在CD患者与黏膜愈合相关的肠道mRNA表达谱分析研究^[32]中也表明HMGCS2对CD有良好的预测效果。

免疫细胞浸润分析显示特征基因与静止记忆CD4⁺T细胞、活性记忆CD4⁺T细胞、活性树突状细胞、中性粒白细胞、CD8⁺T细胞等相关。研究^[18]发现树突状细胞与巨噬细胞及其亚群在CD的病理机制中起着重要作用。肠道树突状细胞能够获得管腔抗原，识别和启动免疫反应，在维持肠道免疫耐受方面起着至关重要的作用。研究^[33]表明CD的发病机制基于组织炎症，由针对管腔细菌抗原的不受抑制的免疫反应引起, CD4⁺T细胞、CD8⁺T细胞、B细胞、CD14⁺单核细胞和自然杀伤细胞等免疫细胞在浸润CD患者肠道时参与这一过程。研究^[34]发现CD主要由T淋巴细胞，尤其是CD4⁺T细胞的慢性炎症引起。研究^[35-36]表明中性粒细胞活化与CD疾病活动相关，中性粒细胞活性增加会损害肠黏膜屏障功能，损害肠上皮，并加重疾病症状^[37]。

综上所述，线粒体相关基因SOD2、MTHFD2、BPHL、PXMP2、RMND1、AGXT2、MAOA、HMGCS2、MAOB与CD发生发展相关。CD的发生发展可能与静止记忆CD4⁺T细胞、活性记忆CD4⁺T细胞、活性树突状细胞、中性粒白细胞、CD8⁺T细胞等免疫细胞有关。针对上述免疫细胞进行研究并制定措施可能有助于治疗CD。