Gene module identification concerning MGUS and MM progression was executed through independent WGCNA analyses on the GSE6477 and GSE5900 datasets. The methodology of WGCNA demonstrates exceptional proficiency in processing large-scale datasets, enabling precise identification of gene clusters substantially linked to pathological progression via sophisticated clustering and modularization methodologies [15]. Quality control and normalization procedures were first applied to ensure data consistency. Crucially, to guarantee statistical robustness, the appropriate soft-thresholding power ($\beta$) was determined based on the scale-free topology criterion (scale-free topology fit index R² $>$0.8). Utilizing this parameter, a weighted adjacency matrix was developed for each dataset and transformed into a topological overlap matrix. Module identification was achieved through hierarchical clustering coupled with dynamic tree cutting, whereby genes exhibiting similar expression profiles and potentially shared biological functions were grouped together. Finally, module–trait associations were assessed to identify modules potentially relevant to disease progression, establishing a foundation for subsequent candidate gene selection.

By integrating transcriptomic data and employing gene co-expression network analytical techniques, a comprehensive strategy was utilized to delineate pivotal genes linked to MM and MGUS progression. The investigative approach encompassed executing set-based operations across four distinct gene collections: DEGs and WGCNA module genes for both MM and MGUS. Intersection analyses were subsequently conducted, initially between MM-related DEGs and MM-related WGCNA module genes (designated as set A), and subsequently between MGUS-related DEGs and MGUS-related WGCNA module genes (designated as set B). The ultimate candidate gene collection (referred to as set C) emerged from the convergence of sets A and B. This methodological strategy was graphically represented through a Venn diagram, wherein individual closed curves depicted various gene collections, with intersecting regions emphasizing the critical genes of significant interest.

The clinical relevance of pivotal DEGs was investigated through an initial evaluation of their diagnostic potential utilizing receiver operating characteristic (ROC) analysis within the GSE6477 dataset. Correlations between gene expression patterns, patient survival outcomes, disease staging, and chromosomal aberrations were subsequently examined. Within the GSE136337 dataset, univariate Cox regression analysis identified genes demonstrating significant links to MM prognosis, while LASSO regression was employed to refine these candidates into the most informative prognostic markers. A risk score model (MGUS-related risk score, MGUSscore) was constructed based on this refined gene set, incorporating two hub genes, DAP3 and UBE2S. Patients were classified into high- and low-risk cohorts based on the median risk score, with overall survival (OS) differences evaluated through Kaplan–Meier (K–M) analysis. The prognostic value of the MGUSscore was evaluated in conjunction with clinical variables (including age, gender, International Staging System [ISS] stage, Lactate Dehydrogenase (LDH), and cytogenetics) through multivariable Cox regression analysis. Predictive performance was evaluated through time-dependent ROC curve analysis, the Concordance index (C-index), and calibration curves. A prognostic nomogram for predicting 1-, 3-, and 5-year OS was constructed utilizing the rms R package (version 6.7-0; Frank E. Harrell Jr., Vanderbilt University School of Medicine, Nashville, TN, USA). The TCGA-MM cohort was jointly analyzed with GSE136337 to facilitate core prognostic gene identification and to further evaluate the robustness and predictive performance of the MGUSscore. The detailed clinical characteristics of the GSE136337 and TCGA-MM cohorts are summarized in Table 1. External validation was subsequently performed across two independent cohorts, GSE4581 and GSE57317.

DEGs between the two cohorts were determined utilizing the ‘limma’ package, with selection criteria of AveExpr $>$5, $|$logFC$|$ $>$0.5, and p 0.05 being applied. Subsequent analyses involved Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment investigations to elucidate potential biological mechanisms and pathways. The GO enrichment approach annotated hub genes’ biological characteristics across three sub-ontological domains: molecular functions, cellular components, and biological processes. Complementarily, KEGG pathway analysis was implemented to identify potential signaling mechanisms linked to the target genes. The application of the ‘clusterProfiler’ [16] package within R facilitated comprehensive GO/KEGG enrichment investigations of intersecting genes, potentially elucidating fundamental disease progression pathways. MM subjects were systematically split into distinct expression clusters per the median MGUSscore threshold. Subsequently, gene set enrichment analysis (GSEA) executed across these stratified subgroups offered nuanced insights into gene expression variation implications.

