$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
本研究仅使用了来自基因表达综合数据库（Gene Expression Omnibus, GEO）的公开且去标识化的数据集。由于本工作涉及对现有公开数据的二次分析，不包含对受试者的直接接触、干预或可识别个人信息的获取，因此无需额外的伦理委员会批准和知情同意。
数据来源与预处理
所有基因表达和单细胞数据集均来自GEO数据库24。对于重度抑郁症，使用了数据集GSE98793，该数据集包含来自128名患者和64名健康对照者的外周血样本。对于皮肌炎，数据集的选择基于预定义的标准，包括智人（Homo sapiens）表达谱分析、可明确区分的疾病组与对照组、可用于探针到基因映射的平台注释信息，以及适用于发现或验证分析的特性。当某个GEO系列包含多种炎症性肌病亚型时，本研究仅提取皮肌炎和正常对照样本。GSE1551、GSE46239和GSE128470被用作发现/训练数据集，而GSE5370、GSE39454和GSE11971被用作独立验证数据集。本研究分析的皮肌炎数据集主要来源于受累的肌肉或皮肤组织，而非外周血。皮肌炎的单细胞数据来源于数据集GSE190510。
原始表达矩阵及其对应的平台注释文件从GEO数据库下载。根据制造商提供的GPL注释文件，将探针ID映射到官方基因符号。无法明确映射到单一官方基因符号的探针被剔除。当多个探针映射到同一基因时，使用limma软件包中的`avereps`函数取平均表达值，将数据在基因水平上合并，从而生成基因-样本表达矩阵。
为了减少强度依赖性偏差并稳定方差，根据表达值的分布情况，适当时采用log2转换。随后使用limma软件包中的`normalizeBetweenArrays`函数进行芯片间标准化。当存在缺失值时，采用K近邻法进行填补。对于整合的皮肌炎训练数据集，使用sva软件包中的`ComBat`函数进行批次效应校正，将数据集/平台来源作为批次变量，并在设计矩阵中包含样本分组（皮肌炎 versus 健康对照），以在批次校正过程中保留感兴趣的生物学变异。
所有分析均在桌面操作系统中使用 R 的集成开发环境进行。采用 limma 软件包进行探针信号汇总和标准化。采用 sva 软件包进行 ComBat 批次效应校正。缺失值通过 K 近邻插补法填补，k 值设为 10。
加权基因共表达网络分析
分别使用WGCNA R软件包对重度抑郁症和皮肌炎数据集进行加权基因共表达网络分析（Weighted Gene Co-expression Network Analysis, WGCNA）25,26。利用flashClust对样本进行层次聚类以识别离群值；剔除树状图高度超过100的样本以及方差位于最低25%的基因。对于每个网络，使用pickSoftThreshold选择软阈值幂（β），以达到近似的无标度拓扑结构（R2 > 0.8）。将邻接矩阵转换为拓扑重叠矩阵（Topological Overlap Matrix, TOM），并通过动态树切割方法识别模块，设定最小模块大小为60，合并切割高度为0.2527。WGCNA R软件包与flashClust联合用于层次聚类分析。为保证结果可重复，随机种子设为12345。使用Pearson相关性分析将模块特征基因（module eigengenes）与疾病状态相关联，并通过Benjamini–Hochberg方法对P值进行校正。对于每种疾病，保留与疾病状态相关性最强且具有统计学意义的模块，作为关键的疾病相关模块。将重度抑郁症数据集中的关键模块基因与皮肌炎数据集中的关键模块基因之间的交集定义为候选的共有基因集，用于后续分析。对整合后的皮肌炎队列单独进行差异表达分析，以表征与皮肌炎相关的转录变化。
功能富集分析
使用 R 进行基因本体（Gene Ontology, GO）富集分析。利用 org.Hs.eg.db 将基因符号转换为 Entrez ID，并通过 clusterProfiler 中的 enrichGO 函数鉴定显著富集的 GO 条目（p < 0.05）。为了对结果进行多维度可视化，使用 enrichplot 包生成柱状图和气泡图，同时使用 circlize 包构建环形图以展示 GO 分类、基因数量及富集因子。图例由 ComplexHeatmap 包添加。差异表达基因的京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）通路富集分析也在 R 中进行。基于 org.Hs.eg.db 数据库将基因符号转换为 Entrez ID，并使用 clusterProfiler 包中的 enrichKEGG 函数鉴定显著富集的通路（FDR < 0.05）28,29,30,31。富集结果通过柱状图和气泡图进行可视化。
基于 GeneMANIA 的功能关联网络分析
基于先前鉴定的共有基因，利用 GeneMANIA 构建了功能关联网络，以探究这些基因及其相关伙伴基因之间的相互作用背景。将基因列表提交至 GeneMANIA，并以智人（Homo sapiens）作为参考物种。GeneMANIA 整合了多种证据类型，包括共表达、物理相互作用、通路、共定位、遗传相互作用以及共享的蛋白质结构域。所得网络被导出并导入网络可视化平台进行可视化与分析。随后，在网络可视化平台中对该网络进行拓扑学分析，以识别高度连接的候选节点32,33,34。
基于机器学习的诊断模型构建
采用多种机器学习算法进行诊断分类，包括随机森林（Random Forest, RF）、支持向量机（Support Vector Machine, SVM）、线性判别分析（Linear Discriminant Analysis, LDA）、朴素贝叶斯（Naive Bayes）、梯度提升机（Gradient Boosting Machine, GBM）、XGBoost、glmBoost、弹性网络（Elastic Net, Enet）、岭回归（Ridge）、最小绝对收缩与选择算子（Least Absolute Shrinkage and Selection Operator, LASSO）、逐步广义线性模型（Stepwise Generalized Linear Model, Stepglm）以及偏最小二乘回归广义线性模型（Partial Least Squares Regression Generalized Linear Model, plsRglm）35。应用两阶段建模框架生成了113种候选模型组合。在第一阶段，使用初始算法在训练队列中进行变量筛选；在第二阶段，利用保留的变量拟合诊断分类模型。所选变量数 ≤5 的模型被排除在进一步比较之外。合并的皮肌炎数据集作为训练队列，标签定义为皮肌炎与健康对照，而独立验证队列则用于外部性能评估。内部重采样与调参策略因算法而异：基于glmnet的模型（LASSO、Ridge和Elastic Net）采用10折交叉验证选择lambda.min；GBM使用10折内部交叉验证确定最优树数量；XGBoost采用5折重采样，根据最小测试对数损失选择最终提升轮次；glmBoost使用基于cvrisk的内部交叉验证确定停止迭代点；LDA则在caret交叉验证框架下拟合。对于当前实现中无明确调参步骤的算法，采用固定参数或软件包默认设置。为减少信息泄露，特征选择、模型拟合和内部调参均仅使用训练队列完成，验证队列仅用于独立预测和基于AUC的性能评估。机器学习工作流管理使用caret包，各算法分别调用glmnet、randomForest、e1071、gbm、xgboost、mboost、plsRglm和MASS包实现。SHAP分析使用shapviz包进行。每次模型拟合前将随机种子设为12345。所选特征少于5个的模型被排除。进一步采用SHapley加性解释（SHapley Additive exPlanations, SHAP）评估模型可解释性与基因水平贡献，并将最具信息量的基因优先作为候选模型筛选特征，用于下游生物学解释。
诊断性能评估
采用“pROC” R 软件包生成受试者工作特征（ROC）曲线，以评估候选生物标志物的诊断性能。在独立数据集（GSE5370、GSE11971 和 GSE39454）中验证了候选标志物的表达水平和预测准确性。进一步使用混淆矩阵评估模型性能。关键模块基因的差异表达通过火山图和箱线图进行可视化，并构建 ROC 曲线以评估单个基因的诊断价值。
基因集富集分析
为了探索与候选共有转录组信号相关的协调性功能变化，使用 clusterProfiler36,37 进行了基因集富集分析（Gene Set Enrichment Analysis, GSEA）。根据差异表达对皮肌炎和对照样本的基因表达数据进行排序。采用对应 KEGG 通路的预定义基因集（c2.cp.kegg.Hs.symbols.gmt）评估各通路内的基因是否表现出一致的上调或下调趋势。统计学显著性定义为 P < 0.05。
免疫细胞浸润分析
采用经标准化、log2转换及批次校正的皮肌炎矩阵进行免疫细胞去卷积分析。使用LM22参考矩阵通过CIBERSORT算法估算免疫细胞亚群的相对丰度38。保留去卷积 P < 0.05 的样本用于后续分析。通过箱线图可视化各组间推断的免疫细胞比例差异，并采用Spearman相关性分析评估免疫细胞亚群与候选共有基因之间的关联性。
用于细胞定位的单细胞RNA测序分析
使用 R 语言中的 Seurat 进行单细胞 RNA 测序分析。采用 Harmony 进行批次校正，DoubletFinder 用于双细胞检测，celda/decontX 用于环境 RNA 估计，Monocle 用于拟时轨迹分析，CellChat 用于细胞间通讯分析，AUCell 用于基因集活性评分，GSVA 用于 ssGSEA 评分。原始计数矩阵以 min.cells = 5 和 min.features = 300 的参数导入 Seurat 对象。为每个细胞计算质量控制指标，包括线粒体、核糖体和血红蛋白基因的比例。仅当细胞满足以下所有标准时才予以保留：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分析结果进行人工注释39。
为了进行下游功能背景分析，在单细胞水平上评估了候选基因的活性，并对相关的免疫细胞亚群进行了轨迹分析和细胞间通讯分析。使用 Monocle 进行拟时序分析，采用基于 DDRTree 的降维方法，随后进行细胞排序。使用 CellChat 进行细胞间通讯分析，数据库限定为人类配体-受体数据库中的“分泌信号通路”类别，并过滤掉涉及细胞数少于 10 个的通讯事件。
对于每个细胞，候选基因活性通过三种互补的方法进行量化：AUCell、ssGSEA 和 AddModuleScore。AUCell 得分基于基因排序矩阵计算，ssGSEA 得分使用 GSVA 框架生成，AddModuleScore 则利用 Seurat 内置函数计算。随后，将得到的 AUCell、ssGSEA 和 AddModuleScore 值整合为一个综合得分矩阵。每种得分类型首先通过 Z 分数变换进行标准化，然后使用最小-最大归一化将其重新缩放到 0–1 范围内。每个细胞的最终综合得分（“Scoring”）定义为三个归一化得分之和：
评分 = 标准化的 AUCell + 标准化的 ssGSEA + 标准化的 AddModuleScore。
在后续亚组分析中，提取CD8⁺ T细胞亚群，并根据该亚群内Scoring值的中位数将细胞分为两组。Scoring值高于中位数的细胞被归入High_Hub_genes组，其余细胞则被归入Low_Hub_genes组。
