A total of 421 patients with BM tumors, whose data included somatic mutation, CNV information and clinicopathological details ( Supplementary Table S1 ), were enrolled in this study from the Memorial Sloan Kettering Cancer Center (MSKCC). These data were sourced from three previously published studies (18–20). The primary tumors of the BM patients included non-small cell lung cancer (NSCLC), melanoma (SKCM), breast cancer (BRCA), and other tumor types with relatively small proportions, amounting to a total of 21 types. Somatic mutation data of all the above-mentioned patients were obtained using the MSK-IMPACT targeted sequencing approach (21). In this study, the Oncotator software was employed to conduct standardized annotation of all the mutation data (22), aiming to eliminate sequencing biases and facilitate subsequent analyses.

This study retrospectively integrated pretreatment whole-exome sequencing (WES) data of somatic mutations, clinical characteristics, and immune checkpoint inhibitor (ICI) treatment follow-up information from 109 NSCLC (23, 24) and 631 SKCM (7, 25–31) patients. NSCLC patients were treated with anti-PD-1 monotherapy or combination immunotherapies. SKCM patients received anti-CTLA-4, anti-PD-1/PD-L1, or combination therapies. Somatic mutation profiles from pretreatment samples were also annotated using the Oncotator. For subsequent analyses, we focused on nonsynonymous mutations encompassing missense mutations, frameshift insertions/deletions (indels), inframe indels, nonsense mutations, and splice site mutations. Detailed clinical characteristics and ICI treatment regimens, including therapeutic responses, were systematically documented in Supplementary Table S2 (NSCLC patients) and Supplementary Table S3 (SKCM patients). The detailed workflow of this study was shown in Figure 1 .

The mutational process signature extraction workflow was performed using non-negative matrix factorization (NMF) algorithm (32). Somatic mutation data of BM patients, formatted as Mutation Annotation Files (MAF), were preprocessed with the R ​maftools package (33) to generate a 96-dimensional trinucleotide context matrix, where rows represented samples and columns corresponded to the 96 possible single-base substitution types within their trinucleotide contexts. Low-frequency mutations (occurring in < 1% of samples) were excluded to minimize noise.

The NMF algorithm, executed via the ​NMF package, decomposed the normalized mutation count matrix into two non-negative matrices: a signature matrix (W, k × 96 dimensions, defining mutational patterns) and an exposure matrix (H, sample × k dimensions, quantifying signature contributions per sample). The optimal number of signatures (k) was determined through cross-validation (10 iterations per k) using the ​nmfEstimateRank function, with stability and reconstruction error metrics. Extracted signatures were annotated by cosine similarity matching against the COSMIC v2.0 database (34) via the ​cosmicsig package to identify known etiological processes. In this study, all mutational signatures were defined as binary variables [following the method described in a previous study (35)]: if the activity of a certain mutational signature exceeded 25%, it was considered that this signature was present in the sample.

Molecular subtyping of BM patients was performed via consensus clustering analysis using the R package ​ConsensusClusterPlus (36). Mutational signature activity data, quantified as the normalized exposure of each signature per sample. Signatures with standard deviation < 0.1 (normalized exposure) were excluded to retain biologically informative patterns, as low-variance signatures may reflect technical noise rather than true molecular heterogeneity. The consensus clustering algorithm was applied with the following parameters: maximum cluster number (k= 8), 1000 iterations, resampling 80% of samples per iteration (pItem = 0.8), Pearson correlation distance, and hierarchical clustering with Ward’s linkage. Cluster stability was evaluated by consensus cumulative distribution function (CDF) curves and delta area analysis, with the optimal k determined by the plateau in consensus index scores and minimal inter-cluster cross-talk. Final molecular subtypes were defined based on the highest cluster consensus (> 0.9).