Progressive advancements in disease research have rendered the substantial impact of immune responses on disease progression increasingly apparent. Immune scores, stromal scores, and ESTIMATE scores were computed for high and low-risk cohorts utilizing the ESTIMATE algorithm. The proportions of immune cell infiltration were determined through the ‘CIBERSORT’ [17] analytical tool. Subsequently, linear regression analysis was implemented to elucidate the associations between immune cell infiltration patterns and risk scores.

Utilizing data procured from the Genomics of Drug Sensitivity in Cancer (GDSC) database, recognized as the most comprehensive pharmacogenomics resource available, chemotherapy sensitivity predictions were generated for individual tumor specimens through implementation of the ‘oncoPredict’ [18] R package. The half-maximal inhibitory concentration (IC₅₀) values for respective pharmaceutical agents were determined via regression analysis, with predictive accuracy being assessed through 10-fold cross-validation methodology employing the GDSC training dataset.

Single-cell transcriptome datasets were procured from the GEO database, specifically the GSE176131 [19] dataset. The repository comprised data collected from 2 healthy control subjects and 9 MM patients, facilitating comprehensive single-cell correlation investigations. Single-cell RNA-seq data processing was implemented through the ‘Seurat’ [20] R package, with cellular expression matrix normalization conducted via the ‘LogNormalize’ methodology. Dimensionality reduction was achieved by principal component analysis, where the top 50 principal components were selected for cellular clustering based on marker gene identification using the elbow plot criterion. Subsequently, the expression profiles of DAP3 and UBE2S across varied cell populations were systematically evaluated.

Single-cell ATAC-seq profiling was conducted on primary specimens derived from three independent individuals: one newly diagnosed (primary) multiple myeloma patient, one relapsed/refractory multiple myeloma patient, and one healthy donor as a control. The final high-quality dataset comprised 9016 cells, consisting of 827 normal plasma cells from the healthy control, 5175 primary malignant cells from the newly diagnosed patient, and 3014 recurrent malignant cells from the relapsed patient. Reads were aligned to the hg38 genome, and peaks were identified via MACS2. Quality control retained cells with a transcription start site (TSS) enrichment score $>$6 and $>$1500 fragments. Dimensionality reduction and UMAP visualization were performed using iterative latent semantic indexing. Clusters were defined using Seurat, while gene activity scores and chromatin co-accessibility were analyzed to assess transcriptional potential and regulatory interactions. Differential accessibility was determined using ArchR (Wilcoxon test, FDR 0.05).

This investigation encompassed three newly diagnosed, treatment-naïve MM patients alongside three untreated MGUS patients serving as controls. Tissue specimens underwent processing for both immunohistochemistry (IHC) and multiplex immunofluorescence (mIF) staining to evaluate UBE2S and DAP3 expression patterns. For IHC analysis, tissue specimens were deparaffinized and rehydrated, subsequently subjected to EDTA-mediated antigen retrieval, followed by serum-based blocking procedures. Primary antibodies targeting UBE2S (1:100, ab197945, Abcam, Cambridge, UK) and DAP3 (1:100, ab302889, Abcam, Cambridge, UK) were applied, with 3,3^′-Diaminobenzidine (DAB) staining employed for visualization. Hematoxylin counterstaining was implemented, and images were acquired using light microscopy (BX53, Olympus, Tokyo, Japan).

For mIF analysis, specimens underwent EDTA antigen retrieval treatment and 3% hydrogen peroxide exposure in darkness for 25 minutes, succeeded by serum blocking. Primary antibodies underwent overnight incubation at 4 °C, after which fluorescent secondary antibodies were incubated for 1 hour. Nuclear staining was accomplished using 4^′,6-diamidino-2-phenylindole (DAPI), and the slides were subsequently mounted employing antifade mounting medium. mIF signals for UBE2S (red), DAP3 (green), and nuclei (blue) were visualized under fluorescence microscopy. Tissue specimens were subjected to IHC and mIF staining procedures utilizing commercially available reagents and antibodies, with comprehensive details provided in Supplementary Tables 1,2.

