公共转录组的再分析揭示了重度抑郁症与皮肌炎在单细胞背景下共享的免疫特征

摘要

本研究旨在结合公共GEO数据集的综合生物信息学再分析，结合单细胞情境化，识别重度抑郁症与皮肌炎之间的候选共享基因，并表征其在与皮肤肌炎相关的单细胞数据集中各免疫细胞子集的分布。

摘要

本研究旨在通过对公共GEO数据集的整合生物信息学再分析，识别重度抑郁症与皮肤肌炎之间的候选共享转录组信号。分析工作流程包括加权基因共表达网络分析（WGCNA）用于关键模块识别，基因本体（GO）和京都基因组百科全书（KEGG）富集分析用于功能表征，GeneMANIA及基于网络可视化平台的网络分析用于候选基因优先排序，以及结合SHapley加法解释（SHAP）进行诊断特征选择的113个机器学习模型的评估。随后进行了基因集富集分析（GSEA）、免疫浸润分析以及基于单细胞RNA测序的情境分析，以进一步表征已识别信号的免疫相关细胞环境。整合了与皮肤肌炎相关的GEO数据集，识别出570个差异表达基因，其中33个候选共享基因通过WGCNA获得了。功能富集和网络分析强调了免疫防御、细胞毒性及PPAR、IL-17和抗原加工等通路，ELANE、PPBP和CTSG作为高度连接的节点出现。基于机器学习的特征优先级保留了8个候选基因，分别是KIF4A、OLR1、KIR2DL4、KRT23、KIR3DS1、AZU1、SCG5和LRRC37E。免疫浸润分析将这些共享基因与调控T细胞（Tregs）、静息肥大细胞、静息树突状细胞以及经典活化（M1）和替代活化（M2）巨噬细胞联系起来。单细胞RNA测序情境分析进一步表明，具有不同候选基因评分状态的CD8⁺ T细胞亚组表现出不同的细胞间通信模式。其中，MIF–（CD74+CD44）轴和来自初真/中央存储T细胞的信号是需要进一步验证的显著特征。总体而言，本研究鉴定了重度抑郁障碍与皮肌炎之间候选共享的转录组信号，并强调了需要在真实共病队列中进一步验证的免疫相关细胞环境。

引言

皮肌炎是一种慢性系统性自身免疫性疾病，其特征是皮肤和骨骼肌的炎症受累，临床表现为对称的近端肌无力和独特的皮肤病变，严重情况下还会出现多^{器官功能障碍}1。临床证据的积累表明，皮肌炎患者常伴有精神疾病，最显著的是重度抑郁障碍^2,3,4。皮肌炎相关的重度抑郁障碍的发病机制是多因素的，源于心理社会困扰、神经内分泌失调和系统性免疫炎症之间的复杂相互作用。持续的疼痛、疲劳和逐渐的肌肉无力会严重损害身体功能和生活质量。这些负担可能导致社会角色的中断、长期的心理压力和自主权的^降低。此外，长期的糖皮质激素暴露可能扰乱下丘脑-垂体-肾上腺（HPA）轴，损害海马体可塑性，从而增加对抑郁相关表现的易感^性6。此外，皮肤肌炎中持续的免疫激活和全身性炎症越来越多地被认为是抑郁^{相关症状的}潜在促成因素7。这些外周介质可能通过调节神经递质代谢和神经可塑性来影响中枢神经系统，从而将系统性自身免疫与神经精神表现联系起来 7,8,9,10。

重要的是，这一理由并不意味着这两种疾病在生物学上都是同质性的。皮肌炎包括临床和血清学上不同的亚组，包括抗MDA5和抗TIF1-γ相关表型，具有不同的炎症和临床特征11,12,13。重度抑郁障碍同样越来越被认可为一种异质性疾病，现有证据支持存在一种免疫炎症亚型，而非单一的通用炎症特征（^14,15）。因此，本研究并非设计为一刀切的共享分子程序，而是筛查可在独立公开数据集中队列层面检测到的候选重叠免疫相关转录组信号。