Somatic mutation co-occurrence and mutual exclusivity analyses were performed using the R package ​maftools (33). For BM patients with available MAF, recurrently mutated genes (altered in ≥ 5% of samples) were first identified, followed by filtering low-penetrance variants (germline polymorphisms excluded by dbSNP v155). Gene pairs exhibiting significant co-occurrence or mutual exclusivity were assessed using the somaticInteractions() function, which applies Fisher’s exact test (two-tailed) with Benjamini-Hochberg correction (FDR < 0.1) to evaluate deviation from random mutation patterns. Prognostic relevance was determined by integrating survival data into a Cox proportional hazards model, where gene pairs with significant interaction terms (P < 0.05) were prioritized. Heatmaps of co-occurrence/exclusivity networks were generated using the ​ComplexHeatmap package (37).

Statistical analyses and data visualization were performed using R software. Pie charts, box plots, and heatmaps were constructed with the ggplot2 package. Lollipop plots visualizing amino acid alteration patterns induced by gene mutations were generated using the maftools package. Survival curves were plotted via Kaplan-Meier analysis, and between-group differences were evaluated by the Log-rank test. Multivariable logistic regression and Cox regression models incorporating confounding factors were implemented with the forestmodel package. In the multivariate model, samples with missing values were deleted. We performed Schoenfeld residual analysis (using the survival R package) to validate that the Cox model proportional hazards (PH) assumptions were met. Wilcoxon rank-sum test and Kruskal-Wallis H test were applied to compare non-normally distributed variables across two and three groups, respectively. Two-sided P values less than 0.05 were deemed statistically significant.

A total of 421 tumor patients with BM and 5882 somatic mutations were included in this study. The median overall survival (OS) for the entire cohort was 21.4 (95% CI: 17.4–28.4) months ( Figure 2A ). The mutational profile of patients with BM is predominantly characterized by C > T single-nucleotide polymorphism missense mutations ( Supplementary Figure S1 ). Primary tumors of BM encompassed 21 histological subtypes ( Figure 2B ), with NSCLC (46.6%), SKCM (15.0%), and BRCA (9.7%) representing the three most prevalent subtypes, followed by colorectal cancer, esophagogastric cancer, and rare tumor types (e.g., uterine sarcoma, pancreatic cancer, and cervical cancer).

Kaplan-Meier analysis revealed that patients with BM originating from NSCLC, SKCM, or BRCA exhibited significantly prolonged survival compared to those with other primary tumors (Log-rank test, P = 0.004; Figure 2C , Supplementary Figure S2 ). This study further explored the prognostic roles of the indicators related to tumor genomic instability, namely fraction genome altered (FGA), microsatellite instability (MSI), and tumor mutational burden (TMB). The results demonstrated that a higher TMB was associated with a better prognosis (HR: 0.89, 95% CI: 0.82-0.97, P = 0.007, Figure 2D ). Meanwhile, we found that BM patients originating from NSCLC, SKCM, or BRCA harbored a significantly elevated TMB (Wilcoxon rank-sum test, P < 0.001; Figure 2E ). These findings underscore TMB and tumor origin as dual determinants of BM outcomes.

Based on the NMF algorithm applied to the somatic mutation matrix of BM patients, four distinct mutational process signatures were extracted ( Figure 3A ). Through cross-referencing with the COSMIC database of annotated mutational signatures, these extracted signatures were systematically classified as COSMIC signatures 4, 6, 7, and 13 according to standardized nomenclature ( Figure 3B ). Their characteristic mutational patterns were illustrated in Figure 3C and detailed mutational activities were shown in Supplementary Table S4 . Signature 4 has been implicated in tobacco exposure-related carcinogenesis, while signature 6 is associated with defective DNA damage repair mechanisms. Signature 7 was identified as the ultraviolet light-induced mutational pattern, and signature 13 was linked to APOBEC enzyme activity. Notably, the APOBEC-associated mutational signature demonstrated widespread prevalence across multiple tumor types. In this BM cohort, patients carrying this signature exhibited significantly elevated mutation burden ( Supplementary Figure S3A ) and tCw motif burden (a specific nucleotide context preferentially targeted by APOBEC enzymes) ( Supplementary Figure S3B ). Furthermore, the presence of APOBEC signature showed significant correlations with mutations in multiple cancer-related genes ( Supplementary Figure S3C ).