Statistical analyses were conducted utilizing R language (version 4.3.3, Foundation for Statistical Computing, Vienna, Austria), Comparisons between two groups were performed using the Wilcoxon rank-sum test, while comparisons among three or more groups were conducted using the Kruskal-Wallis test. A p 0.05 was designated as the threshold for statistical significance. The comprehensive flowchart delineating the entire bioinformatics analytical procedure is illustrated in Fig. 1. Semi-quantitative evaluation of immunofluorescence outcomes was executed through ImageJ software (National Institutes of Health, Bethesda, MD, USA). Data visualization was accomplished via GraphPad Prism (version 10; GraphPad Software, Boston, MA, USA), whereby bar graphs incorporating scatter plots were generated to depict the findings.

Fig. 1.

Flowchart illustrating the comprehensive study design, bioinformatic analysis, and experimental validation workflow. GEO, Gene Expression Omnibus; WGCNA, weighted gene co-expression network analysis; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; KM, Kaplan–Meier; IHC, immunohistochemistry; mIF, multiplex immunofluorescence; scATAC-seq, single-cell assay for transposase-accessible chromatin with sequencing.

Analysis of the GSE6477 dataset identified 961 DEGs between patients with MM and healthy controls (Fig. 2A). Conversely, the GSE5900 dataset revealed 355 DEGs between MGUS and normal controls (Fig. 2B). Heatmaps in Fig. 2C,D depict the 20 DEGs. Using WGCNA on the GSE6477 dataset, 12 gene modules significantly associated with MM were identified, each uniquely color-coded (Fig. 3A,B). The green, yellow, blue, and black modules demonstrated strong correlations with MM, with correlation coefficients and statistical significances as follows: green (r = –0.61, p = 4 $\times$ 10^-10), yellow (r = –0.86, p = 2 $\times$ 10^-27), blue (r = 0.67, p = 1 $\times$ 10^-12), and black (r = –0.46, p = 7 $\times$ 10^-6) (Fig. 3C,D). Analysis of the MGUS dataset GSE5900 also identified 12 modules, with the green module showing a significant positive correlation with MGUS (r = 0.56, p = 1 $\times$ 10^-6) (Fig. 3E,F). Focusing on the blue module linked to MM and the green module associated with MGUS, these modules were intersected with the previously identified DEGs, resulting in the identification of 12 hub genes (Fig. 3G).

Fig. 2.

DEGs in MM and MGUS cohorts. (A,B) Volcano plots depicting the differential expression analysis in the consolidated cohorts of GSE6477 (A) and GSE5900 (B). Red dots indicate significantly upregulated genes, blue dots indicate downregulated genes, and gray dots represent genes with no significant difference. (C,D) Heatmaps displaying the expression profiles of the top 20 DEGs in MM (C) and MGUS (D). Columns represent samples (Control vs. Treatment) and rows represent genes. The color scale ranges from blue (low expression) to red (high expression). DEGs, differentially expressed genes; MGUS, monoclonal gammopathy of undetermined significance.

Fig. 3.

WGCNA analysis conducted across MM and MGUS patient cohorts. (A,B) Analysis of network topology for soft-thresholding power selection in MM (A) and MGUS (B) specimens. Left panels show the scale-free fit index (y-axis) versus soft-thresholding power (x-axis); right panels display mean connectivity. (C) Hierarchical clustering dendrogram depicting MM gene organization, where distinct chromatic representations indicate separate co-expression modules. (D) Heat map visualization illustrating module-trait correlations within MM datasets. (E) Hierarchical clustering dendrogram depicting MGUS gene organization, where distinct chromatic representations indicate separate co-expression modules. (F) Heat map visualization illustrating module-trait correlations within MGUS datasets. (G) Venn diagram illustrating the intersection of genes between identified key WGCNA modules (MM_blue, MGUS_green) and differentially expressed genes (MM_diff, MGUS_diff).

Diagnostic performance of the 12 intersection genes was assessed through ROC curve analysis utilizing the GSE6477 dataset. The area under the curve (AUC) serves as a metric for evaluating diagnostic accuracy. Proximity to 1.0 in AUC values signifies superior diagnostic precision, whereas values approximating 0.5 reflect diminished discriminatory power. An AUC of 0.5 represents a lack of meaningful diagnostic differentiation. The graphical representation in Fig. 4 reveals that each of the 12 genes exhibits AUC values exceeding 0.8, thereby indicating substantial diagnostic potential for MM.

Fig. 4.

ROC analysis of 12 candidate genes in GSE6477. The AUC reflects predictive accuracy. AUC, area under the curve.