除了心理社会压力和治疗暴露外，重度抑郁症与皮肌炎之间更可检验的生物学联系是共享的免疫炎症失调。重度抑郁症是一种异质性疾病，不应假设它具有单一的通用转录组特征。然而，越来越多的证据支持与炎症相关的重度抑郁障碍亚型，外围转录组研究已发现患者中先天免疫、中性粒细胞相关、干扰素和补体通路的失调（14,16,17）。与此同时，MAPK相关的应激信号也被认为与抑郁表型有关¹⁸。相比之下，皮肌炎是一种广为人知的干扰素驱动自身免疫疾病，血液及受影响组织中的转录组研究持续显示I型干扰素及更广泛的免疫炎症程序被激活;近期多组学分析进一步强调了ERK和p38 MAPK相关通路在皮肤肌炎19,20,21中的活性。综合来看，这些发现为探讨免疫相关转录组信号子集是否在独立公开数据集中重度抑郁症与皮肌炎之间存在重叠提供了生物学上的合理依据。

尽管有这些观察，重度抑郁障碍与皮肌炎重叠的分子基础仍不够明确。重要的是，目前公开的数据集并未提供同时被诊断为重度抑郁障碍和皮肌炎患者的真实队列。因此，本研究旨在通过整合生物信息学再分析框架，识别不同公开数据集中重度抑郁症患者抑郁症的抑郁症，而非抑郁症患者的抑郁症，而非抑郁症患者的候选共享转录组^信号。具体而言，通过差异表达分析、加权基因共表达网络分析（WGCNA）、功能富集分析、基于网络的分析以及基于机器学习的特征优先排序，分析了公开可用的转录组数据集，以识别跨疾病候选基因和通路²³。此外，分析了与皮肤肌炎相关的单细胞数据集，以在免疫细胞层面为这些候选基因提供背景。如 图1所示，整体分析工作流程以逐步流程图总结。本研究旨在建立一个以假设为导向的框架，以识别重度抑郁障碍与皮肤肌炎之间的候选共享分子信号。

因此，本研究采用了逐步优先排序框架。疾病相关共表达模块最初分别在重度抑郁障碍和皮肤肌炎数据集中被单独识别，并利用其重叠部分定义候选的跨疾病共享信号。这些候选者随后通过富集和基于GeneMANIA的网络分析进行功能性情境化，在以皮肤肌炎为中心的机器学习分类框架内优先排序，最终在与皮肤肌炎相关的单细胞数据集中进行分析，以提供细胞情境化。

图1：数据收集和分析流程图。请点击此处查看该图的更大版本。

已有多种替代方法用于研究跨病菌分子重叠。差异表达基因（DEG）列表的简单交叉计算简单，但缺乏共表达分析提供的模块级上下文信息，且对任意折叠变化和P值阈值敏感。传统荟萃分析会在同一疾病研究间进行效应值的综合分析，但并不设计用来识别两种不同疾病的共享信号。当前工作流程整合了多个互补的分析层——共表达模块重叠、功能丰富、网络分析、基于机器学习的特征优先级、免疫反卷积和单细胞情境化——每个层在顺序优先级框架内各自承担着不同的任务。这种多层设计有助于逐步减少候选基因数量，并在多个层面提供交叉验证的生物学背景。该方案适用于任何存在公开可得的大规模转录组数据集（可选单细胞数据集）的疾病组合，尤其是在无法获得真正共病队列的情况下。然而，该工作流程是观察性的，不包含正式的因果推理框架;所有发现应被解释为假设生成，并需要独立的实验验证。

总体而言，分析工作流程被设计为顺序优先级策略，而非直接的因果推断框架。每个步骤都有不同的目的：基于WGCNA的模块重叠用于候选共享信号识别，富集/网络分析用于生物情境化，机器学习用于皮肤肌炎相关分类中的特征优先级，以及单细胞分析用于细胞类型级情境化。

方案

本研究仅使用了来自基因表达综合数据库（GEO）的公开可得去标识化数据集。由于该工作涉及对现有公共数据的二次分析，且不包括直接参与者联系、干预或获取可识别个人信息，因此无需额外的伦理委员会批准和知情同意。