We subsequently investigated the prognostic relevance of the four identified mutational signatures. Survival analysis revealed no significant association between the presence of signatures 4, 7, or 13 and clinical outcomes in BM patients ( Supplementary Figures S4A–S4C ). Strikingly, patients harboring signature 6 (DNA damage repair associated) exhibited significantly worse overall survival (Log-rank test, P = 0.003; Figure 3D ). This association remained robust in a multivariable Cox regression model adjusted for age, sex, and tumor type (HR: 1.31, 95% CI: 1.01–1.70, P = 0.036, Figure 3E ). Moreover, when signature 6 was analyzed as a continuous variable in the Cox regression model, its correlation with the inferior prognosis retained statistical significance (HR: 1.39, 95% CI: 1.05–2.04, P = 0.041; Figure 3F ). Subsequently, we conducted a stratified analysis of the prognostic ability of signature 6 based on the primary tumor type. The results showed that in NSCLC/SKCM/BRCA, patients carrying signature 6 still exhibited poorer prognosis (Log-rank test, P = 0.047; Supplementary Figure S5A ). However, in other tumor types with a smaller proportion, there was no significant survival difference between patients with and without this mutational signature (Log-rank test, P = 0.238; Supplementary Figure S5B ), but a consistent trend was observed.

Intriguingly, BM patients carrying signature 6 demonstrated significantly reduced TMB compared to non-carriers (Wilcoxon rank-sum test, P < 0.001, Figure 3G ). Emerging evidence suggests that TMB serves as a surrogate marker for tumor immunogenicity (38, 39), which critically influences antitumor immune responses. Our findings propose a potential mechanistic link wherein signature 6 may potentially compromise prognosis through diminished immunogenicity, though further functional validation is warranted.

While most molecular classification studies in oncology rely on transcriptomic profiling, such approaches are inherently challenged by high-dimensional data complexity and technical noise. This study innovatively leveraged the activity profiles of extracted mutational signatures to investigate potential molecular subtypes in BM patients through consensus clustering analysis ( Supplementary Figure S6A ). The clustering optimization analysis revealed maximal consensus stability with three clusters, as evidenced by the inflection point (maximum slope) in the CDF curve ( Figure 4A ). Clustering consensus of three clusters was further validated by the plateauing consensus index in the clustering bar plot ( Figure 4B ), with detailed clustering heatmap visualizing in Supplementary Figure S6B .

Kaplan-Meier survival analysis demonstrated that cluster 3 exhibited significantly improved survival compared to clusters 1 and 2 (Log-rank test, P = 0.009, Figure 4C ). This prognostic advantage persisted in the multivariable Cox regression model adjusted for clinicopathological covariates (HR = 0.60, 95% CI: 0.41–0.90, P = 0.011, Figure 4D ). Subsequently, we conducted a stratified analysis of the prognostic ability of identified clusters based on the primary tumor type. The results showed that in NSCLC/SKCM/BRCA, patients from cluster 3 still exhibited favorable prognosis compared with cluster 1 and 2 (Log-rank test, P = 0.038, Supplementary Figure S7A ). However, in other tumor types with a smaller proportion, there was no significant survival difference between identified three clusters (Log-rank test, P = 0.382, Supplementary Figure S7B ), but a consistent trend was observed.

BM patients from the cluster 3 also displayed an elevated TMB relative to other clusters (Kruskal-Wallis H test, P < 0.001; Figure 4E ), suggesting enhanced tumor immunogenicity as a potential mechanism underlying its favorable prognosis. Furthermore, the four mutational signatures exhibited differential activity distribution across clusters. Cluster 3 was characterized by heightened signature 7 activity alongside suppressed activities of signatures 4, 6, and 13 ( Figure 4F ), indicating distinct biological etiologies among molecular subgroups. This signature activity distribution provides mechanistic insights into the observed survival disparity and TMB heterogeneity.

We systematically investigated the prognostic implications of all somatic mutations in BM patients. Initial profiling revealed high-frequency mutations (≥ 5% mutation rate) in key driver genes including TP53, KRAS, EGFR, PTPRT, PTPRD, APC, et.al., with TP53 (58% of patients) emerging as the predominant mutational driver in BM ( Supplementary Figure S8 ). Subsequently, Cox regression identified 12 recurrently mutated genes significantly associated with clinical outcomes (P < 0.05), namely PTPRT, ARID1A, FAT1, PAK7, EPHA7, SPEN, NOTCH3, ATM, TET2, PREX2, CARD11, and PIK3CG ( Figure 5A ). Unexpectedly, while RB1 mutations correlated with inferior prognosis ( Supplementary Figure S9 , Figure 5A ), all other prognostic mutations demonstrated protective effects ( Figure 5A ).

