Prognostic value of integrated cytogenetic, somatic variation, and copy number variation analyses in Korean patients with newly diagnosed multiple myeloma

Background To investigate the prognostic value of gene variants and copy number variations (CNVs) in patients with newly diagnosed multiple myeloma (NDMM), an integrative genomic analysis was performed. Methods Sixty-seven patients with NDMM exhibiting more than 60% plasma cells in the bone marrow aspirate were enrolled in the study. Whole-exome sequencing was conducted on bone marrow nucleated cells. Mutation and CNV analyses were performed using the CNVkit and Nexus Copy Number software. In addition, karyotype and fluorescent in situ hybridization were utilized for the integrated analysis. Results Eighty-three driver gene mutations were detected in 63 patients with NDMM. The median number of mutations per patient was 2.0 (95% confidence interval [CI] = 2.0–3.0, range = 0–8). MAML2 and BHLHE41 mutations were associated with decreased survival. CNVs were detected in 56 patients (72.7%; 56/67). The median number of CNVs per patient was 6.0 (95% CI = 5.7–7.0; range = 0–16). Among the CNVs, 1q gain, 6p gain, 6q loss, 8p loss, and 13q loss were associated with decreased survival. Additionally, 1q gain and 6p gain were independent adverse prognostic factors. Increased numbers of CNVs and driver gene mutations were associated with poor clinical outcomes. Cluster analysis revealed that patients with the highest number of driver mutations along with 1q gain, 6p gain, and 13q loss exhibited the poorest prognosis. Conclusions In addition to the known prognostic factors, the integrated analysis of genetic variations and CNVs could contribute to prognostic stratification of patients with NDMM.


Introduction
Multiple myeloma (MM) is the second most common hematological malignancy in Korea and the USA [1][2][3]. The treatment outcomes of patients with MM have markedly improved owing to the development of novel therapeutic agents and the technical advances in molecular diagnostics and detection of minimal residual diseases [4,5]. There is a consensus on the risk factors that aid in stratifying patients with MM. The methods to detect genetic variations include conventional karyotyping and fluorescent in situ hybridization (FISH) [6,7].
Recently, the genetic landscape of MM has expanded owing to the use of high-throughput technologies, such as next-generation sequencing for analyzing targeted genes and wholeexome or genome [8][9][10][11]. Previous studies have revealed novel and frequently mutated genes, such as KRAS, NRAS, BRAF, and FAM46C in MM. Additionally, the chronological evolution of multiple driver events has been demonstrated using serial patient specimens [8]. However, the prognostic effect of most mutated genes, which have a low recurrence rate, has not been clearly identified [12]. There are limited studies on the comprehensive analysis of various cancer driver events and the correlation between structural variants and genomic events.
Recent studies have performed clustering analysis that integrates various prognosis-related genomic variations to classify patients with MM based on the characteristics of genomic alterations [8][9][10]. The heterogeneous nature of MM has warranted further studies on the multifaceted interpretation of subgroups incorporating various genetic variations and their prognostic relevance. Several studies have utilized emerging technologies to examine the prognostic implications of genetic variations. However, there are no consensus guidelines that include somatic variations and small structural variations.
In this study, the genomic profile of Korean patients with newly diagnosed MM (NDMM) was examined using whole-exome sequencing (WES) and new prognosis subgroups were identified through integrated analysis of copy number variations (CNVs) and somatic variants.

Patients
In this study, 67 patients with NDMM exhibiting more than 60% plasma cells in the bone marrow (BM) aspiration were recruited at the Seoul National University Hospital between July 2004 and January 2017. Patients were diagnosed to NDMM based on BM aspiration and biopsy according to the international myeloma working group (IMWG) 2016 guideline [13]. A total of three laboratory hematologists reviewed the slide for diagnosis and differential counting of plasma cells. Patients with >60% plasma cells were included in this study to maximize representativeness of variants for genetic changes of plasma cells. In particular, since the revised IWMG guideline designated plasma cell count over 60% as biomarkers of active myeloma, 60% was established as the selection criteria. The clinical characteristics, including the age of disease onset, sex, CRAB (hypercalcemia, renal impairment, anemia, and bone disease) symptoms, chemotherapy regimens, and survival, as well as the laboratory findings, including complete blood count, blood urea nitrogen/creatinine, albumin, lactate dehydrogenase levels, BM histological findings, and cytogenetic findings (FISH and conventional cytogenetics [CG]), of each patient were recorded. This study was approved by the Institutional Review Board of Seoul National University Hospital (IRB No. 1312-102-544). All study subjects provided their written informed consent to participate in the study.