数据源与预处理

所有基因表达和单细胞数据集均来自GEO数据库²⁴。对于重度抑郁障碍，使用了数据集GSE98793，包括128名患者和64名健康对照组的外周血样。针对皮肌炎，数据集的选择基于预设标准，包括智人表达分析、可明确识别的疾病和控制组、可用于探针到基因定位的平台注释，以及适合发现或验证分析。当GEO系列包含多种炎症性肌病亚型时，本研究仅提取了皮肌炎和正常对照样本。GSE1551、GSE46239和GSE128470被用作发现/训练数据集，而GSE5370、GSE39454和GSE11971则作为独立的验证数据集。本研究分析的皮肤肌炎数据集主要来自受影响的肌肉或皮肤组织，而非外周血液。皮肌炎的单细胞数据来源于GSE190510数据集。

原始表达矩阵及相应的平台注释文件均已从GEO数据库下载。探针ID根据制造商提供的GPL注释映射到官方基因符号。无法明确定位到单一官方基因符号的探针被移除。当多个探针映射到同一基因时，它们会在基因层面使用Limma包中的“avreps”函数实现的平均表达值进行合并，从而生成基因按样本表达矩阵。

为了减少强度依赖偏差并稳定方差，根据表达式值分布在适当情况下应用了log2变换。数组间归一化随后通过limma包中的“normalizeBetweenArrays”函数完成。缺失值（如存在）通过K最近邻补补来推算。对于集成的皮肤肌炎训练数据集，批量校正使用SVA包中的“ComBat”功能完成，数据集/平台来源作为批次变量，样本组（皮肌炎与健康对照组）纳入设计矩阵，以保留批次调整过程中关注的生物学变异。

所有分析均通过桌面操作系统上的 R 集成开发环境在 R 中完成。Limma封装用于探针的摘要和归一化。SVA包用于ComBat批次校正。缺失值通过k=10的K最近邻补值进行补值。

加权基因共表达网络分析

加权基因共表达网络分析（WGCNA）分别针对重度抑郁症和皮肌炎数据集，使用WGCNA R^25,26进行。样本通过 flashClust 进行层级聚类以识别异常值;排除树轮高度超过100且方差最低25%的基因样本。对于每个网络，使用pickSoftThreshold选择软阈值功率（β），以实现近似无尺度拓扑（R² > 0.8）。邻接矩阵被转换为拓扑重叠矩阵（TOM），并通过动态树切割识别模块，最小模块大小为60，合并切割高度为0.25²⁷。WGCNA R 包与 flashClust 一起使用用于层级聚类。为了可重复性，随机种子设置为12345。模块特征基因与疾病状态相关性使用皮尔逊相关系，P值通过Benjamini–Hochberg方法进行调整。对于每种疾病，显示与疾病状态关联最强且最显著的模块被保留为关键疾病相关模块。重度抑郁症数据集中的关键模块基因与皮肌炎数据集中的关键模块基因重叠被定义为下游分析的候选共享基因集。对整合皮肤肌炎队列进行了差异表达分析，以表征与皮肤肌炎相关的转录变化。

功能富集分析

基因本体（GO）富集分析使用R进行。基因符号通过org转换为Entrez ID。Hs.eg.db，以及显著富集的GO项（p < 0.05）均通过enrichGO在clusterProfiler中识别。为了多维可视化结果，使用enrichplot软件生成条形图和气泡图，同时用circlelize软件构建圆图，显示基因类别、基因计数和富集因子。图例是通过ComplexHeatmap软件包添加的。京都基因与基因组百科全书（KEGG）对差异表达基因的途径富集分析也在R中进行了。基因符号基于该组织转换为Entrez IDs。Hs.eg.db数据库及显著富集通路（FDR < 0.05）通过clusterProfiler包28,29,30,31中的enrichKEGG函数识别。富集结果通过条形图和气泡图可视化。

基于GeneMANIA的功能关联网络分析