The GSE136337 and TCGA datasets, encompassing OS and status information of MM patients, were initially retrieved for analytical purposes. Within GSE136337, two patients exhibiting a survival duration of 0 were eliminated from the analytical framework. GSE136337 was designated as the training cohort for risk model development, whereas TCGA-MM, in conjunction with GSE136337, was employed to identify pivotal prognostic genes and subsequently utilized for model validation. Within the training cohort, univariate Cox regression analysis was initially executed to ascertain genes demonstrating significant association with OS (p 0.05). To mitigate potential overfitting, these survival-associated genes were subsequently processed through LASSO Cox regression incorporating 10-fold cross-validation, establishing the optimal penalty parameter ($\lambda$) and ultimately isolating two pivotal prognostic genes, DAP3 and UBE2S (Fig. 5A–D). Based on their expression profiles, an MGUSscore was computed as follows: MGUSscore = (0.211 $\times$ DAP3 expression) + (0.364 $\times$ UBE2S expression). Patients were stratified into distinct risk cohorts per the median MGUSscore. K–M survival analysis demonstrated statistically significant disparities in OS, with individuals in the high-risk cohort demonstrating markedly diminished OS versus those within the low-risk cohort (p 0.001; Fig. 5E).

Fig. 5.

Development of the prognostic model utilizing candidate genes. (A,B) Forest plots of univariate Cox analysis in GSE136337 (A) and TCGA-MM (B). Squares represent Hazard Ratios (HR) and horizontal lines indicate 95% Confidence Intervals (CI). (C,D) LASSO regression analysis in GSE136337 (C) and TCGA-MM (D). Left panels: coefficient profiles; Right panels: parameter selection via 10-fold cross-validation. (E) K-M survival curves of MM patients stratified by MGUSscore in GSE136337. Patients were stratified by median risk score into high-risk (red) and low-risk (blue) groups.

Significant differences were observed regarding ISS stage, LDH, $\beta$2-microglobulin, albumin, del1p, del13q, and del1p32 (Fig. 6A,B). The predictive performance of the risk model was evaluated using time-dependent ROC curve analysis based on the MGUSscore (Fig. 6C). Furthermore, multivariate Cox regression analysis demonstrated that the MGUSscore serves as an independent prognostic factor for patients with MM (hazard ratio: 2.251, 95% confidence interval: 1.570–3.228, p 0.001) (Fig. 6D). A comprehensive nomogram integrating the risk score, patient demographics, ISS stage, LDH, and del1p was constructed to predict OS at 1-, 3-, and 5-year intervals (Fig. 6E). The concordance index (C-index) of the model was 0.673. Internal validation using 1000 bootstrap resamples confirmed the model’s reliability and discriminative ability. Calibration curves at 1, 3, and 5 years indicated good agreement between predicted and observed survival, validating the model’s stability and accuracy (Fig. 6F).

Fig. 6.

Clinical relevance and prognostic performance of MGUSscore. (A) Correlation between disease stage, $\beta$2-microglobulin levels, albumin concentrations, LDH activity, and MGUSscore. *p 0.05; **p 0.01. (B) Association between MGUSscore and chromosomal aberrations, including del1p, del13q, and del1p32. **p 0.01. (C) Time-dependent ROC curves for 1-, 3-, and 5-year OS based on MGUSscore in GSE136337. (D) Multivariable Cox proportional hazards regression analysis incorporating MGUSscore and clinical parameters (MGUSscore, age, gender, ISS stage, LDH, and del1p). (E) Predictive nomogram incorporating MGUSscore, age, gender, ISS stage, LDH, and del1p for estimating 1-, 3-, and 5-year OS probabilities. (F) Calibration plots of the nomogram for predicting 1-, 3-, and 5-year OS accuracy.

The prognostic effectiveness of the MGUSscore received subsequent validation through three independent cohorts. K–M survival analyses uniformly demonstrated that patients classified in the high-risk cohort exhibited notably diminished OS versus their low-risk counterparts across all datasets (Fig. 7A–C). Time-dependent ROC curve analyses further substantiated the model’s generalizability. In the TCGA-MM cohort, AUC values for 1-, 3-, and 5-year OS were ascertained as 0.660, 0.690, and 0.690, respectively (Fig. 7D). For the GSE4581 cohort, AUC values for 1-, 3-, and 5-year OS were procured as 0.661, 0.678, and 0.697, respectively (Fig. 7E). Likewise, in the GSE57317 cohort, AUC values for 1-, 2-, and 3-year OS all exceeded 0.7 (Fig. 7F).