DNA extraction and exome sequencing
DNA isolated from the frozen BM mononuclear cells was used in the exome capture protocol. The SureSelectHuman All Exon V5+UTR probe set included 359,555 exons of 21,522 genes, and the size of the total targeted region was 75 Mb. To generate the standard exome capture libraries, the Agilent SureSelect Target Enrichment protocol was used for generating the Illumina paired-end sequencing library (ver. B.3, June 2015) with 3 μg input genomic DNA (gDNA). The quantity and quality of DNA were examined using PicoGreen (Molecular Probes, Eugene, OR) and Nanodrop (ThermoFisher Scientific, Mississauga, ON). DNA concentration � 50 ng/uL, purity � 1.7 (A260/A280), volumn � 20 uL and total amount � 1 ug was passed for quality control criteria, and determined to acceptable for evaluation. The gDNA (1 μg) was fragmented using adaptive focused acoustic technology (AFA; Covaris). The fragmented DNA was repaired with an adenine ligated to the 3 0 end, followed by ligation of the Agilent adapters. The adapter-ligated product was subjected to polymerase chain reaction (PCR). The final purified product was then quantified using quantitative real-time PCR (qRT-PCR) following the qPCR Quantification Protocol Guide and subjected to quality control using the Caliper LabChipHigh Sensitivity DNA (PerkinElmer). For exome capture, 250 ng of all exon capture libraries were mixed with hybridization buffers, blocking mixes, RNase block, and 5 μL of SureSelect, following the standard Agilent SureSelect Target Enrichment protocol. Hybridization to the capture baits was performed at 65˚C using a heated thermal cycler with the lid temperature maintained at 105˚C for 24 h in a PCR machine. The captured DNA was amplified and the final purified product was quantified using qRT-PCR according to the qPCR Quantification Protocol Guide. The amplified product was subjected to quality control using the TapeStationDNAscreentape (Agilent). The pooled DNA libraries were sequenced using the NovaSeq6000platform (Illumina, San Diego, USA), conducted by outsourced service from Macrogen Inc. (Seoul, South Korea).