基于先前识别的共享基因，构建了一个基于GeneMANIA的功能关联网络，以探索这些基因及其相关伙伴之间的相互作用环境。基因列表提交给GeneMANIA，参考物种为智人。GeneMANIA整合了多种证据类型，包括共表达、物理相互作用、途径、共定位、遗传相互作用和共享的蛋白质结构域。最终生成的网络被导出并导入到网络可视化平台中，用于可视化和分析。随后，在网络可视化平台上对网络进行了拓扑分析，用于可视化和分析，以识别高度连通的候选节点^32、33、34。

基于机器学习的诊断模型构建

诊断分类使用了多种机器学习算法，包括随机森林（RF）、支持向量机（SVM）、线性判别分析（LDA）、朴素贝叶、梯度增强机（GBM）、XGBoost、glmBoost、弹性网（Enet）、脊线、最小绝对收缩与选择算子（LASSO）、逐步广义线性模型（Stepglm）和偏最小二乘回归广义线性模型（plsRglm）³⁵.采用了两阶段建模框架，生成了113个候选模型组合。第一阶段，初始算法用于训练队列的变量筛选;第二阶段，保留变量用于拟合诊断分类模型。包含≤5个选定变量的模型被排除在进一步比较之外。合并后的皮肌炎数据集作为训练队列，标签定义为皮肤肌炎与健康对照组，独立验证队列则用于外部性能评估。内部重采样和调优是针对特定算法的：基于glmnet的模型（如LASSO、Ridge和Elastic Net）采用10折交叉验证来选择lambda.min;GBM采用10折内部交叉验证来确定树的最优数量;XGBoost采用5重采样法根据最小测试日志损耗选择最终加振轮;glmBoost 采用基于 cvrisk 的内部交叉验证来确定停止迭代;LDA则采用了筛槽交叉验证框架。对于当前实现中没有显式调优步骤的算法，采用了固定或包默认设置。为减少信息泄露，特征选择、模型拟合和内部调优仅使用训练队列进行，而验证队列仅用于独立预测和基于AUC的性能评估。Caret包用于机器学习工作流管理，针对单个算法支持glmnet、randomForest、e1071、gbm、xgboost、mboost、plsRglm和MASS。SHAP分析使用shapviz软件包进行。随机种子在每次模型配对前设为12345。选定特征少于5个的模型被排除在外。通过SHapley加法解释（SHAP）进一步评估模型可解释性和基因层面贡献，并优先将最具信息量的基因作为候选模型选择特征进行下游生物学解释。

诊断性能评估

使用“pROC”R软件包生成受试者工作特征（ROC）曲线，以评估候选生物标志物的诊断性能。候选标记的表达水平和预测准确性在独立数据集（GSE5370、GSE11971和GSE39454）中得到了验证。模型表现进一步使用混淆矩阵进行评估。通过火山图和箱形图可视化关键模块基因的差异表达，并构建了ROC曲线以评估单个基因的诊断价值。

基因集富集分析

为探索与候选共享转录组信号相关的协调功能变化，使用clusterProfiler^36,37进行了基因集富集分析（GSEA）。皮肤肌炎和对照样本的基因表达数据按表达差异进行排序。使用对应KEGG通路的预定义基因集（c2.cp.kegg.Hs.symbols.gmt）来评估各通路内的基因是否表现出协调的上调或下调趋势。统计显著性定义为P < 0.05。

免疫细胞浸润分析

归一化、log2转制和批次校正的皮肤肌炎基质被用于免疫解卷。CIBERSORT算法被应用到LM22参考矩阵³⁸中估计免疫细胞亚型的相对丰度。去除卷积 P <0.05的样本被保留用于后续分析。通过箱形图可视化了各组间推断免疫细胞比例的差异，并进行了Spearman相关分析以评估免疫细胞亚组与候选共享基因之间的关联。

单细胞RNA测序分析用于细胞情境化

