We retrieved the single cell RNA sequencing (scRNA-seq) data from Gene Expression Omnibus via the access number GSE189460 (8), which includes pre-treatment bone marrow specimens from 18 patients with multiple myeloma (MM) who underwent bortezomib–melphalan–prednisone therapy. Based on clinical response, 10 of these 18 samples were optimal responders (OR), while the rest 8 samples were suboptimal responders (SOR).

We also retrieved bulk mRNA-seq data from three studies, including two from GEO and one from MMRF database. Specifically, mRNA-seq data, together with complete survival records, were retrieved from a study of 113 responders (R) and 126 non-responders (NR) who were enrolled in bortezomib (PS-341) clinical trial for MM treatment (8). Two sets of transcriptome data are available given that two Affymetrix platforms (GPL96 [HG-U133A] and GPL97 [HG-U133B]) were applied for mRNA sequencing two sets of transcriptome data are available. Meanwhile, the mRNA-seq data from 558 cases who enrolled in Total Therapy 2 (TT2) and Total Therapy 3 (TT3) was retrieved GEO (GSE2658) (14). In addition, we also obtained mRNA-seq profiles for 764 MM patients from CoMMpass (MMRF) database deposited in the Multiple Myeloma Research Foundation (https://research.themmrf.org).

Single-cell preprocessing was performed in Seurat v4.1.1. Cells exceeding 10% mitochondrial gene content or 5% hemoglobin gene expression, or expressing fewer than 200 or more than 5,000 genes, were filtered out in line with previous studies (15–17). Doublets were detected and removed using the DoubletFinder package. Batch effects were corrected via the RunHarmony function in the harmony package (18). Prior to differential, enrichment and statistical analyses, data normalization (i.e., log-transformation), clustering, and dimensionality reduction analyses were performed to all scRNA-seq data using Seurat. We use FindVariableFeatures to select genes (n=2000) and conducted principal component analysis (PCA, via RunPCA), in which the top 20 PCs were retained for the following analysis. Cell clustering (resolution: 0.6) was conducted using FindClusters, with the identified clusters being annotated according to marker genes from CellMark2.0 and well-characterized lineage markers.

Differentially expressed genes (DEGs) for each cluster were identified using Seurat’s FindMarkers function. In our study, we focused on upregulated genes in downstream analyses. The potential biological pathways of these genes were identified via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using clusterProfiler v4.8.2 (19).

We applied Infercnvpy with default parameters to distinguish malignant from non-malignant plasma cells. Normal cells from the tumor microenvironment (TME), specifically T/NK cells, were used as the reference population.

We utilized the “Scissor” (v 2.0.0) (20) algorithm to link the survival status of 764 bulk transcriptomic samples with complete transcriptome data (expression data using log2(fpkm+1)) and complete survival status in the MMRF_COMMPASS dataset to single-cell data of multiple myeloma. Patients with OS status of death were classified as worse status, while those with OS status of alive were classified as good status. The Scissor function was run on epithelial cells with the following parameter settings: alpha=0.05, family = “cox”. Scissor+ cells were associated with worse status, while Scissor- cells were associated with good status.

SCENIC analysis was performed using the pySCENIC v0.12.1 pipeline (21) to infer regulon activity scores (RAS) in malignant plasma cells. GRNBoost2 was used to infer co-expression modules of transcription factors (TFs) and candidate targets. RcisTarget identified enriched DNA motifs within these modules, defining each TF and its direct targets as a regulon. Regulon activity per cell was quantified using AUCell.

Monocle v2.28.0 (22) was used to reconstruct differentiation trajectories among malignant plasma cells. Following dimensionality reduction and cell ordering, cells were mapped onto branched trajectories. Branch Expression Analysis Modeling (BEAM) was then used to identify genes exhibiting branch-dependent expression dynamics, shedding light on fate decision mechanisms.

The CellChat v1.6.1 algorithm (23) was applied to a merged Seurat object containing malignant plasma cells and other TME populations. After constructing the CellChat object with a curated ligand–receptor database, we used computeCommunProb and computeCommunProbPathway to infer the interaction probabilities at both individual receptor–ligand and signaling pathway levels.

We scored malignant plasma cells for “cancer stemness” and T cells for “cytotoxicity” using AUCell. Gene sets were sourced from the Molecular Signatures Database (MsigDB v2023.1; https://www.gsea-msigdb.org/gsea/index.jsp), with the T-cell cytotoxicity signature detailed in Supplementary Table S1 .

We dichotomized samples into PTPRG-high and PTPRG-low groups based on the determined cut-offs (via surv_cutpoint in survminer package). The survival rate of these two groups and comparisons were visualized using Kaplan–Meier curves.

The multiple myeloma cell lines U266 and NCI-H929 were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (Gibco), 100 U/mL penicillin G, and 100 μg/mL streptomycin at 37°C in a humidified incubator containing 5% CO₂. Cells in the logarithmic growth phase (~4 × 10⁵) were transfected with 5 nM gene-specific siRNA or negative control siRNA (si-NC; GenePharma, Shanghai) using Lipofectamine 3000 (Invitrogen), following the manufacturer’s protocol. Knockdown efficiency was confirmed by quantitative reverse transcription PCR (qRT-PCR); Primer sequences are listed in Supplementary Table S2 .

Total RNA was extracted from MM cells using TRIzol reagent (Thermo Fisher Scientific) and reverse-transcribed using PrimeScript™ RT (Takara). Quantitative reverse transcription PCR (qRT-PCR) was performed using HiScript II Q RT SuperMix (TRANS, AU341) to assess the expression of target genes, with β-actin used as the internal control. Each experiment included three biological replicates, each with technical triplicates. Primer sequences used in this study are listed in Supplementary Table S3 .

We used the CCK8 assay to detect the viability of cells in accordance with the manufacturer’s protocol. U266 and NCI-H929 cells transfected with siNC, siPTPRG-1, and siPTPRG-2 were seeded into 96-well plates at a density of 5,000 cells per well. At 0, 24, 48, 72, and 96 hours, 10 μL of CCK-8 solution (Biosharp, Shanghai, China) was added to each well, followed by incubation for 1 hour at 37°C. Absorbance was then measured at 450 nm using a microplate reader (BD Biosciences, USA). The data were analyzed and visualized using GraphPad Prism software.

Apoptosis was assessed using an Annexin V-FITC/PI apoptosis detection kit (BD Biosciences, Franklin Lakes, NJ, USA). U266 and NCI-H929 cells were harvested 24 hours after siRNA treatment and resuspended in 1× binding buffer. A total of 100 μL of the cell suspension was incubated with 5 μL Annexin V-FITC and 2.5 μL propidium iodide (PI) for 30 minutes at 37°C in the dark. Samples were then analyzed using a BD FACSCanto II flow cytometer (BD Biosciences). Cells positive for Annexin V but negative for PI were considered early apoptotic, while double-positive cells were considered late apoptotic or necrotic. Flow cytometry data were analyzed using FlowJo software (BD Biosciences).

In order to analyze the effect of PTRPG on downstream pathway proteins, cells and tissues were lysed with RIPA lysis buffer containing 1% PMSF. The lysates were then centrifuged, and the supernatant was collected. The quantified protein supernatant was supplemented with 4× protein loading buffer proportionally, boiled for 10 min to denature the protein, and stored at −80°C. Proteins were then separated by 10% SDS-PAGE and electrophoretically transferred onto polyvinylidene fluoride membranes where they were blocked with 5% skim milk and incubated with β-actin (Proteintech, 20536-1-AP, 1:10,000), anti-PTRPG (Abclonal, A14253, 1:1000), caspase-3 (CST, 24232, 1:800), and cleaved caspase-3 (Proteintech, 68773-1, 1:3,000). overnight at 4°C. Next, the membranes were incubated with horseradish peroxidase-conjugated anti-rabbit IgG. Antigen–antibody complexes were then detected with enhanced chemiluminescence reagent. The resulting Images were processed and analyzed using ImageJ software.

We performed a minimum of three independent biological replicates for all experiments. Continuous variables were compared using the Mann–Whitney U test (two groups) or Kruskal–Wallis test (more than two groups); categorical variables were assessed by χ² test. All analyses were conducted in R v4.0.5. Two-tailed p < 0.05 was considered statistically significant otherwise stated.

A total of 18 samples with single-cell data were obtained from the GSE189460 dataset ( Figure 1 , Figure 2A ), including 10 OR tissue samples and 8 SOR tissue samples ( Figure 2B , Supplementary Table S4 ). After data quality control, a total of 103,171 single cells were retained. Dimensionality reduction and clustering analysis identified 27 cell clusters ( Supplementary Figure S1A ). Based on the expression of cell marker genes, 11 cell types were annotated as plasma cells, T/NK cells, CD14+ monocytes, B cells, CD16+ monocytes, Proliferative cells, Hematopoietic Stem and Progenitor Cells (HSPCs), Megakaryocytes, myeloid Dendritic Cells (mDCs), plasmacytoid Dendritic Cells (pDCs), and Pre-B cells ( Figure 2C, D ). The three most abundant cell types were Plasma cells, T/NK cells, and CD14+ monocytes ( Figure 2C , Supplementary Figure S1B ). The calculated proportion of each cell type was consistent across examined samples indicating a good data integration ( Supplementary Figure S1C ). Analysis of the changes in the proportion of various cell types between the two groups revealed that the proportion of plasma cells in the SOR group was significantly higher than in the OR group, while the proportions of T/NK cells, CD14+ monocytes, pre-B cells, and HSPCs were significantly lower in the SOR group ( Figure 2E , Supplementary Figure S1D ). Furthermore, the copy number variation (CNV) score for each plasma cell was calculated using inferCNV ( Supplementary Figure S1E ), leading to the identification of 35,944 malignant plasma cells ( Figure 2F , Supplementary Figure S1F ).

We then focused on malignant plasma cells to explore the potential heterogeneity. By conducting clustering analyses, the identified malignant plasma cells were reclassified into 5 cell subclusters ( Figure 3A ). By correlating with prognostic phenotypes from bulk transcriptome data, we identified a subcluster containing 1,331 malignant plasma cells associated with poor prognosis (i.e., worse survival) while a subcluster containing 3,520 malignant plasma cells shows good prognosis (i.e., good survival) ( Figure 3B, D ). Cell proportion analysis revealed that in malignant plasma cell subcluster 3 (MalPlasma3), the proportion of cells from the SOR group was higher than that from the OR group ( Figure 3C, E ). Furthermore, 93.1% of the worse survival malignant plasma cells originated from the MalPlasma3 subcluster ( Figure 3F ), indicating that the MalPlasma3 subcluster might be a core driver of poor prognosis, being enriched in the SOR group and highly correlated with treatment resistance. Additionally, worse survival malignant plasma cells were exclusively present in the SOR group ( Figure 3G, H ), suggesting these cells might be the direct cause of treatment failure.

We then conducted differential gene expression analysis to identify genes differing across malignant plasma cell subclusters ( Figure 4A ) and between worse survival and good survival malignant plasma cells ( Figure 4E ) (padj<0.05, |log2FC|>0.25) and conducted functional enrichment analysis to reveal potential pathways these differential genes involved ( Figure 4B, C, F, G ). KEGG results indicated that the Stemness up pathway was significantly activated in the MalPlasma3 subcluster and worse survival malignant plasma cells. HALLMARK results showed significant activation of E2F targets, MYC targets V1/V2, and G2M checkpoint pathways in the MalPlasma3 subcluster and worse survival malignant plasma cells. These pathways are closely related to cell stemness, cell cycle, and cell proliferation and differentiation. These findings collectively suggest that the MalPlasma3 subcluster and worse survival malignant plasma cells play important roles in promoting cell proliferation, differentiation, and invasion, which directly leads to patient drug resistance and poor prognosis. Furthermore, transcription factor analysis ( Figure 4D ) showed that in the malignant cell MalPlasma3 subcluster, the top 5 transcription factors with the highest activity were E2F8, E2F7, FOXM1, E2F1, and TIMELESS, a group of regulators known to be related to cell proliferation, invasion, and cell cycle. These results suggested that worse survival malignant plasma cells exhibit significant characteristics in terms of cell stemness and cycle regulation, leading to poor prognosis. Therefore, AUCELL analysis was further used to assess the activity of the Cancer stemness pathway signature in worse survival, good survival, and background malignant plasma cells. The results showed that the Cancer stemness signature activity was significantly higher in worse survival malignant plasma cells than in good survival and background malignant plasma cells ( Figure 4H ), whereas the differences between good survival and background malignant plasma cells were nonsignificant ( Figure 4H ).

We then investigated the dynamic evolution process of malignant plasma cells through psuedotime analysis. As shown in Fig4, malignant plasma cells exhibit 5 different states: cells in state 1 were considered potential starting points, followed by a bifurcation at branch point 1, where cells in state 2 developed towards the left of the trajectory, and cells in state 5 developed towards the right ( Figure 5A, B ). Furthermore, we found that more SOR group malignant plasma cells and worse survival malignant plasma cells were located at the terminal end of the differentiation time after branch point 1 as visually documented in cell trajectory plot ( Figure 5C, D ) and ridge plot ( Figure 5E ). Using the branched expression analysis modeling (BEAM), we identified 50 branch-dependent genes that play key roles in regulating cell differentiation from pre-branch to post-branch (Cell fate 1, Cell fate 2). Based on expression similarity, these genes were further divided into 6 modules (clusters). As shown in Figure 5F , Cluster 1 exhibited an overall upward trend in gene expression levels during differentiation from pre-branch towards the left of the trajectory after branch point 1 (Cell fate 1). Furthermore, HALLMARK enrichment analysis showed that Cluster 1 genes were mainly enriched in pathways related to cell cycle, proliferation, and differentiation, such as G2M checkpoint and E2F targets ( Figure 5G ). Consistent with the previous differential gene enrichment results, these results further indicate that malignant plasma cells gradually differentiate into a worse survival phenotype, which is accompanied by dramatically enhanced abilities in proliferation, differentiation, and invasion, thus affecting drug efficacy and leading to patient drug resistance.

To further investigate the potential role of tumor microenvironment (TME) in MM, we further reclustered T/NK cells using unsupervised dimensionality reduction and clustering analysis. We found that T/NK cells were grouped into 12 cell clusters ( Supplementary Figure S2A ). Based on gene expression of cell marker, a total of 7 cell subclusters were identified, including Natural Killer cells (NK), Effector T cells (Effect T), Helper T cells (Th), Naïve T cells, Memory T cells, Regulatory T cells (Treg), and Interferon T cells (IFN T) ( Figure 6A, C ). Notably, all T/NK cell subclusters were shared between the OR and SOR groups, but they exhibited heterogeneous cell proportions ( Supplementary Figure S2B, D-E ). Analysis of cell counts and proportions between these two groups revealed that the top 3 most abundant cell types were NK cells, T-effect cells, and T naïve cells, respectively ( Supplementary Figure S2C ). Furthermore, Naïve T and Th cells showed a significant increase in proportion in the OR group, while Memory T cells showed a significant decrease in proportion in the OR group ( Figure 6B, D , Supplementary Figure S2F ). This finding suggests that remodeling of the immune microenvironment may be a key factor for good treatment response. It has been shown that Naïve T cells are unactivated T cells with high proliferative potential and differentiation capacity, while Th cells can enhance anti-tumor immune responses. In the OR group, the increased proportion of Naïve T and Th cells may reflect stronger anti-tumor immune potential and indicate that patients are more sensitive to the bortezomib-melphalan-prednisone regimen. Given that memory T cells have long-term survival and rapid response capabilities to antigens, the significant decrease in the proportion of these cells in the OR group may reflect their effectively activation and differentiation into effector T cells while reduced proportions of Memory T cells during treatment. Further scoring analysis of the cytotoxic features of T/NK cell subclusters ( Figure 6E ) showed that compared to the SOR group, the cytotoxicity feature scores of NK cells, Effect T cells, and Naïve T cells were all significantly elevated in the OR groups, indicating that these cells have stronger cytotoxicity and immune killing effects in the OR group, enabling them to effectively attacking and clearing tumor cells, and exhibiting stronger sensitivity to drugs.

To better interpret the communication between cellular components in TME, we constructed cell interaction networks of potential receptor-ligand pairs for the OR and SOR groups, respectively ( Figure 7A, B ). We observed that compared to the OR group, the communication between different cellular components in SOR group samples varied considerably. Given the observed significant activation of pathways related to cell stemness, cell cycle, and differentiation in the MalPlasma3 subcluster and worse survival malignant plasma cells, we thus focused on the interaction of cells in this subcluster with HSPCs. We found that MalPlasma3 subcluster, acting as signal-sending cells, presented significant larger number of interactions with HSPCs in OR groups as compared to the interactions in SOR group ( Figure 7C, D ). Further comparative analysis of ligand-receptor pairs between the two groups revealed that in the OR group, the MalPlasma3 subcluster specifically regulated HSPCs via IGF1-IGF1R, a signaling pathway regulating cell growth, proliferation, and apoptosis by activating signaling pathways such as PI3K/Akt and MAPK ( Figure 7E ). These results were echoed with previously observed increased proportion of HSPCs in the OR group, reflecting a relatively healthier bone marrow microenvironment, but stronger hematopoietic/immune reconstruction ability, and potential anti-tumor immune regulatory effects in this group. In addition, we found that in the OR group, Memory T, NK, and Effect T cells specifically regulated the MalPlasma3 subcluster via IFNG-(IFNGR1+IFNGR2) ( Figure 7F ). Furthermore, in the OR group, Pre-B cells were found to specifically regulate the MalPlasma3 subcluster via CCL28-CCR10.

We further leveraged bulk transcriptomics data from GSE9782 (parallel detection using GPL96 and GPL97 platforms) to assess the relationship between tumor stemness-related genes and prognosis. As shown in Fig8A, a total of 2298 differentially expressed genes (DEGs) were identified in the GPL96 platform data in tumor samples, of which 785 genes were upregulated and 1513 genes were downregulated. In the GPL97 platform data, 1119 DEGs were identified in tumor samples, with 613 genes upregulated and 506 genes downregulated ( Figure 8A ) (P<0.05). Furthermore, we mapped the upregulated DEGs from the MalPlasma3 subcluster (compared to other malignant plasma cell subclusters, a total of 1485 upregulated DEGs) and worse survival malignant plasma cells (compared to good survival malignant plasma cells, a total of 2679 up DEGs) with Cancer stemness pathway genes. We found that Cancer stemness pathway related DEGs were upregulated, including 449 genes in the MalPlasma3 subcluster and 693 genes in worse survival malignant plasma cells. Finally, these genes were overlapped with the upregulated DEGs from the two GSE9782 datasets, resulting in 7 tumor stemness-related genes, including NONO, CBX3, SLC25A3, PTPRG, NPM1, HINT1, HNRNPA1 ( Figure 8B ). Among the 7 genes, PTPRG was specifically highly expressed only in MalPlasma3 and MalPlasma5 subclusters ( Figure 8D ). Considering factors such as whether candidate genes were highly expressed in the responder group (i.e., R group) of the bulk dataset (GSE9782) and in worse survival malignant plasma cells, as well as poor prognosis, PTPRG was ultimately selected. Notably, we found that PTPRG was highly expressed in the NR group of the GSE9782 dataset (GPL96 platform) ( Figure 8C ) and in worse survival malignant plasma cells ( Figure 8E ). Further prognostic validation of PTPRG using bulk datasets GSE9782(best cutoff:6.8011457193872), GSE2658 (best cutoff:9.417325114), and MMRF-CoMMpass (best cutoff: 2.119787759) revealed that samples with high PTPRG expression had significantly poorer prognosis ( Figure 8F , Supplementary Figure S3 ).

To better reveal potential functional role of PTPRG in MM, we performed in vitro knockdown studies using two distinct siRNAs in U266 and NCI-H929 cell lines. Effective suppression of PTPRG expression was confirmed by qRT-PCR ( Figure 9A, B ), and subsequent CCK-8 assays revealed a significant reduction in cell viability upon PTPRG knockdown ( Figure 9C, D ). Consistently, both flow cytometric analysis and Western blotting demonstrated an increased apoptotic fraction in PTPRG-depleted cells ( Figure 9E–H ), indicating that PTPRG loss inhibits proliferation and promotes apoptosis in MM.