Bioinformatic evaluation of sequencing data
FASTQ files were aligned to the reference human genome (hg19; GRC37) using the Burrows-Wheeler aligner (BWA, v0.62) [14]. Duplicate PCR reads were removed using Picard 1.98. Variant calling was performed using the "HaplotypeCaller" in Genome Analysis Toolkit 2.7-2 [15]. To detect the candidate gene mutations, a filtering strategy was used. Low-quality variants with a low total depth (<20) and a low altered allele count (<10) were filtered out. Synonymous and noncoding variants were discarded. The variants with an allele frequency of more than 0.01 when compared with those in the 1000 Genomes Project (http://browser. 1000genomes.org/), the Exome Aggregation Consortium (http://exac.broadinstitute.org/), and the NHLBI exome sequencing project (ESP6500, http://evs.gs.washington.edu/EVS/) databases were excluded. As a matched control sample was not included in this study, a stringent variant selection pipeline was applied to prioritize the high-confidence set of somatic mutations. The driver genes were selected based on two previous studies [8,10]. Genes that were identified as driver genes in at least one of these two studies were selected as driver genes for this study (S1 Table in S1 File). Additionally, the in silico prediction algorithms, SIFT [16,17], CADD [18] and PolyPhen2 [19] were used, and the clinical significance was interpreted according to the ACMG guidelines [20].

Copy number analysis
Copy number alterations were analyzed using a CNVkit [21] and Nexus software version 5 (Biodiscovery, El Segundo, CA). The copy number in the NDMM dataset was called against an MM-negative karyotype FISH panel. Heatmap plots were drawn with the "heatmap" command in CNVkit. Data were loaded into Nexus 5.0 and the copy number calls were generated genome-wide for each sample based on the fixed thresholds for deletions and duplications specified in the settings. SNP-Rank Segmentation Algorithm [22], a statistics-based algorithm similar to the circular binary segmentation, applies both the copy number value and the Ballele frequency to the segmentation (S1 Fig in S1 File). GISTIC algorithms were used to determine the significance of focal somatic copy number alterations.

Statistical analysis
All statistical analyses were performed using the statistical R-project program (version 3.6.2), PASW statistics version 18 (SPSS Inc., Chicago, USA), and MedCalc version 12.0 (MedCalc Software, Mariakerke, Belgium). A pairwise correlation analysis between somatic mutations and CNVs was performed using the Fisher test. Cumulative overall survival (OS) curves for the groups with or without genomic variations were calculated using the Kaplan-Meier (KM) method and compared using the log-rank test and Breslow test. The Cox proportional hazards model was used to evaluate the prognostic impact of CNVs and mutated genes on OS. A multivariate analysis was performed on the full set of significant variables in the univariate analysis. The differences were considered significant at P < 0.05.

Characteristics of the study population
The demographic and clinical characteristics of the study population, including the laboratory tests performed on the day of MM diagnosis, are described in Table 1. The numbers of male and female patients were 34 and 33, respectively. The median age of the study population was 65 years. Among the total nucleated cells in the BM aspiration, the median percentage of BM plasma cells was 75.8%. The numbers of patients with ISS stage I, II, and III tumors were 7, 25, and 33, respectively ( Table 1).

Identification of CNVs in MM
The gain and loss of the p arm, q arm, and both arms of chromosomes 1 to 22 were analyzed. The most common chromosomal gain was 1q. Among the odd-numbered chromosomes, the most frequent whole-arm gain was observed in chromosomes 5, 7, 9, 15, and 19. The most common chromosome loss was 13q loss, followed by the losses of 16q, 22q, 8p, 14q, 1p, and 6q ( Table 2; Fig 1). The frequency of the gain and loss of each chromosome is represented in S2 Table in S1 File.

Analysis of selected driver gene mutations
Driver gene mutations were detected in 59 of the 83 genes selected based on two previous studies [8,10]. Variations classified as pathogenic, likely pathogenic, and variants of unknown significance (VUS) were separately selected based on the ACMG guideline to analyze the frequency and mutation types of each gene (Fig 2). The most frequently mutated driver gene was IGLL5, which was detected in eight patients. Additionally, seven patients had mutations in ATM, six patients had mutations in NRAS, KRAS, and DIS3, five patients had mutations in MAN2C1 and BHLHE41, four patients had mutations in MAML2, DUSP2, and BRAF, and three patients had mutations in TRAF2, TP53, TET2, and KDM6A. The diagram including the variation and domain for each gene is shown in S2 Fig in S1 File. The p.V600E variation in BRAF was detected in three patients, while the p.Q61R(L) variation in KRAS was detected in three patients. In NRAS, the p.G12D(V) and p.Q61L(K) variations were detected in two patients. Mutations, such as p.S403L, p.L581V, and p.Q584X, were detected in ATM (S2 Fig in S1 File).

Correlation analysis of genomic variants
CNVs were frequently detected with moderate or high correlation between each other. In patients with hyperdiploid MM, CNVs in odd-numbered chromosomes were highly correlated with each other. The CNVs between the following chromosome pairs that occurred simultaneously were highly correlated: chromosomes 5 and 9 (Pearson correlation coefficient (PCC) = 0.745; P < 0.001), chromosomes 5 and 15 (PCC = 0.791; P < 0.001), chromosomes 5 and 19    gain were significant independent prognostic markers (Table 3). Additionally, KM analysis revealed that patients with 1q gain (P = 0.03; Fig 4A), 6p gain (P = 0.003; Fig 4B), 6q loss (P = 0.009; Fig 4C), 8p loss (P = 0.03; Fig 4D), and 13q loss (P = 0.046; Fig 4E) exhibited poor prognosis. In patients with driver gene mutations, univariate and multivariate analyses revealed that mutations in MAML2 (HR = 3.32; 95% CI = 1.34-8.27; P = 0.010) and BHLHE41 (HR = 5.16; 95% CI = 1.91-13.9; P = 0.001) were significantly associated with poor OS. The KM plot of patients with driver gene mutations is shown in Fig 5. Patients with BHLHE41 (P < 0.001, Fig 5A) and MAML2 (P = 0.016, Fig 5C) mutations exhibited a significantly poor prognosis in the log-rank test, while those with BRAF mutations exhibited a significantly poor prognosis in the Breslow test (P = 0.036, Fig 5B). Additionally, CREBBP (P < 0.001), HIST1H1D (P = 0.044), PRDM1 (P = 0.025), KLHL6 (P < 0.001), TRAF3 (P = 0.049), and UBR5 (P = 0.003) mutations, which were detected only in one or two patients, were significantly associated with poor prognosis. Prognostic impact of CNVs identified through FISH, CG, and WES analyses. Quantitative results of FISH and CG are shown in Table 4. Compared with the CG analysis of metaphase chromosomes, FISH and WES analyses revealed a higher frequency of abnormalities. CNV and FISH analyses identified 1q gain in 52.2% and 60.0% of the patients, respectively. Moreover, CNV and FISH analyses identified 13q deletion in 49.3% and 48.1% of the patients, respectively. Survival analysis of patients with or without copy number alterations identified through FISH, CG, and WES analyses was performed (Fig 6). KM analysis revealed that 1q gain identified through FISH was not significantly correlated with OS (P = 0.831; Fig 6A). Similarly, 1q gain identified through FISH and CG analyses (FISH+CG) was not significantly associated with OS (P = 0.174; Fig 6B). Meanwhile, patients with 1q gain identified through WES Ideogram with regions of chromosomal gain and loss. Copy number gains and deletions are marked as blue and red bars to the right and left side of the ideogram, respectively.
Mutational burden (CNVs and driver gene mutations) as a prognostic factor. The distribution of the number of mutated driver genes in each patient is shown in Fig 7A. One or more mutations were detected in 63 patients (94%). The median number of total mutations was 2.0 (95% CI = 2.0-3.0; range = 0-8). The distribution plot of only VUS or higher mutations (ACMG classification) is shown in Fig 7B. The median number of mutations was 1.0 (95% CI = 1.0-2.0; range = 0-8). The CNVs were detected in 56 patients (83.6%). The median total CNV number in each patient was 6.0 (95% CI = 5.7-7.0; range 0-16; Fig 7C).
KM analysis revealed that both CNVs and driver gene mutations were associated with poor prognosis in the group with a high mutational burden. Patients with � 4 gene mutations exhibited lower OS than those with < 4 gene mutations (P = 0.035; Fig 7D). Additionally,    patients with � 3 mutant (VUS) driver genes exhibited poorer prognosis than patients with < 3 mutant driver genes (P = 0.033; Breslow test; Fig 7E). The survival rate of patients with � 4 CNVs was lower than that of patients with < 4 CNVs (P = 0.035; Fig 7F).

Factors affecting OS in patients with hyperdiploid NDMM.
Among patients with hyperdiploid NDMM, patients with � 5 trisomies exhibited better prognosis than those with < 5 trisomies (P = 0.037; Fig 8A). Patients with < 4 driver gene mutations exhibited a more favorable prognosis than those with � 4 driver gene mutations (P = 0.004; Fig 8B). The survival analysis of patients with hyperdiploid NDMM with or without driver gene mutations revealed that patients without BRAF mutations exhibited a better prognosis than those with BRAF mutations (P < 0.001; Fig 8C).
OS of clusters classified based on K-means analysis. Several factors affecting the survival rate of patients with NDMM identified in this study were classified according to the K-means clustering method. The following factors were included in the analysis: True-hyperdiploid (T-HRD) (>5 trisomy), 1q gain, 6q loss, 6p gain, 8p loss, 13q loss, 17 loss, t(4:14), t (14:16), and the total number of driver gene mutations for each patient. Upon classification into five clusters, 1q gain, 6p gain, 8p loss, 13q loss, and the number of driver gene mutations were identified as the significant clustering factors. Among them, the cluster 3 group with 1q gain, 6p gain, 13q loss, and driver gene mutation number of 7.5 as the clustering center was associated with poor prognosis. In contrast, the cluster 4 group with no CNVs and a low centering value with a driver gene mutation number of 1.8 exhibited the best prognosis. The cluster 5 group with driver gene mutation number of 5.6, 1q gain, and 13q loss exhibited significantly poorer OS than the cluster 4 group (P = 0.033; Fig 9).

Discussion
In this study, the CNV and somatic variant profiles of Korean patients with NDMM were identified using WES. A comprehensive analysis revealed that the gain of 1q and 6p chromosome arms and the loss of 6q, 8q, and 13q chromosome arms were associated with poor OS. Additionally, MAML2 and BHLHE41 mutations were shown to have an adverse prognostic impact on OS. Patients with a high frequency of CNVs or a high number of mutations exhibited poor prognosis. Cluster analysis revealed that patients with the highest number of total driver gene mutations along with 1q gain, 6p gain, and 13q loss were associated with the poorest prognosis. These findings suggested the utility of integrated analysis of CNVs and somatic mutations in predicting prognosis. The differences in the genomic profile between different ethnic groups were examined. The genomic profile of Korean patients with NDMM was compared with that reported previously [23][24][25][26]. Among the CNVs, 1q gain (52.2%) and 13q loss (49.3%) were most frequently detected, followed by losses of 8p, 14q, 1p, and 6q. The prevalence of 6p gain and 1q gain (26-45%) in this study was higher than that reported in previous studies. The frequency of 1q gain detection using FISH was high in Korean patients with myelodysplastic syndrome [27]. This indicated that 1q gain is a candidate CNV that is specific for Korean patients with myeloma. Among the somatic variants, the frequency of mutations in the following genes was higher than 6%: IGLL5, ATM, NRAS, KRAS, DIS3, MAN2C1, BHLHE41, MAML2, and BRAF. The frequency of IGLL5 mutation was high (18%), which is reported to occur exclusively with KRAS/NRAS mutation [28]. In contrast, the frequency of mutated genes in the MAPK pathway, such as KRAS, NRAS, and BRAF, was lower than that reported in previous studies on the Caucasian population (20-36%) [9,25,29,30]. This low frequency might be owing to the characteristics of the enrolled patients. Previously, Cifola et al. reported that the frequency of KRAS/NRAS mutations was low in plasma cell leukemia, which is the aggressive and high-risk form of plasma cell dyscrasia [31]. As this study performed integrative analysis on patients with � 60% plasma cells in the BM aspirate, the characteristics of patients with advanced disease could contribute to the low frequency of KRAS/NRAS mutations in this study. Meanwhile, the incidence of other mutated genes in patients with NDMM varied in different studies. This variability might result from the heterogeneous composition of patients with myeloma rather than the ethnic difference.
CNVs are considered to be one of the most important drivers of cancer development and progression [12,24,25,30,[32][33][34][35]. Recently, various methods have been developed for analyzing CNV through WES, which has enabled the integrated analysis of CNVs and somatic mutations [36] in addition to overcoming the limitations of CNV detection through CG and/or FISH. In the case of NDMM, the prognostic relevance was not clearly revealed except for the well-known variants such as t(4;14) and del (17). The effects of these abnormalities, especially for 1q gain and 13q loss on prognosis are controversial, and only a few studies have included them for risk stratification criteria [37][38][39]. The prognostic discrepancy could be caused by method sensitivity in addition to several unmeasuring confounding factors. A previous study showed that 5-53% of patients who were normal in FISH and cytogenetics had abnormality in chromosomal microarray test [40]. CNV detection through WES can overcome the limitations associated with low proliferative activity of cells or a low number of neoplastic plasma cells in the BM when compared with CNV detection through CG. Moreover, the detection of CNVs through single nucleotide polymorphism array or WES, which can encompass the entire chromosome, would be more helpful in predicting prognosis than FISH, which only detects small target chromosomal areas. In this study, 1q gain and 13q loss identified through WES analysis adversely affected OS, whereas 1q gain and 13q loss identified through FISH and/or CG methods did not have prognostic relevance. The deletion of whole arm of 13q (N = 29; 87.9%) and the gain of whole arm of 1q (N = 30; 85.7%) were observed in a majority of patients with 13q deletion and 1q gain in this study. Previously, Binder et al. analyzed the prognosis of patients with monosomy 13 and del(13q) separately in 1,181 NDMM patients, and reported that only patients with monosomy 13 had adverse prognosis [41]. Both of these studies suggest that the 13q loss, which had previously shown neutral or favorable results, may have included many patients with interstitial deletion rather than whole chromosome loss. Furthermore, it provides additional rationale for enabling the integration of WES and FISH tests to provide enhanced genetic information for risk stratification and improved prediction of outcomes.
Additionally, WES could explore the regions that are not covered by routine FISH probes. Recently, CNVs including 6p gain, 6q loss, and 8p loss have been reported to emerging recurrent alterations in multiple myeloma [40]. As CNV detection technologies have been developed, interests in the clinical implications of these new and emerging marker have raised. In this study, patients with 6p gain, 6q loss, and 8p loss were statistically significantly associated with poor prognosis. Further studies are needed to validate these CNVs as adverse prognostic markers in a large cohort.
In this study, BHLHE41 and MAML2 somatic mutations were associated with poor prognosis. BHLHE41, which is located at 12p12.1, is known as a helix-loop-helix superfamily domain that is involved in various cellular functions, such as proliferation, differentiation, tumorigenesis, and circadian rhythms [42]. The expression of BHLHE41 is reported to be upregulated in patients with Waldenstrom macroglobulinemia [43], and might be associated with poor prognosis in NDMM. MAML2 is located at 11q21. The prognostic relevance of MECT1-MAML2 and CRTC1-MAML2 fusion oncogenes has been reported in mucoepidermoid carcinomas [44,45]. In this study, MAML2 mutations exhibited positive correlation with various CNVs including 17q loss and 22q loss, and these variations may aid in predicting the prognosis. TP53 mutations are a well-known adverse prognostic factor in myeloma. Niccolo et al. demonstrated that TP53 mutations were associated with adverse progression-free survival (PFS) and OS and that a rare mutation in DNAH11 affected OS [9]. Another study reported that TP53, ATR, ATM, and ZFHX4 mutations are associated with poor PFS or OS [30]. Although TP53 did not show statistical significance in OS in this study (P = 0.613), it is presumed to be due to too small number of mutation-positive patients (N = 3). However, it could be owing to the differences in inclusion criteria, such as selecting only patients with >60% PC or targeting the Korean population. Most previous studies report the heterogeneous mutational landscape of MM. Some patients exhibited redundancy in gene mutations as two or more mutations were detected in genes involved in the same pathway [12,46,47]. Consistent with these findings the current study suggested the heterogeneity of genomic variants in MM. We suggest that the prediction power of mutational burden in structural and somatic variants is higher than that of a single variation.
In this study, we presented novel prognostic subgroups in patients with NDMM through K-means clustering analysis. The clusters were divided according to CNVs and mutations that affect OS. The "Cluster 3" exhibited the highest number of driver gene mutations and CNVs, including 1q gain, 13q loss, and 6q loss. Although efforts to stratify potential subgroups according to OS were attempted previously [9,10], the characteristics of each cluster group have not been distinguished. However, this study revealed that prognosis-related CNVs, such as 1q gain, 13q loss, and 6q loss have clustering relevance, whereas somatic gene mutations have no clustering importance. The overall number of mutations in the driver gene panel played a significant role in stratifying each group.
Recently, prognosis-related classifications are being reconsidered for patients with hyperdiploid MM [48,49]. In this study, patients with hyperdiploid MM and more than five trisomies exhibited favorable OS, which was based on a review of hyperdiploid variations that affected prognosis. Moreover, patients with hyperdiploid MM and an increased number of driver gene mutations exhibited poor prognosis. These results suggest that the integrated analysis of mutation burden (both CNVs and somatic mutations) using WES could aid in precisely predicting the clinical outcomes of patients with MM.
This study has several limitations; CD138 purification and sorting steps have not been implemented. Instead, only patients with >60% plasma cells in the BM aspirate were selected for the study. Most patients showed high numbers of plasma cells in the BM. Additionally, there were no control samples to remove the germline background. To overcome this limitation, the gene variants observed in healthy Korean individuals (n = 2,000) were removed. And, we focused on the driver genes selected by multiple algorithms in several previous studies. Based on the ACMG guidelines, variants with VUS, pathogenic, and likely pathogenic were separately analyzed. Furthermore, CNVs under 2,500 kb were excluded in CNV analysis. However, this study demonstrated that large CNVs with size > 2,500 kb were poor prognostic factors. As CNV analysis cannot detect chromosomal translocations, the results of translocations, such as t(4;14) and t(14;16) were referenced for clustering integrated analysis.
In conclusion, this study comprehensively analyzed the somatic mutations and CNVs of patients with NDMM using WES. Additionally, this study proposed a new method for classifying patient groups with poor prognosis and predicting OS. MM is a heterogeneous disease comprising several subclones, and multiple driver gene mutations are detected in one patient. Therefore, an integrated analysis through the application of WES in the future would aid in predicting the prognosis in a clinical setting.