使用修拉在R中进行了单细胞RNA测序分析。Harmony用于批量校正，DoubletFinder用于双重检测，celda/decontX用于环境RNA估计，Monocle用于伪时间轨迹分析，CellChat用于细胞间通信分析，AUCell用于基因组活性评分，GSVA用于ssGSEA评分。原始计数矩阵导入Seurat对象，参数最小单元=5，最小特征=300。计算了每个细胞的质量控制指标，包括线粒体、核糖体和血红蛋白基因比例。只有当细胞满足以下所有条件时才被保留：nFeature_RNA >500、nCount_RNA <5,000、percent_mito <25、percent_ribo >3和percent_hb <1。排除检测到少于3个细胞的基因。此外，在下游分析前还去除了MALAT1和线粒体基因。初步过滤后，使用DoubletFinder在每个样本中识别出双重态，PCs=1：30，pN=0.25;预期双重蛋白比例根据样本特异性细胞数设定（<4,000个细胞：2.5%;4,000–8,000个细胞：5%;>8,000个细胞：6.5%）。仅保留单人单发。利用decontX进一步估算环境RNA污染情况，污染评分<0.2的细胞被保留。

过滤数据通过LogNormalize方法进行归一化，比例因子为10,000，随后进行变量基因的识别、数据标度和主成分分析。批次效应通过 Harmony 校正，批次变量为 orig.ident。前15个Harmony维度用于UMAP可视化和邻居图构建。聚类使用FindNeighbors和FindClusters进行，最终聚类结果定义分辨率为0.05。细胞类型根据典型标记基因及FindAllMarkers结果进行人工注释，³⁹。

在下游功能情境化中，候选基因活性在单细胞层面进行评估，相关免疫细胞亚集则接受轨迹分析和细胞间通信分析。采用基于DDRTree的降维单片镜进行伪时间分析，随后进行细胞排序。使用CellChat进行细胞间通信分析，利用人类配体-受体数据库（限于分泌信号类别）进行，且过滤掉少于10个细胞的通信。

对于每个细胞，候选基因活性通过三种互补方法量化：AUCell、ssGSEA和AddModuleScore。AUCell得分基于基因排序矩阵计算，ssGSEA分数则使用GSVA框架生成。AddModuleScore是通过Seurat内置函数计算的。所得的AUCell、ssGSEA和AddModuleScore值随后合并为单一评分矩阵。每种评分类型首先通过Z分数变换标准化，随后通过最小-最大归一化重新标度至0–1范围。每个单元格的最终综合得分（“评分”）定义为三个归一化分数的总和：

评分 = 归一化 AUCell + 归一化 ssGSEA + 归一化 AddModuleScore。

在下游亚组分析中，提取CD8⁺ T细胞亚组，并根据该子集内的中位评分值进行二分。得分值高于中位数的单元格被分配到High_Hub_genes组，其余单元格则分配给Low_Hub_genes组。

结果

重度抑郁障碍与皮肌炎之间共享候选基因的鉴定

经过对与皮肤肌炎相关的GEO数据集的数据合并、归一化和批量校正（图2A，B），共识别出570个差异表达基因（图2C，D），其中包括517个上调基因和53个下调基因。该皮肌炎差异表达分析用于表征疾病相关的转录变化。

与此同时，WGCNA被单独应用于重度抑郁症和皮肌炎数据集，以识别与疾病相关的共表达模块。重度抑郁障碍数据集的最佳软阈值为4，皮肌炎数据集为9（见图3A，B）。每个数据集检测到五个模块。模块-性状相关分析显示，灰色模块与重度抑郁障碍的关联最为强烈（Cor = −0.28，p = ...

讨论

皮肌炎是一种慢性系统性自身免疫性疾病，皮肤和肌肉受累显著，临床观察显示，皮肌炎患者可能还会出现显著的精神负担，包括与重度抑郁障碍相符的症状。在此背景下，本研究采用整合生物信息学再分析框架，识别重度抑郁障碍与皮肌炎之间的候选共享转录组信号，并结合了额外的免疫浸润分析和单细胞情境化。整合了WGCNA、功能富集分析、基于GeneMANIA的网络分析以及基于机器学习的特征优先排序，以筛选候选的跨病基因和通路。在此框架下，重点突出了包括KIF4A、OLR1、KIR2DL4、KRT23、KIR3DS1、AZU1、SCG5和LRRC37E在内的八个模型选择基因，供后续解释。采用免疫浸润分析来考察这些基因观察到的免疫相关环境，并在与皮肤肌炎相关的数据集中进行单细胞分析，特别是针对CD8⁺ T细胞中的KIR2DL4和AZU1，提供了额外的细胞情境信息。

