高松林,覃 雁,张 鹏,谢兰芳,徐 哲,梁 菲,黄贵华
1.广西中医药大学,广西 南宁 530000
2.柳州市中医医院,广西 柳州 545000
3.广西中医药大学第一附属医院,广西 南宁 530000
溃疡性结肠炎(ulcerative colitis,UC)是一种以结肠、直肠炎症为特征的慢性非特异性炎症性疾病[1]。UC 是炎症性肠病(inflammatory bowel disease,IBD)的一种,在过去的30 年里,IBD 在我国的年发病率逐渐上升,达到3.01/10 万,年患病率为47.06/10 万[2]。2023 年全世界UC 的患病率估计为500 万例,并且发病率呈上升趋势[3]。UC 根据内镜下表现和临床症状可分为活动期和缓解期,活动期UC 患者常有腹痛、腹泻、黏液脓血便等典型临床症状,是临床治疗的重点和难点[4]。UC 是结直肠癌的高风险致病因素之一,UC 患者发生结直肠癌的风险是普通人群的2~3 倍[5-6]。目前,活动期UC 的治疗方案旨在调节免疫系统和肠道炎症以缓解症状,主要使用柳氮磺吡啶、糖皮质激素、免疫抑制剂和其他生物制剂进行治疗,但是疗效有限、副作用明显、复发率高或价格昂贵[7-10]。因此,迫切需要发现新的、更安全的和更有效的药物。
铜是维持人体健康的重要微量元素,细胞内铜浓度保持在一个相对较低的范围内,适度增加可引起细胞毒性,甚至导致细胞死亡[11-12]。铜死亡是一种由铜过量引发的新型细胞死亡方式,在这个过程中铜离子直接与线粒体呼吸中三羧酸循环中的脂肪酰化组分结合,导致脂肪酰化蛋白聚集和铁硫簇蛋白丢失,引起蛋白质毒性应激,最终导致细胞死亡[13]。UC 患者中血清铜较正常人明显升高[14]。铜代谢中含结构域Murr1 蛋白(copper metabolism Murr1 domain,COMMD)1 是COMMD 蛋白家族中第1 个确定的成员,对铜稳态起着重要的作用,COMMD1 的缺失会导致肝细胞铜累积[15-16]。COMMD1 的表达或功能降低可能与UC 的发病机制有关[17]。COMMD1 是核因子-κB(nuclear factor kappa B,NF-κB)的负性调节因子,在控制肠道炎症和抑制结肠炎相关癌症的进展方面具有重要作用[18]。多篇研究报道铜死亡在类风湿性关节炎、系统性红斑狼疮、银屑病等多种免疫系统疾病中发挥着重要的作用,UC 与这些疾病具有类似的发病机制[19-22]。目前尚缺乏活动期UC 中铜死亡相关的研究,本研究通过生物信息学方法探讨活动期UC 中铜死亡相关基因的潜在调控机制和诊断生物标志物,并预测可以通过调节铜死亡治疗活动期UC 的中药,旨在为活动期UC 的诊断和个性化治疗提供新思路。
从GEO 数据库(http://www.ncbi.nlm.nih.gov/geo)下载活动期UC 的表达数据集(GSE11223、GSE87466、GSE87473 和GSE92415)。GSE11223 数据集中有132 个结肠样本,包括63 例活动期UC 患者样本和69 例健康受试者样本(对照),作为训练集。GSE87466 和GSE92415 数据集均有87 例活动期UC 患者样本和21 例健康受试者样本,GSE87473有106 例活动期UC 患者样本和21 例健康受试者样本,三者作为验证集。使用“limma” R 包对GSE11223、GSE87466、GSE87473 和GSE92415 数据集中的数据进行标准化处理,包括注释和校正。从已发表的文献中获取与铜死亡相关的基因(cuproptosis-related genes,CRG),获得了19 个基因[23-25]。技术路线流程如图1 所示。
图1 技术路线流程Fig.1 Technical route flow
以P<0.05 为筛选条件,运用“limma” R 包和Wilcox 检验分析GSE11223 数据集中活动期UC 和对照组中差异表达的铜死亡相关基因(differentially expressed cuproptosis-related genes,DECRGs)。运用“pheatmap”和“ggpubr” R 包对结果进行可视化。为了解铜死亡基因调节因子之间的相关性,运用“corrplot”和“circlize” R 包对DECRGs 进行相关性分析并展示。
运用CIBERSORT 算法计算活动期UC 组和对照组中的免疫细胞浸润程度,P<0.05 表示具有统计学意义,应用“ggpubr” R 包构建箱线图。应用“limma”和“corrplot” R 包分析DECRGs 和免疫细胞之间的相关性,并运用“ggplot2” R 包绘制相关性热图。
“ConsensusClusterPlus”软件包用于计算无监督聚类的数量及其稳定性[26]。为了分析活动期UC 患者中与DECRGs 相关的潜在亚型,根据DECRGs 的表达量运用“ConsensusClusterPlus” R 包进行无监督共识聚类分析,使用k-means 算法将63 个活动期UC 样本分类为集群,k值定义为从1 到9 以生成不同的子类型。通过判断不同k值下的共识聚类、矩阵累积分布函数(uniform clustering cumulative distribution function,CDF)曲线和一致性聚类得分,得到最可靠的聚类分析结果。使用主成分分析(principal component analysis,PCA)对分型的结果进行评估,使用R 包“ggpubr”比较不同亚型中DECRGs 的基因表达谱,P<0.05 代表具有显著性,运用R 包作图进行可视化。
应用CIBERSORT 算法和“limma”“reshape2”“ggpubr”等R 包对不同亚型进行免疫浸润分析,比较不同亚型之间免疫细胞浸润的差异,并以热图和箱型图进行展示。
GSVA 是一种通过表达谱计算基因集富集的非参数、无监督方法,用于识别各种样本群体中生物过程和途径之间的差异[27]。运用“GSVA”和“GSEABase” R 包进行GSVA 富集分析,以阐明不同亚型之间基因富集集的差异。从MSigDB 数据库下载“c2.cp.kegg.v7.4.symbols”文件,运用“limma”包比较不同亚型之间的GSVA 评分来识别基因差异表达的途径,P<0.05 作为判断差异有统计学意义的标准。
基于DECRGs 分型后的基因表达文件,使用“WGCNA” R 包对活动期UC 疾病基因表达矩阵中波动最大的前25%的基因进行WGCNA,筛选分型后GSE11223 数据集的核心模块和基因。使用pickSoftThreshold 函数获得相邻函数中加权参数的最佳值,并用作后续网络构建的软阈值。随后建立加权邻接矩阵,并将其转换为拓扑重叠矩阵(topological overlap matrix,TOM),并通过基于1-TOM 相异矩阵的分层聚类创建基因模块,找出模块间的相似性,将相关检验P值最小的模块确定为活动性UC 分型基因的关键基因模块。
结合“caret”“randomForest”“kernlab”“ggplot2”和“xgboost” R 包构建随机森林(random forest,RF)、支持向量机(support vector machine,SVM)、极限梯度上升(eXtreme gradient boosting,XGB)模型和广义线性模型(generalized linear models,GLM)4 种机器学习模型。提取了WGCNA 中核心基因的表达量,并用4 种模型对结果进行了预测。通过绘制残差箱线图、累积残差分布图和接收器工作特性(receiver operating characteristic,ROC)曲线来评价模型的准确性,以确定最佳机器学习模型。对每个模型基因进行重要性评分,筛选重要性排名前5 的基因作为活动期UC 疾病的特征基因。
提取最佳机器学习模型中5 个疾病特征基因的表达量,使用“rms”和“rmda” R 包构建列线图对每个特征基因进行评分,最后基于综合评分评估患者患病的概率。通过校准曲线评估预期概率和实际结果之间的一致程度,使用决策曲线评估列线图的临床实用性和优势。根据最佳机器学习模型中疾病特征基因的表达量,在GSE87466、GSE87473 和GSE92415 3 个数据集上进行验证,绘制ROC 曲线验证预测结果。
运用“limma”“corrplot”和“circlize” R 包分析疾病特征基因和DECRGs 之间的表达相关性。
将 DECRGs 和活动期 UC 特征基因输入Coremine Medical 数据库(https://www.coremine.com/),以P<0.05 为阈值检索对应的中药。将检索到的中药根据《中国药典》2020 年版对中药名称进行规范化,未收录的中药不进行统计,再导入中医传承计算平台(V3.5)中进行分析。
通过Wilcox 检验确定活动期UC 组中有11 个CRG 的表达与对照组中的表达差异显著。NOD 样受体热蛋白结构域相关蛋白3(NOD-like receptor thermal protein domain associated protein 3,NLRP3)在活动期UC 患者样本中显著上调,而核转录因子E2 相关因子2(nuclear factor erythroid 2-like-2,NFE2L2)、溶质载体家族31 成员1(solute carrier family 31 member 1,SLC31A1)、铁氧还蛋白1(ferredoxin 1,FDX1)、硫辛酸合成酶(lipoic acid synthase,LIAS)、二氢硫辛酰胺脱氢酶(dihydrolipoyl dehydrogenase,DLD)、二氢硫辛酰转乙酰基酶(dihydrolipoyl acetyltransferase,DLAT)、丙酮酸脱氢酶E1 α1 亚基(pyruvate dehydrogenase E1 subunit alpha 1,PDHA1)、谷氨酰胺酶(glutaminase,GLS)、二氢硫辛胺支链转酰基酶E2(dihydrolipoamide branched chain transacylase E2,DBT)和甘氨酸裂解系统H 蛋白(glycine cleavage system protein H,GCSH)在活动期UC 患者样本中显著下调(P<0.05)(图2-A、B)。随后,对这些基因进行了相关性分析发现PDHA1与LIAS、DLAT,DBT与PDHA1,DLD与GCSH、DLAT、PDHA1、DBT,NFE2L2与PDHA1、DBT、DLD,FDX1与DLAT、PDHA1、SLC31A1、DBT、DLD、NFE2L2呈显著的正相关(图3)。NLRP3与其他基因呈负相关,其中与DBT负相关性最显著。
图2 活动期UC 中的DECRGsFig.2 DECRGs in active UC
图3 DECRGs 间的相关性分析Fig.3 Correlation analysis of DECRGs
为了研究活动期UC 的免疫微环境,使用CIBERSORT 算法评估活动期UC 组与对照组之间的免疫细胞浸润水平(图4-A)。与对照组相比,活动期UC 组激活的记忆性CD4+T 细胞(T cells CD4 memory activated)、静息NK 细胞(NK cells resting)、M0 巨噬细胞(macrophages M0)、激活的肥大细胞(dendritic cells activated)和中性粒细胞(neutrophils)上调;静息记忆性CD4+T 细胞(T cells CD4 memory resting)、激活的NK 细胞(NK cells activated)、M2巨噬细胞(macrophages M2)、激活的树突状细胞(dendritic cells activated)下调(图4-B)。
图4 DECRGs 免疫浸润相关性分析Fig.4 Correlation analysis of immune infiltration in DECRGs
对免疫细胞与NLRP3、NFE2L2、SLC31A1、FDX1、LIAS、DLD、DLAT、PDHA1、GLS、DBT和GCSH进行相关性分析,发现初始B 细胞(B cells naive)与LIAS,M1 巨噬细胞(macrophages M1)与DLD、GCSH,M2 巨噬细胞与GCSH,激活的肥大细胞、静息NK 细胞与DLD,γδT 细胞(T cells gamma delta)与FDX1、GLS呈正相关;M0 巨噬细胞与DLAT、FDX1、LIAS,中性粒细胞与LIAS,激活的NK 细胞与DLD,浆细胞(plasma cells)与LIAS,静息记忆性CD4+T 细胞与NLRP3呈负相关(P<0.01)(图4-C)。
根据DECRGs 的表达量对活动期UC 患者进行共识聚类分型,通过综合分析共识聚类、CDF 曲线和一致性聚类得分的结果发现,k=2 时分类最佳(图5-A、B、C、E)。因此,活动期UC 组的63 例患者被分为C1 和C2 2 个亚组。PCA 显示C1 和C2样本之间具有明显的差异,分类效果较好(图5-D)。图5-F 显示了活动期UC 组中2 种亚型中DECRGs表达水平的热图。在C1 和C2 之间观察到不同的CRG 表达谱,NFE2L2、SLC31A1、FDX1、LIAS、DLD、DLAT、PDHA1、GLS、DBT和GCSH在C1中的表达高于C2,NLRP3在C2 中的表达高于C1(图5-G)。
图5 共识聚类分型情况Fig.5 Consensus clustering classification
22 种免疫细胞在不同亚组之间的相对百分比见图6-A。与C1 相比,C2 中M0 巨噬细胞数量更高(图6-B)。GSVA 用于评估不同亚组之间的生物学功能差异,高亲和性免疫球蛋白E 受体(highaffinity immunoglobulin E receptor,FcεRI)信号通路、趋化因子信号通路、T 细胞受体信号通路等通路在C1 亚组中上调。此外,丁酸代谢、萜类骨架的生物合成、脂肪酸代谢和组氨酸代谢等途径在C2 亚组中上调(图6-C)。
图6 亚型间DECRGs 相关免疫浸润分析和GSVA 分析Fig.6 Analysis of DECRGs related immune infiltration and GSVA among subtypes
应用WGCNA 算法根据活动期UC 分型基因建立加权基因共表达网络,以识别与DECRGs 聚类密切相关的关键基因模块。共有9 个重要的共表达模块,蓝色模块与C1 呈最强的正相关(相关系数r=0.85,P= 2×10−18),与C2 呈最强的负相关(r=0.85,P= 2×10−18),见图7。以“geneSig Filter=0.5,moduleSig Filter=0.8”为条件,经Wilcox 检验且P<0.05 表示有统计学意义,从蓝色模块筛选出234 个核心基因。
图7 分型WGCNA 确定活动期UC 相关模块和核心基因Fig.7 WGCNA of genotyping to determine active UC related modules and core genes
为了进一步识别具有高诊断价值的关键标志物,基于蓝色模块234 个核心基因,建立了4 个机器学习模型进行训练。从残差箱线图可以看出GLM、XGB和SVM 的残差在0~0.125,RF 的残差在0.125~0.250(图8-A)。从累积残差分布图得出的结论与上述结果一致(图8-B)。ROC 曲线中XGB 模型的曲线下面积为0.911(图8-C),大于GLM、RF 和SVM,故XGB 模型是最佳机器学习模型。中胚层特异性转录基因(mesoderm-specific transcript,MEST)、跨膜蛋白98(transmembrane protein 98,TMEM98)、晶状体蛋白λ1(crystallin lambda 1,CRYL1)、染色体14开放阅读框(chromosome 14 open reading frame 147,C14orf147)和跨膜蛋白141(transmembrane protein 141,TMEM141)是XGB 模型中重要性得分最高的前5 位基因(图8-D)。
图8 活动期UC 机器学习模型的构建Fig.8 Construction of active UC machine learning model
选择XGB 模型重要性得分最高的前5 位基因作为疾病特征基因,并构建列线图以预测活动期UC 的发病风险(图9-A)。校准曲线显示实线和虚线在图中彼此很接近,表明模型预测能力较好(图9-B)。决策曲线红线表示的模型远离全部(all)曲线,再次表明模型效果较好(图9-C)。以3 个验证数据集GSE87466、GSE87473 和GSE92415 从外部来验证模型的预测能力,结果发现模型的准确率可以达到91.7%、95.7%和94.9%(图9-D~F)。
图9 活动期UC 机器学习模型的验证Fig.9 Verification of machine learning model of active UC
相关性分析显示,NLRP3与疾病特征基因存在着较强的负相关性,其他DECRGs 中的1个或者多个基因与5 个疾病特征基因有较强的正相关性,它们之间可能存在着相互调控作用,见图10。
图10 DECRGs 与5 个疾病特征基因间的相关性分析Fig.10 Correlation analysis between DECRGs and five disease characteristic genes
共筛选出338 味中药,其中《中国药典》2020年版收录206 味药,根据中医传承计算平台(V3.5)对这206 味药进行出现频次统计、四气五味、归经统计,根据第10 版《中药学》进行功效类别统计。白术、丹参、莪术、紫草、桑叶、青风藤、桑枝均出现在3 个铜死亡基因中,可见其在调节铜死亡治疗活动期UC 中发挥着重要的作用。四气以温、寒和平为主,五味以苦、甘、辛为主,主归肝、脾经,肺经和胃经次之(图11)。中药功效以清热类、补虚类和活血化瘀类为主,清热类中以清热解毒药和清热燥湿最多(表1)。
表1 中药功效分类统计Table 1 Classified statistics of efficacy of traditional Chinese medicines
图11 中药的四气、五味和归经统计Fig.11 Statistics of four qi, five flavor and meridian tropism of traditional Chinese medicines
UC 是一种炎症性肠病,近年来UC 的发病率和患病率逐渐上升,给我国的医疗系统带来了重大的负担和挑战[28-29]。然而,UC 确切的机制尚不清楚,多种因素参与了UC 的发病机制,包括免疫反应紊乱、肠屏障损伤、遗传易感性、肠道生态失调和环境应激[30-31]。UC 不易治愈,复发率高,与结肠癌的发生密切相关,虽然有许多药物可用于抑制结肠炎,但常规治疗具有一定的局限性和严重的不良反应[32]。铜死亡参与多种免疫系统疾病的发生和进展[33-34],然而,目前尚缺乏铜死亡在活动期UC 中的研究。
本研究通过比较CRG 在活动期UC 组和对照组结肠中的表达,发现与正常人群相比,有11 个CRGs 在活动期UC 结肠中的表达异常,活动期UC患者的NLRP3显著上调,其他的DECRGs 均下调,表明DECRGs 在活动期UC 的发展中起重要作用。活动期UC 患者和小鼠模型中均检测到NLRP3 水平升高[35-36]。Zhang 等[37]通过对56 例UC 患者和274 例对照组进行基因分型发现rs10754558 和rs10925019 在中国汉族人群中与UC 的易感性显著相关,这表明NLRP3可能在UC 的发病机制中起重要作用。通过对这些基因进行了相关性分析发现活动期UC 中DECRGs 之间存在着显著的相关性,NLRP3与其他DECRGs 均呈负相关,其中与DBT负相关性最显著。活动期UC 和对照组之间的免疫细胞浸润结果显示,活动期UC 患者中激活的记忆性CD4+T 细胞、静息NK 细胞、M0 巨噬细胞、激活的肥大细胞和中性粒细胞上调,与以往研究结果一致[38-43]。同时,DECRGs 与免疫浸润细胞的相关性分析显示M1 巨噬细胞与DLD、GCSH呈正相关,M0 巨噬细胞与DLAT、FDX1、LIAS呈显著的负相关,这些结果表明DECRGs 可能与活动期UC 进展和免疫浸润有关。
此外,基于DECRGs 使用无监督聚类分析将63例活动期UC 患者分为2 个亚型,以更好地了解CRG 的表达模式。在不同的亚型中DECRGs 的表达存在差异,其中在C2 亚型中NLRP3的表达高于C1 亚型,其他DECRGs 的表达均低于C1 亚型。NLRP3 是核苷酸结合域样受体家族的一员,被激活后与凋亡相关斑点样蛋白(apoptosis-associated speck-like protein containing a CARD,ASC)和效应分子半胱氨酸蛋白酶-1 前体(pro-caspase-1)组成NLRP3 炎性小体[44]。NLRP3 炎症小体是炎症反应系统的重要组成部分和肠道内稳态的重要调节因子,其异常激活后通过释放炎性细胞因子、损伤肠上皮细胞、破坏肠黏膜屏障,诱导参与UC 的发生发展的炎症级联反应[45-46]。NFE2L2 是一种重要的传感器蛋白,可能通过调节促炎细胞因子和诱导II期解毒酶在预防UC 方面发挥潜在的重要作用[47]。FDX1 是铜离子载体诱导细胞死亡的关键调节因子,也是细胞蛋白脂酰化的上游调节因子,通过结合LIAS 促进其与脂酰载体蛋白GCSH 的功能性结合来直接调节蛋白质脂酰化[48]。FDX1 表达与肿瘤细胞静止和炎症呈正相关,与肿瘤侵袭呈负相关,FDX1 的高表达在抑制结肠肿瘤中起着至关重要的作用[49]。PDHA1 是葡萄糖代谢过程中必不可少的成分,其失活被认为与肿瘤中的代谢重编程和乳酸过量产生有关[50]。NLRP3 炎症小体的激活需要乳酸发酵[51]。不同亚型的免疫学分析显示C2 亚型中M0巨噬细胞数量更高。在不同诱导因子的作用下,巨噬细胞可分为M1 和M2 极化型,在UC 患者中巨噬细胞极化表型以M1 型为主,诱导局部炎症反应,加重肠道炎症和黏膜屏障损伤[52]。可见,C2 亚型可能比C1 亚型的预后更差。GSVA 分析表明,C1 亚型在FcεRI 信号通路、趋化因子信号通路、T 细胞受体信号通路等通路中富集。C2 亚型在丁酸代谢、萜类骨架的生物合成、脂肪酸代谢和组氨酸代谢等途径中显著富集。
近年来,WGCNA 算法和机器学习模型越来越多地用于UC 相关基因的筛选和疾病患病率的预测[53-56],但均未涉及到活动期UC。本研究运用WGCNA 筛选出234 个与DECRGs 聚类密切相关的核心基因,通过4 种机器学习发现,XGB 模型是最佳机器学习模型,MEST、TMEM98、CRYL1、C14orf147和TMEM141是活动期UC 的特征基因。然后,构建了列线图模型,使用MEST、TMEM98、CRYL1、C14orf147和TMEM141诊断活动期UC,GSE87466、GSE87473 和GSE92415 数据集用于评价5 个基因的预测效能,ROC 结果表明该模型具有良好的预测性能和一定的临床应用价值,可作为活动期UC 的诊断生物标志物。
通过预测疾病特征基因和DECRGs 的中药,总结发现调控DECRGs 治疗活动期UC 的中药功效以清热类、补虚类、活血化瘀类为主,清热类中以清热解毒药和清热燥湿药最多,体现了活动期UC 以脾胃虚弱为本,湿热毒瘀蕴结为标,与活动期UC常见的大肠湿热、热毒炽盛等证型相符[57-59]。四气以温、寒和平为主,五味以苦、甘、辛为主,甘温健脾能补益脾胃,苦寒燥湿能清脾胃、大肠湿热,辛能散能行,可调气行血以通肠腑。归经方面主归肝、脾、肺、胃经,与疾病的归经基本相同[60],体现了微观铜死亡基因调节活动期UC 层面与宏观整体调节层面的方向相符。频次统计显示白术、丹参、莪术、紫草、桑叶、青风藤、桑枝在铜死亡基因预测中出现的频次最高,可见它们在其中发挥着重要的作用。白术多糖和挥发油可以通过调节肠道菌群的组成和宿主代谢来治疗UC[61-63]。丹参提取物能够通过抑制Toll 样受体4(toll-like receptor 4,TLR4)/磷脂酰肌醇-3-羟激酶(phosphatidylinositol-3-hydroxykinase,PI3K)/蛋白激酶B(protein kinase B,Akt)/哺乳动物雷帕霉素靶蛋白(mammalian target of rapamycin,mTOR)信号通路改善UC 的炎症[64]。莪术提取物中的姜黄素可以通过抑制炎症介质和细胞因子、清除氧自由基等过程发挥免疫调节和抗炎作用治疗UC[65]。紫草素及其衍生物通过抑制NLRP3 炎症小体和NF-κB 信号通路的激活,减轻结肠上皮紧密连接的破坏,从而改善UC 的炎症和黏膜屏障损伤[66]。
综上所述,本研究通过生物信息学揭示了铜死亡相关基因在活动期UC 中的表达模式。通过机器学习,确定了MEST、TMEM98、CRYL1、C14orf147和TMEM141是参与活动期UC 发病机制的核心基因,可作为活动期UC 的诊断生物标志物。本研究还建立了1 个预测模型,可以准确评估活动期UC的患病风险。本研究系统地评估了活动期UC 与铜死亡之间的关系,并预测了能够调节铜死亡治疗活动期UC 的中药,这不仅有助于提高对活动期UC铜死亡的认识,也有利于活动期UC 的诊断和治疗。
利益冲突所有作者均声明不存在利益冲突