$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
This study used only publicly available, deidentified datasets from the Gene Expression Omnibus (GEO) database. Because the work involved secondary analysis of existing public data and did not include direct participant contact, intervention, or access to identifiable personal information, additional ethics committee approval and informed consent were not required.
Data sources and preprocessing
All gene expression and single‑cell dataset were obtained from the GEO database24. For major depressive disorder, dataset GSE98793 was used, which comprises peripheral blood samples from 128 patients and 64 healthy controls. For dermatomyositis, datasets were selected based on predefined criteria, including Homo sapiens expression profiling, clearly identifiable disease and control groups, available platform annotation for probe-to-gene mapping, and suitability for discovery or validation analysis. When a GEO series contained multiple inflammatory myopathy subtypes, only dermatomyositis and normal control samples were extracted for the present study. GSE1551, GSE46239, and GSE128470 were used as the discovery/training datasets, whereas GSE5370, GSE39454, and GSE11971 were used as independent validation datasets. The dermatomyositis datasets analyzed in this study were derived mainly from affected muscle or skin tissues rather than peripheral blood. Single‑cell data for dermatomyositis were sourced from dataset GSE190510.
Raw expression matrices were downloaded from the GEO database together with the corresponding platform annotation files. Probe IDs were mapped to official gene symbols according to the manufacturer-provided GPL annotation. Probes that could not be mapped unambiguously to a single official gene symbol were removed. When multiple probes mapped to the same gene, they were collapsed at the gene level using the average expression value implemented by the `avereps` function in the limma package, thereby generating a gene-by-sample expression matrix.
To reduce intensity-dependent bias and stabilize variance, log2 transformation was applied when appropriate according to the distribution of expression values. Between-array normalization was then performed using the `normalizeBetweenArrays` function in the limma package. Missing values, when present, were imputed using K-nearest neighbor imputation. For the integrated dermatomyositis training datasets, batch correction was performed using the `ComBat` function in the sva package, with dataset/platform origin treated as the batch variable and sample group (dermatomyositis versus healthy control) included in the design matrix to preserve the biological variation of interest during batch adjustment.
All analyses were performed in R using an Integrated development environment for R on a desktop operating system. The limma package was used for probe summarization and normalization. The sva package was used for ComBat batch correction. Missing values were imputed using K-nearest neighbor imputation with k = 10.
Weighted gene co-expression network analysis
Weighted Gene Co-expression Network Analysis (WGCNA) was performed separately for the major depressive disorder and dermatomyositis datasets using the WGCNA R package25,26. Samples were hierarchically clustered using flashClust to identify outliers; samples exceeding a dendrogram height of 100 and genes in the bottom 25% of variance were excluded. For each network, a soft-thresholding power (β) was selected using pickSoftThreshold to achieve approximate scale-free topology (R2 > 0.8). The adjacency matrix was transformed into a Topological Overlap Matrix (TOM), and modules were identified via dynamic tree cutting with a minimum module size of 60 and a merge cut height of 0.2527. The WGCNA R package was used together with flashClust for hierarchical clustering. The random seed was set to 12345 for reproducibility. Module eigengenes were correlated with disease status using Pearson correlation, with P-values adjusted by the Benjamini–Hochberg method. For each disease, the module showing the strongest and most significant association with disease status was retained as the key disease-associated module. The overlap between the key module genes from the major depressive disorder dataset and those from the dermatomyositis dataset was defined as the candidate shared gene set for downstream analyses. Differential expression analysis of the integrated dermatomyositis cohort was performed separately to characterize dermatomyositis-related transcriptional changes.
Functional enrichment analysis
Gene Ontology (GO) enrichment analysis was performed using R. Gene symbols were converted to Entrez IDs using org.Hs.eg.db, and significantly enriched GO terms (p < 0.05) were identified using enrichGO in clusterProfiler. For multi‑dimensional visualization of the results, bar plots and bubble plots were generated using the enrichplot package, while a circular plot was constructed with the circlize package to display GO categories, gene counts, and enrichment factors. Legends were added with the ComplexHeatmap package. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes was also conducted in R. Gene symbols were converted to Entrez IDs based on the org.Hs.eg.db database, and significantly enriched pathways (FDR < 0.05) were identified using the enrichKEGG function from the clusterProfiler package28,29,30,31. Enrichment results were visualized using bar and bubble plots.
GeneMANIA-based functional association network analysis
Based on the previously identified shared genes, a GeneMANIA-based functional association network was constructed to explore the interaction context among these genes and their related partners. The gene list was submitted to GeneMANIA using Homo sapiens as the reference species. GeneMANIA integrates multiple evidence types, including co-expression, physical interactions, pathways, co-localization, genetic interactions, and shared protein domains. The resulting network was exported and imported into a network visualization platform for visualization and analysis. Topological analysis of the network was then performed in a network visualization platform for visualization and analysis to identify highly connected candidate nodes32,33,34.
Machine learning-based diagnostic model construction
Multiple machine-learning algorithms were used for diagnostic classification, including 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), and Partial Least Squares Regression Generalized Linear Model (plsRglm)35. A two-stage modeling framework was applied to generate 113 candidate model combinations. In the first stage, the initial algorithm was used for variable screening in the training cohort; in the second stage, the retained variables were used to fit a diagnostic classification model. Models with ≤5 selected variables were excluded from further comparison. The combined dermatomyositis datasets served as the training cohort, with labels defined as dermatomyositis versus healthy controls, whereas the independent validation cohort(s) were used for external performance evaluation. Internal resampling and tuning were algorithm-specific: glmnet-based models (LASSO, Ridge, and Elastic Net) used 10-fold cross-validation to select lambda.min; GBM used 10-fold internal cross-validation to determine the optimal number of trees; XGBoost used 5-fold resampling to select the final boosting round according to the minimum test log-loss; glmBoost used cvrisk-based internal cross-validation to determine the stopping iteration; and LDA was fitted under the caret cross-validation framework. For algorithms without explicit tuning steps in the present implementation, fixed or package-default settings were used. To reduce information leakage, feature selection, model fitting, and internal tuning were performed using the training cohort only, while the validation cohort(s) were used solely for independent prediction and AUC-based performance assessment. The caret package was used for machine-learning workflow management, with glmnet , randomForest , e1071 , gbm , xgboost , mboost , plsRglm , and MASS for individual algorithms. SHAP analysis was performed using the shapviz package . The random seed was set to 12345 before each model fitting. Models with fewer than 5 selected features were excluded. Model interpretability and gene-level contribution were further assessed using SHapley Additive exPlanations (SHAP), and the most informative genes were prioritized as candidate model-selected features for downstream biological interpretation.
Evaluation of the diagnostic performance
Receiver operating characteristic (ROC) curves were generated using the “pROC” R package to assess the diagnostic performance of candidate biomarkers. Expression levels and predictive accuracy of the candidate markers were validated in independent datasets (GSE5370, GSE11971, and GSE39454). Model performance was further assessed using confusion matrices. Differential expression of key module genes was visualized using volcano and box plots, and ROC curves were constructed to evaluate the diagnostic value of individual genes.
Gene set enrichment analysis
To explore coordinated functional changes associated with the candidate shared transcriptomic signals, Gene Set Enrichment Analysis (GSEA) was performed using clusterProfiler36,37. Gene expression data from dermatomyositis and control samples were ranked according to differential expression. Predefined gene sets corresponding to KEGG pathways (c2.cp.kegg.Hs.symbols.gmt) were used to evaluate whether genes within each pathway exhibited a coordinated trend of up- or down-regulation. Statistical significance was defined as P < 0.05.
Immune cell infiltration analysis
The normalized, log2-transformed, and batch-corrected dermatomyositis matrix was used for immune deconvolution. The CIBERSORT algorithm was applied to estimate the relative abundance of immune cell subtypes using the LM22 reference matrix38. Samples with deconvolution P < 0.05 were retained for downstream analysis. Differences in inferred immune-cell proportions between groups were visualized using box plots, and Spearman correlation analysis was conducted to assess associations between immune cell subsets and candidate shared genes.
Single-cell RNA sequencing analysis for cellular contextualization
Single-cell RNA-seq analyses were performed in R using Seurat. Harmony was used for batch correction, DoubletFinder for doublet detection, celda/decontX for ambient RNA estimation, Monocle for pseudotime trajectory analysis, CellChat for cell-cell communication analysis, AUCell for gene-set activity scoring, and GSVA for ssGSEA scoring. Raw count matrices were imported into Seurat objects with the parameters min.cells = 5 and min.features = 300. Quality-control metrics, including mitochondrial, ribosomal, and hemoglobin gene proportions, were calculated for each cell. Cells were retained only if they satisfied all of the following criteria: nFeature_RNA > 500, nCount_RNA < 5,000, percent_mito < 25, percent_ribo > 3, and percent_hb < 1. Genes detected in fewer than 3 cells were excluded. In addition, MALAT1 and mitochondrial genes were removed prior to downstream analysis. After initial filtering, doublets were identified in each sample using DoubletFinder, with PCs = 1:30 and pN = 0.25; expected doublet rates were set according to sample-specific cell numbers (<4,000 cells: 2.5%; 4,000–8,000 cells: 5%; >8,000 cells: 6.5%). Only singlets were retained. Ambient RNA contamination was further estimated using decontX, and cells with contamination scores < 0.2 were retained.
The filtered data were normalized using the LogNormalize method with a scale factor of 10,000, followed by identification of variable genes, data scaling, and principal component analysis. Batch effects across samples were corrected using Harmony with orig.ident as the batch variable. The first 15 Harmony dimensions were used for UMAP visualization and neighbor-graph construction. Clustering was performed using FindNeighbors and FindClusters, and the final clustering result was defined at a resolution of 0.05. Cell types were annotated manually according to canonical marker genes together with FindAllMarkers results39.
For downstream functional contextualization, candidate-gene activity was evaluated at the single-cell level, and the relevant immune-cell subset was subjected to trajectory and intercellular communication analyses. Pseudotime analysis was performed using Monocle with DDRTree-based dimensionality reduction followed by cell ordering. Cell-cell communication analysis was conducted using CellChat with the human ligand–receptor database, restricted to the Secreted Signaling category, and communications involving fewer than 10 cells were filtered out.
For each cell, candidate-gene activity was quantified using three complementary approaches: AUCell, ssGSEA, and AddModuleScore. AUCell scores were calculated based on gene-ranking matrices, and ssGSEA scores were generated using the GSVA framework. AddModuleScore was computed using the Seurat built-in function. The resulting AUCell, ssGSEA, and AddModuleScore values were then combined into a single score matrix. Each score type was first standardized by Z-score transformation and subsequently rescaled to a 0–1 range using min–max normalization. The final composite score (“Scoring”) for each cell was defined as the sum of the three normalized scores:
Scoring = normalized AUCell + normalized ssGSEA + normalized AddModuleScore.
For downstream subgroup analyses, the CD8⁺ T-cell subset was extracted, and cells were dichotomized according to the median Scoring value within this subset. Cells with Scoring values greater than the median were assigned to the High_Hub_genes group, whereas the remaining cells were assigned to the Low_Hub_genes group.