需要注意的是，本研究中识别的跨疾病信号在两种疾病之间并不对称。与重度抑郁障碍相关的模块与疾病状态的相关性仅较小，而关键的皮肤肌炎相关模块则表现出显著更强的相关性。因此，这里识别的重叠...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
AddModuleScore	Seurat函数	版本4.4.0	Seurat中模块得分计算 RRID: NA
AUCell	Bioconductor	版本1.32.0	单细胞基因集活动评分 RRID: SCR_021327
caret	CRAN	版本7.0.1	机器学习工作流支持 RRID: SCR_022524
celda / decontX	Bioconductor	版本1.24.0	环境RNA污染估计 RRID: NA
CellChat	GitHub / CellChat	版本2.2.0	细胞间通信分析 RRID: SCR_021946
CIBERSORT / LM22签名矩阵	CIBERSORT	LM22	免疫细胞浸润估计 RRID: NA
clusterProfiler	Bioconductor	版本4.12.6	功能富集分析 RRID: SCR_016884
Cytoscape	Cytoscape Consortium	版本3.10	网络可视化平台用于可视化和分析 RRID: SCR_003032
DoubletFinder	GitHub / McGinnis Lab	版本2.0.4	单细胞数据集中的双重检测 RRID: NA
e1071	CRAN	版本1.7.16	支持向量机和朴素贝叶斯建模 RRID: NA
gbm	CRAN	版本2.2.2	梯度提升机建模 RRID: NA
基因表达综合数据库（GEO）	国家生物技术信息中心（NCBI）	GSE98793	重度抑郁症大体转录组数据集 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE1551	皮肌炎训练数据集；骨骼肌活检样本 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE46239	皮肌炎训练数据集；皮肤活检样本 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE128470	皮肌炎训练数据集；从炎症性肌病队列中提取的皮肌炎样本 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE5370	独立的皮肌炎验证数据集；未经治疗的成人肌肉样本 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE11971	独立的皮肌炎验证数据集 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE39454	独立的皮肌炎验证数据集；从炎症性肌病队列中提取的皮肌炎样本 RRID: NA
基因表达综合数据库（GEO）	NCBI	GSE190510	与皮肌炎相关的单细胞RNA-seq数据集 RRID: NA
GeneMANIA	多伦多大学 / GeneMANIA	本研究访问的Web服务器版本	功能关联网络构建 RRID: RRID:SCR_005709
glmnet	CRAN	版本4.1.8	LASSO、Ridge和Elastic Net建模 RRID: NA
GSVA	Bioconductor	版本2.0.7	ssGSEA评分 RRID: NA
Harmony	CRAN	版本1.2.4	单细胞数据整合的批次校正 RRID: NA
limma	Bioconductor	版本3.60.6	差异表达分析、探针汇总和归一化实用工具 RRID: SCR_010943
MASS	CRAN	版本7.3.61	线性判别分析 RRID: NA
mboost	CRAN	版本2.9.11	glmBoost建模 RRID: NA
Monocle	Bioconductor	版本2.38.0	伪时间轨迹分析 RRID: SCR_016339
org.Hs.eg.db	Bioconductor	版本3.19.1	人类基因注释数据库 RRID: NA
plsRglm	CRAN	版本1.5.1	偏最小二乘广义线性建模 RRID: NA
pROC	CRAN	版本1.18.5	ROC曲线分析 RRID: SCR_024286
R统计软件	R统计计算基金会	版本4.4.2	主要统计计算环境 RRID: SCR_001905
randomForest	CRAN	版本4.7.1.2	随机森林建模 RRID: SCR_015718
RStudio	Posit Software, PBC	版本2024.4.1.748	R的集成开发环境 RRID: SCR_000432
Seurat	CRAN / Satija实验室	版本4.4.0	单细胞RNA-seq预处理、聚类和可视化 RRID: SCR_016341
shapviz	CRAN	版本0.10.2	基于SHAP