Functional enrichment analysis mapped these prognostic genes to critical oncogenic pathways including NOTCH signaling, HIPPO signaling, and cell cycle regulation ( Supplementary Figure S10 ), with their potential therapeutic targetability stratified by druggable categories in Supplementary Figure S11 . Furthermore, we explored the co-occurrence and mutual exclusivity among the prognostically mutated genes ( Figure 5B , Supplementary Table S5 ). The results revealed significant co-occurrence patterns between PTPRT and multiple prognostic genes, with the highest interaction frequency ( Figure 5C ), suggesting a central role for PTPRT mutations in modulating BM outcomes through synergistic genetic interactions.

Previous analyses have indicated that PTPRT played a crucial role in the mutation interaction network of BM. In both univariate survival analysis (Log-rank test, P = 0.002, Figure 6A ) and multivariable Cox regression model (HR = 0.54, 95% CI: 0.36–0.82, P = 0.003, Figure 6B ), we found that patients with PTPRT mutations had significantly improved survival outcomes. The lollipop plot presented the detailed amino acid changes caused by PTPRT mutations, with R364Q and E548K being the hotspot mutations ( Figure 6C ). Further analysis demonstrated that BM patients carrying PTPRT mutations harbored significantly elevated TMB levels (Wilcoxon rank-sum test, P < 0.001, Figure 6D ).

Interestingly, stratified survival analysis by tumor origin revealed preserved prognostic advantage of PTPRT mutations in NSCLC/SKCM/BRCA (Log-rank test, P = 0.012, Figure 6E ), but not in other subtypes. Considering that NSCLC and SKCM are two common tumor types for current immunotherapy, and the TMB was significantly increased in patients with PTPRT mutations, we hypothesized whether this mutation could evaluate the efficacy of immunotherapy. Based on a cohort of 109 NSCLC patients who received ICI treatment (ICI cohort 1), this study found that patients with PTPRT mutations exhibited better ICI prognosis (Log-rank test, P = 0.011, Figure 6F ). Meanwhile, the treatment response rate (the proportion of CR and PR) was significantly increased (mutant vs. wild-type: 58.3% vs. 32.2%, Fisher exact test, P = 0.038, Figure 6G ), and TMB was also enriched in the mutant patients (Wilcoxon rank-sum test, P < 0.001, Figure 6H ). In the validation cohort of 631 ICI-treated patients with SKCM (ICI cohort 2), improved survival outcomes ( Figure 6I ), increased response rates ( Figure 6J ), and elevated TMB ( Figure 6K ) were also observed in PTPRT mutant patients. PTPRT mutations improved survival and enhanced immunotherapy response in NSCLC and SKCM patients, which suggests that PTPRT may serve as a pan-cancer biomarker for immune checkpoint inhibitor efficacy.

We systematically investigated the associations between copy number variations (CNVs) of all included genes and clinical outcomes in BM patients. Two critical genes, PTEN and DUSP4, were identified through this screening, both exhibiting recurrent copy number deletion patterns. Survival analysis revealed that BM patients harboring PTEN deletions demonstrated significantly reduced overall survival (Log-rank test, P = 0.008, Figure 7A ). This association retained statistical significance in multivariable Cox regression models after adjusting for key prognostic confounders including age, sex, and primary tumor histology (HR: 1.71, 95% CI: 1.05–2.78, P = 0.032, Figure 7B ).

Conversely, DUSP4 deletions conferred a protective effect, with carriers showing significantly prolonged survival compared to non-deletion controls (Log-rank test, P = 0.006, Figure 7C ). The favorable prognostic value of DUSP4 deletions remained robust upon multivariable adjustment (HR: 0.31, 95% CI: 0.13–0.76, P = 0.011, Figure 7D ).