Fig. 7.

Prognostic validation of MGUSscore across multiple cohorts. (A–C) Kaplan-Meier overall survival curves for the TCGA-MM (A), GSE4581 (B), and GSE57317 (C) cohorts. Patients were stratified into high-risk (red) and low-risk (blue) groups based on the median MGUSscore. p-values were calculated using the log-rank test. (D–F) Time-dependent ROC curves evaluating the prognostic accuracy of MGUSscore. Panels display the Area Under the Curve (AUC) for predicting: (D) 1-, 3-, and 5-year overall survival in TCGA-MM; (E) 1-, 3-, and 5-year overall survival in GSE4581; (F) 1-, 2-, and 3-year overall survival in GSE57317.

Differential expression validation of DAP3 and UBE2S genes through the GSE6477 dataset demonstrated statistically significant elevation in MM patient populations relative to MGUS subjects and normal individuals (Fig. 8A). Subsequent clinical investigations utilizing the GSE136337 dataset indicated that elevated genetic expression of these specific genes linked to diminished OS and progression-free survival (PFS) parameters (Fig. 8B,C). Moreover, DAP3 and UBE2S exhibited a link to the ISS of MM (Fig. 9A). Additionally, elevated UBE2S expression was linked to chromosomal abnormalities, including 1 del1p, del13q, and del1p32, whereas DAP3 was linked exclusively to del13q (Fig. 9B–D). These observations suggest that DAP3 and UBE2S may contribute to the malignant progression of MM.

Fig. 8.

Expression patterns and prognostic implications of DAP3 and UBE2S in MM. (A) Boxplots comparing the expression levels of DAP3 and UBE2S across Normal (blue), MGUS (red), and MM (green) groups in the GSE6477 dataset. (B) K–M survival curves demonstrating the correlation between DAP3 and UBE2S expression levels and OS in MM patients derived from the GSE136337 cohort. Patients were stratified into high-expression (red) and low-expression (blue) groups based on optimal cut-off values. (C) K–M survival curves depicting the link between DAP3 and UBE2S expression levels and PFS in MM patients derived from the GSE136337 cohort. Color coding and stratification are consistent with (B).

Fig. 9.

Correlation of DAP3 and UBE2S expression levels with clinical and cytogenetic characteristics in MM. (A) Boxplots comparing DAP3 (left) and UBE2S (right) expression levels across International Staging System (ISS) stages: Stage I (red), Stage II (green), and Stage III (blue). (B–D) Boxplots comparing gene expression between patients with and without specific chromosomal deletions: (B) del(1p), (C) del(13q), and (D) del(1p32). Patients are stratified by the absence (red, FALSE) or presence (blue, TRUE) of the deletion. Statistical significance: *p 0.05; **p 0.01; ***p 0.001; NS = not significant.

Comparative expression analysis executed between the two risk cohorts identified a sum of 108 DEGs ($|$log₂FC$|$ $>$0.5, p 0.05), consisting of 70 upregulated and 38 downregulated genes (Fig. 10A,B). KEGG pathway analysis demonstrated enrichment in pathways associated with cell cycle regulation, p53 signaling cascades, and drug resistance mechanisms, encompassing platinum drug resistance and antifolate resistance, thereby emphasizing the potential contributions of these genes to cancer therapeutic response and resistance mechanisms (Fig. 10C). GO enrichment analysis demonstrated that these DEGs participate in critical biological processes like cell division, microtubule cytoskeleton organization, and chromosome segregation, thereby highlighting their involvement in cellular proliferation and mitotic processes, particularly within cell cycle regulation (Fig. 10D). Interestingly, GSEA analysis, revealed substantial enrichment in gene sets related to ribosome function (e.g., ribosomal subunit) and mitochondrial components (e.g., mitochondrial protein-containing complex). These findings suggest a highly active cellular state in the high-risk cohort, highlighting a nexus between protein synthesis machinery and cellular bioenergetics (Fig. 10E).

Fig. 10.

