Identification of an angiogenesis-related risk score model for survival prediction and immunosubtype screening in multiple myeloma

Data preprocessing

Gene expression profiles, somatic mutation data, and associated clinicopathological data of MM patients were retrieved from the Multiple Myeloma Research Foundation (MMRF) program (https://research.themmrf.org and http://www.themmrf.org). GSE24080 and GSE5900 from the GEO repository were utilized to acquire clinical parameters and normalized gene expression data [16, 17]. Patients lacking significant clinicopathological or survival information were excluded from further analyses. Thirty-six AAGs were obtained from the MSigDB Team (Hallmark Gene set) (Supplementary Table 1). The maftools R package was used to analyze and visualize the masked somatic mutation data [18].

Identification of expression patterns of AAGs

According to the expression of 35 AAGs (for which VTN was excluded because it was undetected), an unsupervised clustering analysis was carried out to identify distinct angiogenesis-related patterns. The number and robustness of the clusters were assessed by a consensus clustering algorithm [19]. The R package ConsensusClusterPlus executed the above steps in 1000 iterations to ensure the robustness of the classification [20]. The gene set variation analysis (GSVA) algorithm was applied to the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set (c2.cp.kegg.v7.4) to determine the biological functional differences in the AAG-related clusters [21].

Immunohistochemistry (IHC)

After written informed consent was obtained, IHC targeting of human CD34 and VEGFA proteins was performed on 4-μm-thick bone marrow specimens from 20 MM patients hospitalized at the Affiliated Hospital of Shandong University of Traditional Chinese Medicine from July 2020 to July 2022 following informed and consent-approved protocols from the ethics committee of the Affiliated Hospital of Shandong University of Traditional Chinese Medicine (2020) (ethical review No. (010) - KY). After dewaxing and rehydration, the sections were incubated with 0.3% H2O2 formaldehyde to inhibit endogenous peroxidase activity and subjected to microwave antigen retrieval and cooled at room temperature. After the cells were incubated with 10% goat serum, primary antibodies, including mouse anti-human CD34 (1:200; Servicebio, China) and mouse anti-human VEGFA (1:200; Servicebio, China), were applied. After washing, the sections were incubated with a biotinylated anti-mouse secondary antibody and labeled with streptavidin-peroxidase solution (Servicebio). The presence of antigen was visualized by staining sections with DAB (Servicebio) and counterstaining sections with hematoxylin (Servicebio).

Evaluation of VEGFA expression and microvessel density (MVD)

The determination of the plasma cell percentage in VEGFA-stained bone marrow sections was performed manually. The calculation formula for the H score was as follows: H score = intensity of staining × % positivity. Megakaryocytes strongly express VEGF, so they served as internal positive controls [22].

The degree of angiogenesis was assessed using CD34-labeled microvessels. The regions with the highest number of microvessels were identified at 100× magnification and then converted to 400× magnification to measure the microvessel numbers in each region. The average was obtained as the MVD for each region [23].

Associations between molecular patterns and the clinical characteristics and prognosis of MM patients

To determine the clinical significance of the clusters by consensus clustering, we investigated the associations between molecular patterns and clinical features and between molecular patterns and survival outcomes. The clinical variables included age, gender, type, MRI, LDH, ALB, and HGB. Moreover, the differences in event-free survival (EFS) and overall survival (OS) between patients with different patterns were evaluated via Kaplan–Meier (K–M) analyses via the survival and survminer R packages [24].

Correlations between molecular patterns and immune characteristics

We used single-sample gene set enrichment analysis (ssGSEA) to assess the immune microenvironment based on the immune gene sets (Supplementary Table 2) obtained from the research of Charoentong [25], which can provide a better picture of tumor conditions than CIBERSORT when tumor cells are present in similar proportions. The expression of immune checkpoints (ICPs), including PD-L1, CTLA4, PD-1, PD-L2, LAG3, and HAVCR2, was subsequently compared between the different groups. Additionally, the correlations between the AAG_score and the differentially expressed ICPs were calculated.

For MM, most related studies have focused on CD138-positive cells obtained via magnetic cell sorting (MACS). Due to the limitations of the sorting methods, the remaining nontumor cells were randomly mixed in the extracted samples, still characterizing the matrix environment of MM in a relatively large sample quantity [26].

Identification of immune-related differentially expressed genes (DEGs) and functional annotation

The limma R package was used to determine the DEGs between distinct angiogenesis clusters with the cutoff criteria |logFC| ≥ 1 and p value < 0.05 [27]. To explore the potential biological processes related to the DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the clusterProfiler R package [28]. The immune-related DEGs were determined by the intersection of the immune gene set (Supplementary Table 2) and the DEGs.

Generation of the angiogenesis-associated prognostic AAG_score and nomogram

We performed K–M analysis for immune-related DEGs, and 26 genes associated with both OS and EFS were selected with the criterion of p < 0.01. Based on the combined role of these 26 genes in MM progression, the AAG_Score signature was constructed to comprehensively assess the role of these molecules in patient prognosis, the immune microenvironment, and immunotherapy response. Least absolute shrinkage and selection operator (LASSO) Cox regression was used to reduce the dimensionality, and the 11 most stable molecules were selected to construct the AAG_score signature with the glmnet R package [29]. The formula is as follows:

AAG_score = (0.2874 × expression of IFI16) + (−0.1035 × expression of STAP1) + (0.1546 × expression of GEMIN6) + (−0.2298 × expression of SLC7A7) + (−0.0694 × expression of LST1) + (−0.0011 × expression of IGHM) + (−0.0710 × expression of FUCA1) + (−0.0009 × expression of PD-L2) + (0.1405 × expression of NUF2) + (−0.1112 × expression of CD22) + (−0.0192 × expression of ADAM28).

The AAG_score model was validated in the test set and the entire cohort. Cox regression analyses were subsequently performed to determine the independent prognostic factors. The risk prediction model was subsequently constructed as a nomogram according to the independent prognostic factors. Calibration curve analysis and decision curve analysis (DCA) were performed to evaluate the performance of the constructed nomogram.

Acquisition of immunotherapeutic cohorts

The IMvigor210 immunotherapeutic cohort with RNA-seq data and complete clinical data was included in our study [30]. The IMvigor210 study investigated the efficacy of an anti-PD-L1 antibody in advanced or metastatic urothelial carcinoma patients. The complete transcriptome profile and clinical information were downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies/, and the count was converted to transcripts per kilobase million (TPM) values. The DEseq2 R package was used for normalization.

Statistical analysis

Comparisons between two groups were performed using the Wilcoxon test. One-way ANOVA and the Kruskal–Wallis test were used to analyze differences among three or more groups. Spearman and distance correlation analyses were used for correlation analyses. All the statistical analyses were performed using R software (version 4.1.0) and its relevant packages. P < 0.05 was regarded as statistically significant.

Differential expression and genetic mutation landscape of AAGs in MM

The detailed flowchart of this work is shown in Figure 1. We first evaluated the expression levels of the 35 AAGs in MM samples and normal samples in the GSE24080 and GSE5900 datasets. A total of 17 DEGs were found (Figure 2A). A protein–protein interaction (PPI) network was constructed through Mode to reveal the hub genes (Supplementary Figure 1A); among these genes, only VEGFA and ITGAV were significantly upregulated in MM patients. Immunohistochemistry was used to examine the relationship between MVD and the expression of VEGFA, a key molecule involved in angiogenesis, and survival (Figure 2B). The Surv_cutpoint function was used to determine the optimal cutoff points for VEGFA (0.7993967) and MVD (1.5243902) in the datasets, dividing the patients into high- and low-expression groups. K–M analyses revealed that patients in the groups with a lower MVD and lower expression of VEGFA had better survival (Figure 2C, 2D).

Figure 1. The entire analytical process of the study.

Figure 2. Differential expression of AAGs in MM and generation of angiogenesis subgroups. (A) Differences in the expression of AAGs between MM patients and healthy controls (^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001). (B) Representative images of MM patient BM tissues stained with CD34 and VEGFA antibodies (scale bar, 50 μm). Overall survival analyses of the MVD group (C) and VEGFA group (D) using Kaplan–Meier curves. (E) Consensus matrix heatmap depicting three gene clusters associated with the AAGs. The colors indicate the similarity or dissimilarity among samples at each iteration, with blue indicating higher similarity and white indicating lower similarity. (F) PCA based on the AAGs revealed three disjointed populations. Red, cluster A samples; green, cluster B samples; blue, cluster C samples. (G) OS analysis of patients in the AAG_clusters using Kaplan–Meier curves.

Next, we selected bone marrow samples from primary MM patients from the MMRF to evaluate the degree of variation. The findings revealed that missense mutations, SNPs, and C>T mutations were more common, with the highest mutation frequency being 1883 (Supplementary Figure 1B). Furthermore, we explored somatic mutations in AAGs in MM. As the waterfall diagram depicted (Supplementary Figure 1C), only 78 of the 724 (10.77%) MM samples presented genetic mutations, and the findings suggested that VCAN and COL3A1 are the genes with the highest mutation rates among the AAGs. These findings indicated the potential role of AAGs, especially the hub genes VEGFA and ITGAV, in MM pathogenesis, despite the low mutation intensity of these AAGs.

Generation of angiogenesis subgroups in MM

A total of 559 MM patients from the GSE24080 dataset were enrolled in this study to determine the relationship between angiogenesis and clinical characteristics. To explore the relationship between the expression patterns of AAGs and MM subtypes, we performed a consensus clustering analysis to classify MM patients according to the expression levels of these AAGs. The optimal number of clusters was determined by evaluating the quality of clustering based on consensus values, the cumulative distribution function, and the proportion of ambiguous clustering. Our findings indicated that the optimal clustering variable was 3 (Figure 2E), and MM patients in the entire cohort were well dispersed in cluster A (n = 171), cluster B (n = 260), and cluster C (n = 128). Principal component analysis (PCA) also revealed three distinct populations (Figure 2F). The clinicopathological variables of these three clusters are shown in Supplementary Figure 2A. Furthermore, the EFS and OS times of patients in the AAG clusters were explored, and a significant difference in EFS was observed (Figure 2G); however, the results revealed no significant difference in OS. Additionally, GSVA was performed, as displayed in Supplementary Figure 2B; cluster C was enriched in hematopoiesis-associated pathways, metastasis-associated pathways, and immune-associated pathways, such as hematopoietic cell lineage, cell adhesion molecules (CAMs), ECM receptor interaction, cytokine receptor interaction, complement and coagulation cascade, allograft rejection, and the intestinal immune network. Therefore, these three AAG clusters exhibited substantial differences in clinical characteristics and biological functions.

Characteristics of immunity in different subgroups

To determine the relationship between AAGs and the immune characteristics of MM patients, we analyzed the associations between the three clusters and 28 immune metagenes. As shown in Figure 3A, a significant difference was found in the enrichment of the 28 immune metagenes among the three clusters, most of which exhibited high expression in cluster C cells, especially mast cells, myeloid-derived suppressor cells (MDSCs), macrophages, regulatory T cells, activated dendritic cells, central memory CD8 T cells, neutrophils, activated B cells, immature dendritic cells, T follicular helper cells, effector memory CD8 T cells, natural killer T cells, and eosinophils. Furthermore, we investigated the associations between the three clusters and the expression of the 6 ICPs. The expression levels of CTLA4 and HAVCR2 in cluster C were significantly greater than those in the other clusters, while cluster B presented lower expression levels of PD-L1 and LAG3 (Figure 3B).

Figure 3. Characteristics of immunity in different subgroups. (A) Differences in 28 infiltrating immune cells among the three AAG clusters (^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001). (B) Expression of 6 ICPs in the three subtypes.

Identification of gene subgroups based on DEGs

To investigate the biological reactions in the three angiogenesis subgroups, we obtained DEGs among the three AAG clusters with the “limma” package and conducted GO and KEGG analyses. GO enrichment analysis demonstrated that these genes were enriched mainly in immune-associated biological processes (Figure 4A). KEGG analysis revealed the enrichment of immune- and cancer-associated pathways, indicating that angiogenesis is a key factor in regulating immunity and tumorigenesis. Then, we separated the immune-related genes from the DEGs and utilized K–M analyses to screen genes associated with survival. Finally, 26 genes associated with both OS and EFS were selected with the criterion of p < 0.01. To study the specific regulatory mechanisms involved, the patients were divided into two different gene clusters according to 26 prognostic genes by a consensus clustering method. The clinical characteristics of the two groups are shown in Figure 4B. K–M analyses revealed that patients with the A subtype had shorter OS and EFS than patients with the B subtype did (Figure 4C, 4D). Additionally, a boxplot of angiogenesis-related genes is presented in Figure 4E; only VEGFA, SERPINA5, and MSX1 were expressed at higher levels in the A subtype than in the B subtype.

Figure 4. Identification of gene subgroups based on DEGs. (A) GO and KEGG pathway enrichment analyses of DEGs among the three AAG clusters. (B) Cluster diagram of clinical characteristics between the two gene clusters. Overall survival (C) and event-free survival (D) analyses of patients stratified according to gene cluster analysis using Kaplan–Meier curves. (E) Differences in the expression of AAGs between the two gene clusters (^*P < 0.05; ^**P < 0.01; ^***P < 0.001).

Construction and validation of the prognostic AAG_score model

A prognostic model combining angiogenesis and immune features was developed based on the 26 immune-related DEGs. A total of 559 MM patients were randomly assigned to a training set (n = 275) or a test set (n = 274). The LASSO Cox regression model was used to screen the 11 most representative immune-related molecules concerning patient prognosis (Figure 5A, 5B). The results of univariate Cox analysis of the selected molecules are shown in Supplementary Table 3. With a cutoff value of 0.82 determined by the survminer R package (Figure 5C), patients in the training set were classified into high- and low-risk groups. As the AAG_score increased, patient mortality increased notably. A heatmap was constructed to show the expression of 11 representative genes in the two AAG_score groups (Figure 5D). Additionally, K–M analysis indicated that low-risk patients had a survival advantage over high-risk patients (Figure 5E), and the areas under the receiver operating characteristic (ROC) curves (AUCs) of 1-, 3-, and 5-year OS were 0.74, 0.75, and 0.76, respectively, indicating the perfect performance of the AAG_score in the training set (Figure 5F).

Figure 5. Construction and validation of the prognostic AAG_score model. (A) LASSO coefficient profiles of 26 DEGs related to immunity and prognosis. (B) Tenfold cross-validation for tuning parameter selection in the LASSO model. According to log(λ), the partial likelihood deviation graph was drawn, where λ is the tuning parameter. Partial likelihood deviation values are displayed, and the error bars indicate SEs. Dotted vertical lines are drawn at the optimal values according to the minimum criteria and 1-SE criteria. (C) The optimal cutoff point for dichotomizing patients into low and high AAG_score groups were determined by the survminer R package. The optimal cutoff value was 0.82. (D) Percentage of deaths in the high- and low-risk groups as the AAG_score increased. Expression patterns of 11 selected molecules in different risk groups. (E) Overall survival analysis of risk groups using Kaplan–Meier curves. (F) ROC analysis of the ability of the AAG_score to predict 1-, 3-, and 5-year OS and its specificity. (G) Alluvial diagram showing the changes in the AAG clusters, gene clusters, AAG_scores, and clinical outcomes. (H) Differences in the AAG_scores among the three AAG clusters. (I) Differences in the AAG_scores between the two gene subgroups.

We then investigated the distribution of all patients in three AAG clusters, two gene clusters, and two AAG score groups (Figure 5G). The AAG_score was highest in AAG_cluster B and lowest in AAG_cluster C (Figure 5H). The performance of the gene clusters was similar to that of the AAG clusters (Figure 5I), indicating that the higher the AAG_score was, the worse the survival was. Based on the above results, we concluded that the AAG_score was a good predictor of prognosis in MM patients.

The AAG_score was further tested in the test set and the entire set depending on the cutoff of the training set, and similar results were obtained. K–M survival analyses revealed that, compared with patients in the high-risk group, patients in the low-risk group had superior survival. The AUC still indicated good performance (Supplementary Figures 3 and 4).

Correlations between the AAG_score and clinical features

Univariate and multivariate Cox analyses, which included the clinical factors of age, gender, type, MRI, LDH, ALB, and HGB, demonstrated that the AAG_score could serve as an independent prognostic factor for assessing MM patient outcomes (Supplementary Table 4 and Figure 6A). To develop a clinically relevant quantitative method for predicting the mortality rate of MM patients, we constructed a nomogram that incorporated the AAG_score and independent clinical prognostic factors (Figure 6B). The calibration curves for 1-, 3-, and 5-year OS for the derived nomogram and the ideal model demonstrated that the predictions were highly accurate (Figure 6C). Furthermore, the DCA curves showing the clinical utility of each model indicated that the nomogram had the greatest net benefit for survival prediction, followed by the AAG_score, while other variables were less effective (Figure 6D–6F).

Figure 6. Correlations between the AAG_score and clinical features. (A) Forest plot showing that the AAG_score was an independent prognostic biomarker according to multivariate analyses of the entire cohort. (B) Nomogram for predicting the 1-, 3-, and 5-year OS of MM patients. (C) The 1-, 3-, and 5-year calibration plots of the nomogram. (D–F) The 1-, 3-, and 5-year decision curve analyses of the clinical benefit rate.

Associations of the AAG_score with the immune microenvironment and the immunotherapeutic response

We then investigated the correlation between 28 immune metagenes and the AAG_score. Activated CD4+ T cells, activated CD8+ T cells, memory B cells, effector memory CD4+ T cells, and Gramma delta T cells exhibited markedly increased infiltration in the high-risk groups; CD56 dim natural killer cells and central memory CD4+ T cells showed no significant difference between the two groups, while the other cell types presented the opposite trend (Figure 7A). We also analyzed the correlation between the AAG_score and each immune cell type. In addition to the correlation between the AAG_score and effector memory CD4 T cells, there was a significant negative correlation between the AAG_score and immune cell abundance (Figure 7B). The expression status of 6 ICPs was explored, and the majority of the patients in the low-risk group were highly expressed (Figure 7C).

Figure 7. Association of the AAG_score with the immune microenvironment and immunotherapeutic response. (A) Differences in the levels of 28 infiltrating immune cells between the low-risk and high-risk groups (^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001). (B) The correlation between the AAG_score and the infiltration of 28 immune cell types. Blue, negative correlation; orange, positive correlation. (C) Expression of 6 ICPs in the low-risk and high-risk groups. (D) Correlations between the 11 selected molecules and 28 infiltrating immune cell types. Red, positive; blue, negative. (E) The correlation between the 11 selected molecules and 6 ICPs. (F) OS analyses of the high- and low-AAG_score groups in the anti-PD-L1 immunotherapy cohort using Kaplan–Meier curves. (G) Differences in the AAG scores of patients in different anti-PD-L1 clinical response groups. (H) The proportion of patients in the high- and low-AAG_score groups who responded to anti-PD-L1 immunotherapy.

In addition, considering the significant role of the AAG_score in immune cell infiltration, we correlated 11 molecules with the abundance of immune infiltrating cells. Spearman correlation analyses revealed that LST1, IGHM, PD-L2, CD22, FUCA1, SLC7A7, and ADAM28 exhibited prominent positive correlations with most of the immune cell infiltrates, while NUF2, IFI16, and GEMIN6 presented the opposite trend (Figure 7D). We also found that the expression of 11 key molecules was negatively correlated with that of almost all immune checkpoint molecules, and GEMIN6, LST1, STAP1, and CD22 showed particularly significant correlations (Figure 7E).

Immunotherapy is considered the standard of care for the treatment of MM. Using data from the IMvigor210 cohort receiving anti-PD-L1 immunotherapy, we investigated the ability of the AAG_score to predict the clinical response of patients to immune checkpoint blockade therapy. As depicted in Figure 7F, patients in the low-risk subgroup had a better prognosis than patients in the high-risk subgroup did. Patients with progressive disease exhibited a greater risk than did those with a partial response (Figure 7G). We also observed an enhanced therapeutic response to PD-L1 blockade therapy in low-risk patients compared to high-risk patients, implying that the low-risk group may be more sensitive to immunotherapy (Figure 7H).