《Frontiers in Immunology》：Machine learning meets psoriasis: identifying key lactylation biomarkers as potential targets for diagnosis and therapies
引言
银屑病是一种影响全球约2–3%人口的慢性免疫介导的炎症性皮肤病，其特征是角质形成细胞过度增殖、异常分化以及IL-23/IL-17免疫轴的持续激活。尽管靶向TNF-α、IL-17和IL-23的生物制剂已显著改善了临床疗效，但仍有相当一部分患者反应不完全、停药后复发或随时间推移疗效逐渐丧失。这些临床挑战表明，除了经典的细胞因子信号传导之外，还存在其他调控层级，共同导致了疾病的持续性和异质性。
近年来，免疫代谢领域的进展揭示，银屑病皮损表现出深刻的代谢重编程，呈现“类瓦博格效应”（Warburg-like）表型。角质形成细胞表现出糖酵解增强、乳酸生成增加以及缺氧诱导因子-1α（HIF-1α）信号激活。银屑病皮肤中升高的乳酸不再被视为单纯的代谢副产物，而是一种能够塑造免疫反应和上皮行为的生物活性代谢物。越来越多的证据表明，乳酸有助于维持17型炎症并加强IL-17驱动的转录程序，从而促进慢性疾病进展。
该领域的一个重大概念突破是赖氨酸乳酰化（Kla）的发现，这是一种新型的乳酸衍生翻译后修饰，直接将细胞代谢状态与表观遗传调控联系起来。组蛋白乳酰化已被证明可以调节免疫细胞中的基因表达程序，影响巨噬细胞极化和炎症反应。新出现的证据表明，类似的机制在银屑病中也发挥作用，其中代谢重排和免疫激活紧密相连。近期的机制研究已经开始定义银屑病病理中的特定乳酰化事件。值得注意的是，据报道，IL-17A信号通过KLK8–IL-17受体–HAT1调控轴促进角质形成细胞中的组蛋白乳酰化，建立了一个放大炎症基因转录和角质形成细胞增殖的正反馈回路。乳酸可用性的增加也与抗IL-17A疗法反应性降低有关，而乳酸转运的药理学抑制在实验模型中增强了治疗效果。
除了角质形成细胞，乳酸还充当表皮细胞和浸润免疫细胞群之间的代谢信使。表皮来源的乳酸可以通过单羧酸转运蛋白（MCTs）转运到巨噬细胞和其他免疫细胞中，重塑它们的代谢程序和炎症表型。因此，银屑病可能代表一个代谢协调的炎症生态系统，其中乳酸驱动的表观遗传重塑维持了免疫与上皮的相互作用。
尽管取得了这些进展，但银屑病中全面的乳酰化图谱仍未完全阐明。关键乳酰化蛋白的身份、组蛋白与非组蛋白乳酰化的相对贡献以及乳酰化特征的个体间异质性仍知之甚少。传统的实验方法虽然是必不可少的，但可能无法完全捕捉免疫-代谢调控网络的多维复杂性。
在此背景下，机器学习（ML）提供了一个变革性的机会，可以加速生物标志物的发现和机制推断。通过整合转录组学、表观基因组学、蛋白质组学和临床数据集，机器学习算法可以识别潜在模式，优先考虑候选的乳酰化相关调节因子，并揭示炎症回路中的非线性相互作用。本研究利用生物信息学和机器学习，识别了银屑病中与乳酰化相关的潜在生物标志物。我们分析了来自GEO数据库的银屑病相关数据集，通过差异表达分析和加权基因共表达网络分析（WGCNA）识别关键基因。然后将这些基因与先前研究中发现的乳酰化相关基因（LRG）进行比较，以发现乳酰化相关关键基因（其编码的蛋白质可能发生乳酰化）。使用两种机器学习算法，我们确定了关键的诊断生物标志物，并开发了列线图模型以提高预测准确性和临床相关性。此外，我们进行了双样本孟德尔随机化（TSMR）和基于汇总数据的孟德尔随机化（SMR）分析，以探索这些生物标志物与银屑病之间的因果关系。
材料与方法
研究流程如图1所示。分析使用了来自基因表达综合库（GEO）的三个独立的银屑病基因表达数据集：GSE13355、GSE14905和GSE121212。仅包括健康对照和皮损皮肤样本，排除非皮损样本。GSE13355（64个正常样本，58个银屑病样本）和GSE14905（21个正常样本，28个银屑病样本）是基于GPL570平台的微阵列数据，而GSE121212（38个正常样本，28个银屑病样本）是基于GPL16791平台的RNA测序数据。GSE13355作为训练集，GSE14905和GSE121212作为验证集。单细胞RNA测序数据来自GSE151177（GPL18573），包括五个正常和五个皮损样本。共327个乳酰化相关基因来源于先前文献。
在GSE13355数据集中，使用“limma” R包识别银屑病中表达显著改变的基因，阈值设为|logFC| > 0.5且p值 < 0.05。使用“ggplot2”创建火山图，并对|logFC| > 3且p值 < 0.005的基因进行标记。使用“pheatmap”生成前50个显著差异表达基因的层次聚类热图。通过“clusterProfiler”进行DEGs的GO和KEGG富集分析功能注释，并使用“enrichplot”可视化结果。
为了识别参与银屑病发病机制的关键基因模块，使用R中的WGCNA包进行加权基因共表达网络分析（WGCNA）。使用“cutreeStatic”函数去除离群样本，设定“cutHeight”为76，“minSize”为10。“pickSoftThreshold”函数确定网络构建的软阈值。使用此阈值，“blockwiseModules”函数识别模块，设置“minModuleSize”为60，“mergeCutHeight”为0.25。使用“plotDendroAndColors”可视化层次聚类树状图。皮尔逊相关分析评估基因模块与临床特征之间的关系，通过热图显示。“基因显著性”（GS）评估基因与性状的关联，而“模块成员资格”（MM）描述连接模式。
我们实施了两种机器学习算法，通过交集分析来识别共识基因。使用“randomForest” R包进行随机森林分析，以处理高维数据，采用基尼指数减少来确定变量重要性，并选择前5个基因进行进一步分析。使用“glmnet”包进行LASSO回归分析，采用L1惩罚逻辑回归，并通过10折交叉验证选择正则化参数lambda。该方法通过将不具信息性的变量系数收缩为零来自动选择特征。两种方法共同选中的基因被确定为共识关键基因。然后，我们使用“rms” R包基于这些关键诊断基因构建了临床列线图。
使用R中的“pROC”包通过ROC曲线分析评估候选基因的诊断性能，AUC值表示每个基因对疾病的分类能力。使用“ggpubr”包的箱线图可视化基因表达模式，以便跨样本组比较转录谱。
从长春益思实验动物技术有限公司获得雄性BALB/c小鼠，饲养在受控环境（22 ± 2 °C，45-55%湿度）中，具有12小时光/暗周期，自由获取食物和水。开始处理前，无菌剃除一块2.5×3 cm的背部皮肤区域。实验涉及两组（每组n=5）。通过每日涂抹62.5 mg的5%咪喹莫特乳膏连续七天诱导银屑病模型，对照组不做处理。每日拍摄疾病进展照片，并根据红斑、鳞屑和皮肤增厚评估皮损严重程度，使用改良的银屑病面积和严重程度指数（PASI），每项评分从0（无）到4（非常严重）。这些分数的总和构成PASI评分。
石蜡包埋的组织切片脱蜡、复水，苏木精染色1分钟，然后冲洗。用1%盐酸酒精处理20秒，冲洗，再用1%氨水处理30秒，再次冲洗。切片在85%和95%乙醇中各脱水5分钟，伊红染色30秒，用无水乙醇冲洗。用无水乙醇重复脱水5分钟，接着用二甲苯处理并用中性树胶封片。然后在光学显微镜下检查病理变化。
使用TRIzol从组织中提取RNA，并通过Prime-Script RTase进行逆转录合成cDNA。然后使用cDNA和特异性引物（详见Supplementary Table 1）进行实时定量PCR（RT-qPCR）。使用2-ΔΔCt方法计算相对mRNA表达量，以GAPDH作为内参基因。
我们使用CIBERSORT方法量化组织样本中的免疫细胞，仅关注具有统计学显著性的标本（P值 < 0.05）。小提琴图说明了患病组和健康组之间免疫细胞浸润的差异。我们还使用皮尔逊分析探讨了关键诊断生物标志物与免疫细胞之间的相关性，结果以热图显示。
通过NetworkAnalyst构建了miRNA-基因和TF-基因调控网络。靶向基因的microRNA来自TarBase，转录因子来自JASPAR。使用Enrichr上的DSigDB（
https://maayanlab.cloud/Enrichr/）预测针对关键生物标志物的潜在药物。
使用“Seurat” R包处理单细胞测序数据。质量控制保留了满足以下条件的细胞：（1）200–6000个基因，（2）检测到的计数 > 500，（3）线粒体基因 ≤ 10%，（4）血红蛋白基因 ≤ 1%。使用“NormalizeData”对基因表达数据进行归一化，使用“vst”方法识别前2000个可变基因。使用“RunPCA”进行降维之前，使用“ScaleData”进行数据缩放。使用“RunHarmony”去除批次效应。使用“FindNeighbors”和“FindClusters”进行细胞聚类。初始细胞注释使用“singleR” R包，随后使用CellMarker数据库进行手动校正以确保准确性。
我们从IEU OpenGWAS数据库获得了银屑病GWAS数据集（GWAS ID: ukb-b-10537），包括5314个病例和457619个对照。相关生物标志物的eQTL数据来自eQTLGen数据库。使用TwoSampleMR R包进行双样本孟德尔随机化（MR）分析。我们选择与基因强相关的SNP（p < 5.0×10-8），并去除处于连锁不平衡（R2 < 0.001，10000 kb窗口）中的SNP。排除F统计量低于10的SNP以避免弱工具变量偏倚。我们的主要分析使用逆方差加权（IVW）估计来探索基因表达与银屑病风险之间的因果关系，并通过MR Egger、加权中位数、简单模式和加权模式分析进行进一步验证。敏感性分析包括用于异质性的Cochrane‘s Q检验和用于水平多效性的MR-Egger截距检验。
使用SMR软件进行分析，该分析采用top cis-eQTL SNV作为工具变量，结合eQTL和GWAS数据来研究基因表达与疾病之间的因果关系。HEIDI检验有助于区分因果关联与连锁不平衡。当pSMR< 0.05且pHEIDI> 0.05时，表明基因表达与银屑病之间存在因果关系。
使用R（版本4.4.3）进行统计分析。使用学生t检验评估正态分布连续数据的差异，而对非正态数据使用Wilcoxon秩和检验。报告双尾p值，显著性水平设为p < 0.05。使用皮尔逊相关系数评估变量关联。
结果
差异表达基因的鉴定与功能分析
火山图显示在GSE13355数据集中有1483个上调基因和1772个下调基因。层次聚类突出了前50个显著基因之间不同的表达模式。功能注释将上调基因与细胞周期和免疫调控联系起来，GO分析显示在“染色体分离”和“趋化因子活性”方面富集，KEGG分析将其与“细胞周期”通路相关联。下调基因在GO术语如“含胶原的细胞外基质”中富集，并与KEGG通路如“黏着斑”和“Rap1信号通路”相关联。
识别银屑病相关基因模块并阐明其功能
我们对GSE13355进行了WGCNA分析，去除了一个离群样本（cutHeight为76）。选择软阈值为6以确保构建无标度网络。我们识别出22个基因共表达模块，其中与银屑病显著相关。青绿色模块与疾病状态呈最高正相关（r = 0.96, p = 1e-66），而绿色模块呈最强负相关（r = -0.78, p = 6e-26）。这些发现表明这些基因簇可能在银屑病发展中发挥功能性作用。青绿色模块（0.98）和绿色模块（0.85）中的“基因显著性”（GS）和“模块成员资格”（MM）之间存在强相关性，凸显了它们在银屑病发病机制中的重要性。
这些模块的GO和KEGG分析显示，青绿色模块与细胞增殖和免疫反应调节相关。GO分析显示在诸如“有丝分裂染色体分离”、“核分裂”和“趋化因子介导的信号传导”等术语中富集，而KEGG通路分析则突出了“细胞周期调节”和“IL-17介导的炎症反应”。相反，绿色模块在与细胞-基质相互作用和细胞外组织相关的GO类别中富集，以及与KEGG通路如“Rap1依赖性信号级联”相关联。值得注意的是，青绿色模块的基因功能与上调基因一致，而绿色模块的基因功能与下调基因一致。分析结果的一致性增强了我们研究的可信度，并为银屑病发病机制提供了新的见解。
乳酰化相关关键基因的鉴定及列线图模型的构建
我们观察到青绿色模块中的基因与上调基因之间，以及绿色模块中的基因与下调基因之间的功能相似性。为了识别银屑病中的关键基因，我们将青绿色模块基因与上调基因、绿色模块基因与下调基因取交集，识别出1623个关键基因。这些基因再与文献中的乳酰化相关基因取交集，最终得到26个乳酰化相关关键基因（LRKG），其编码的蛋白质可能发生乳酰化。随后，我们使用随机森林和LASSO回归从这些LRKG中筛选关键诊断基因。随机森林算法根据重要性确定了前5个基因，而LASSO回归找到了8个非零系数的基因。四个基因——MKI67、MPHOSPH6、ENO1和FABP5——是两种方法共有的。为了评估它们的综合诊断能力，使用这四个生物标志物创建了多变量逻辑回归模型，如列线图所示。
使用表达分析和ROC曲线验证诊断效能
为了评估四个生物标志物的诊断潜力，我们在训练集和验证集中检查了它们的表达。箱线图显示，在银屑病皮损中MKI67、MPHOSPH6、ENO1和FABP5的水平显著高于健康对照组。ROC曲线分析证实，每个基因以及组合模型的AUC值均超过0.80，表明具有强大的诊断性能。包含所有四个基因的逻辑回归模型在一个数据集中达到了0.999的AUC值，在另外两个数据集中达到了1.000的AUC值，强调了这四基因特征对于银屑病的稳健诊断能力。
在小鼠银屑病模型中验证关键基因表达
为了验证乳酰化相关生物标志物，我们创建了咪喹莫特（IMQ）诱导的银屑病小鼠模型。IMQ处理的小鼠表现出明显的银屑病皮炎特征，并且与对照组相比，具有更高的红斑、鳞屑、厚度和PASI评分。H&E染色证实模型诱导成功，IMQ处理的小鼠出现明显的表皮增生。RT-qPCR显示，与对照组相比，IMQ处理小鼠的皮损皮肤中MPHOSPH6、ENO1、FABP5和MKI67的表达显著增加。
免疫浸润分析
为了描述银屑病皮肤与健康对照相比的免疫细胞浸润特征，我们使用CIBERSORT算法分析了GSE13355的组织样本。小提琴图显示，与对照组相比，银屑病中活化的CD4+记忆T细胞、滤泡辅助T细胞、调节性T细胞、活化的NK细胞、单核细胞、M0巨噬细胞、M1巨噬细胞和活化的树突状细胞的水平更高。相比之下，记忆B细胞、静息树突状细胞和静息肥大细胞在银屑病皮损中显著降低。此外，关键诊断基因大多与幼稚B细胞和滤泡辅助T细胞等免疫细胞呈正相关，而与静息肥大细胞和记忆B细胞呈负相关。这些发现突出了连接基因表达和免疫细胞活性的调控网络，暗示了控制炎症的潜在治疗靶点。
miRNA-基因和TF-基因调控网络及银屑病生物标志物的药物靶点预测
微小RNA（miRNA）和转录因子（TF）在调节基因表达中至关重要。为了理解关键诊断基因的调控机制，我们使用计算预测构建了miRNA-基因和TF-基因网络。miRNA-基因网络包括381个节点（377个miRNA和4个诊断基因）和595个相互作用，其中MKI67受281个miRNA调控，ENO1受247个调控，MPHOSPH6受49个调控，FABP5受18个调控。值得注意的是，有9个miRNA调控所有四个基因。TF-基因网络有33个节点（29个TF和4个基因）和39个相互作用，其中ENO1受14个TF调控，FABP5受12个，MKI67受7个，MPHOSPH6受6个。FOXC1和NFIC调控ENO1、FABP5和MKI67，凸显了它们作为潜在关键转录因子进行进一步研究的重要性。我们在DSigDB数据库中识别出103种药物（校正后p值 < 0.05）靶向四个生物标志物，统计上最显著的前五种药物已列出。
在单细胞测序数据中验证生物标志物的表达谱和细胞分布
聚类分析识别出17个不同的细胞群体。经过SingleR注释和手动校正后，识别出五种主要细胞类型：“T细胞”、“角质形成细胞”、“树突状细胞”、“巨噬细胞”和“黑色素细胞”。比较分析显示，与健康对照相比，皮损皮肤中T淋巴细胞和树突状细胞显著增加，突出了它们在银屑病中的作用。这与CIBERSORT算法的结果一致。点图显示MKI67在所有细胞簇中表达较低，而其他基因在细胞类型间的表达存在差异。UMAP可视化进一步详细展示了诊断标志物在不同细胞簇中的表达谱和分布。
双样本孟德尔随机化分析揭示MPHOSPH6和ENO1是银屑病的风险生物标志物
为了检查关键诊断基因是否与银屑病存在因果关系，我们使用MPHOSPH6、ENO1、MKI67和FABP5的eQTL数据作为暴露变量，银屑病GWAS数据作为结局变量，进行了双样本孟德尔随机化（TSMR）分析。主要使用的方法是IVW。我们的发现表明MPHOSPH6、ENO1与银屑病之间存在显著的因果关系，其比值比（OR）大于1，表明它们的过表达可能增加银屑病的风险。图表显示了五种方法评估MPHOSPH6与银屑病因果关系的OR和p值，以及SNP位点效应的森林图。散点图显示了每个SNP位点对银屑病和基因表达的影响，五种方法一致地将MPHOSPH6与增加的银屑病风险联系起来。漏斗图显示了对称的SNP分布，表明无异质性。留一法敏感性分析证实了我们研究结果的稳定性，因为移除任何SNP都不会显著改变整体效应估计，从而加强了我们的结论。补充材料详细列出了ENO1的结果。Cochran‘s Q检验显示没有显著的异质性（p > 0.05），MR-Egger回归截距接近于零且p值大于0.05，表明没有水平多效性。总之，我们的TSMR分析有力地表明，过表达MPHOSPH6和ENO1会增加银屑病的风险。
SMR分析进一步表明MPHOSPH6是银屑病的风险生物标志物
我们进行了SMR分析，以bSMR为4.280e-03、pSMR为2.660e-03确认了MPHOSPH6与银屑病之间的因果关系，两者均具有统计学显著性。Top SNP是rs34638657，pHEIDI值为0.061表明没有显著的异质性。这些结果与TSMR分析一致，支持较高的MPHOSPH6表达与增加的银屑病风险之间存在正向因果关系，从而加强了我们的结论。
讨论
银屑病是一种慢性炎症性皮肤病，其特征是免疫反应失调和表皮增殖加速。了解银屑病背后的分子机制对于识别生物标志物和潜在治疗靶点至关重要。在本研究中，我们对GSE13355数据集的基因表达谱进行了全面分析，以探索参与银屑病发病机制的关键基因。差异表达分析揭示了银屑病皮损中存在大量上调基因和下调基因。上调基因的功能注释强调了它们在诸如“染色体分离”和“趋化因子活性”等生物学过程中的富集，突出了异常细胞分裂和免疫细胞激活在银屑病发展中的重要性。下调基因在与“含胶原的细胞外基质”相关的过程富集，表明组织稳态和细胞外基质重塑的破坏，这是银屑病皮损的特征。这些发现提供了银屑病的详细分子图谱，强调了免疫失调和异常细胞周期