Differential expression and functional enrichment analysis between high- and low-MGUSscore cohorts. (A) Volcano plot displaying the DEGs between the high and low MGUSscore cohorts (Red: upregulated; Blue: downregulated; Gray: not significant). (B) Heatmap of DEGs. The color gradient ranges from blue (low expression) to red (high expression). (C,D) Bar charts displaying KEGG pathway (C) and GO term (D) enrichment results. (E) Gene Set Enrichment Analysis (GSEA) plot showing enriched pathways in the high-risk group.

Spearman correlation analysis demonstrated a notable link between the two core genes and infiltrating immune cells present in the tumor microenvironment (TME), with particular emphasis on plasma cells (Fig. 11A,B). Subsequent analysis demonstrated that elevated MGUSscore was markedly associated with enhanced infiltration of plasma cells, and resting memory CD4⁺ T cells, while exhibiting decreased infiltration of CD8⁺ T cells and M1 macrophages (Fig. 11C). The high-MGUS score cohort demonstrated reduced stromal, immune, and ESTIMATE scores, accompanied by an elevated tumor purity score, suggesting enhanced tumor purity and deteriorated prognosis (Fig. 11D).

Fig. 11.

Immune infiltration and drug sensitivity analysis in high- and low-MGUSscore cohorts. (A,B) Lollipop charts showing correlations between immune cells and DAP3 (A) or UBE2S (B). Dot size represents correlation coefficient; color indicates p-value. (C) Boxplots comparing immune cell infiltration between high-risk (red) and low-risk (blue) groups. (D) Boxplots of tumor microenvironment scores (ESTIMATE, Stromal, Immune, Tumor Purity) in high-risk (red) and low-risk (blue) groups. p-values were calculated by Kruskal-Wallis test. (E) Violin plots of predicted IC₅₀ values for Venetoclax, Bortezomib, and Cyclophosphamide. Red: high-risk; Blue: low-risk. (F) Heatmap showing correlations between gene expression and drug sensitivity. Color gradient ranges from blue (negative correlation) to red (positive correlation). *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns, not significant.

Drug sensitivity predictions for individual patients were generated utilizing sensitivity data derived from the GDSC database. Within the high-risk cohort, reduced IC₅₀ values were observed for Venetoclax, Bortezomib, and Cyclophosphamide (Fig. 11E). The graphical representations demonstrate differential sensitivity patterns to Bortezomib across distinct risk stratifications, wherein diminished IC₅₀ values typically indicate enhanced drug susceptibility. Such comparative analysis enables the evaluation of therapeutic response variations between high- and low-risk patient populations, thereby establishing a foundation for personalized treatment protocols. Furthermore, the correlation between the two core genes and frequently employed anti-neoplastic agents was investigated. The findings revealed that expression levels of both genes exhibited negative correlations with the IC₅₀ values of Bortezomib (Fig. 11F), indicating that higher expression of these genes is associated with increased sensitivity to the drug.

Single-cell data procured from GSE176131 underwent analysis utilizing the Seurat package. Cell clustering was executed through the UMAP algorithm, with subsequent annotation into B cells, CD4⁺ T cells, CD8⁺ T cells, myeloid cells, immature red cells, and plasma cells to facilitate gene expression observation [19] (Fig. 12A). Predominant expression of DAP3 and UBE2S was observed within plasma cells and CD8⁺ T cells (Fig. 12B). Subsequently, plasma cells and CD8⁺ T cells were extracted for additional analysis, with core genes being visualized through UMAP plots (Fig. 12C). Expression levels of DAP3 and UBE2S within plasma cells and CD8⁺ T cells demonstrated elevation corresponding to MM grading progression (Fig. 12D–G). Trajectory analysis was conducted using Monocle3 [21, 22, 23] within the identical low-dimensional space to represent potential developmental pathways of these cellular populations, thereby illustrating DAP3 and UBE2S expression patterns throughout cellular development (Fig. 12H). Validation through the Human Protein Atlas [24] demonstrated substantial expression of DAP3 and UBE2S within plasma cells, thereby corroborating previously obtained findings (Fig. 12I,J).

Fig. 12.

Single-cell analysis of DAP3 and UBE2S in MM. (A) Annotated UMAP plot of cell clusters. (B) Dot plot of gene expression (Dot size: % expressed; Color: average expression). (C) UMAP visualization depicting core gene expression within isolated plasma cells and CD8⁺ T cells. (D–G) Bar plots of UBE2S and DAP3 levels in CD8⁺ T cells (D,E) and Plasma cells (F,G) across Normal and MM groups (***p 0.001). (H) Trajectory analysis of expression dynamics. (I,J) Validation UMAP plots from the Human Protein Atlas for DAP3 (I) and UBE2S (J).

Single-cell ATAC-seq profiles comprising 827 normal, 5175 primary, and 3014 recurrent plasma cells were analyzed to delineate epigenetic landscapes. Dimensionality reduction utilizing UMAP revealed significant heterogeneity, with malignant populations clustering distinctly from normal controls, indicating profound chromatin remodeling (Fig. 13A,B). Assessment of chromatin-inferred gene activity demonstrated sustained DAP3 accessibility across malignant clusters (Fig. 13C), whereas UBE2S exhibited marked upregulation in myeloma cells relative to controls (Fig. 13D). Quantitatively, a distinct “Normal Recurrent Primary” accessibility trajectory was delineated, evidencing stage-dependent epigenetic modulation of UBE2S throughout disease progression. Genomic track analysis was subsequently performed to elucidate regulatory architectures. The DAP3 locus exhibited a consistent chromatin architecture, characterized by synchronized open chromatin signals at the DAP3 and adjacent YY1AP1 promoters (Fig. 13E). This conserved co-accessibility implies that these genes are governed by shared cis-regulatory elements or reside within a common active topologically associating domain (TAD). Conversely, the UBE2S locus displayed a coordinated regulatory unit involving a persistently open upstream element and the UBE2S promoter (Fig. 13F). While the upstream element remained stable, specific signal enrichment at the UBE2S promoter was observed in malignant samples, indicating an epigenetic shift favoring enhanced transcriptional potential.

Fig. 13.

Single-cell chromatin accessibility landscape of DAP3 and UBE2S. (A) UMAP plot colored by sample group: Normal (blue), Primary (red), and Recurrent (green). (B) UMAP plot colored by annotation: Plasma (red) and Malignant Plasma (blue). (C,D) Feature plots displaying chromatin accessibility for DAP3 (C) and UBE2S (D); color gradient from blue (low) to yellow (high) indicates signal intensity. (E) Normalized signal tracks at the DAP3 locus (chr1:155,637,966–155,737,967) for Normal (red), Primary (blue), and Recurrent (green) samples. (F) Normalized signal tracks at the UBE2S locus (chr19:55,357,776–55,457,777), formatted as in (E).

Both IHC and mIF analyses corroborated the observations that UBE2S and DAP3 expression levels were markedly elevated in MM tissues versus MGUS controls, demonstrating statistical significance. IHC analyses demonstrated intensified DAB staining patterns for UBE2S and DAP3 within MM specimens, whereas mIF exhibited more pronounced red (UBE2S) and green (DAP3) fluorescent signals in MM tissues (Fig. 14A,B and Fig. 15A,B). Subcellular localization analysis revealed that DAP3 was predominantly located in the cytoplasm, whereas UBE2S showed diffuse expression in both the nucleus and cytoplasm. Semi-quantitative assessment of mIF additionally validated that UBE2S and DAP3 expression levels were substantially upregulated in MM versus MGUS (p 0.01, respectively) (Fig. 15C). These findings indicate that UBE2S and DAP3 may contribute essential functions in MM progression and merit additional investigation.

Fig. 14.

IHC validation of UBE2S and DAP3 in MGUS and MM tissues. (A,B) Representative immunohistochemical staining patterns for UBE2S (A) and DAP3 (B) within MGUS and MM tissue specimens. In each panel, the upper image shows a low-magnification view, and the lower image displays a higher-magnification view of the boxed region. Brown staining indicates positive protein expression. Scale bars: 50 µm (main images) and 20 µm (magnified inserts).

Fig. 15.

mIF validation of DAP3 and UBE2S expression in MGUS and MM tissues. (A,B) Representative immunofluorescence images of MGUS (A) and MM (B) tissues. The left panels show individual channels for UBE2S (red), DAP3 (green), and DAPI (blue); the center panels show the merged view; the right panels show high-magnification insets of the merged view. DAP3 exhibits predominantly cytoplasmic localization, while UBE2S displays diffuse expression in both the nucleus and cytoplasm. (C) Quantitative assessment of relative fluorescence intensity for UBE2S (left) and DAP3 (right). *p 0.05. Scale bars: 20 µm (center merged panels; applies to left single-channel panels) and 5 µm (right magnified insets).