Identification of the green-yellow module involved in CD8 + T cells and CD274 expression in LUAD BrMs via WGCNA

The CD8+ T cell level and PD-L1 expression are essential factors in determining the tumor microenvironment classification and the responsiveness of ICIs in NSCLC [37]. Thus, we performed weighted gene co-expression network analysis (WGCNA) by integrating the infiltration level of CD8+ T cells, cytotoxic lymphocytes, and CD274 expression. In our study, the top-5000 variation genes in RNA-seq data and power of β = 4 (scale free R2 = 0.89) as the soft threshold were selected to build the co-expression network. The analysis indicated that 13 co-expression modules were identified by applying the dynamic tree cut method (Fig. 2A, Figure S1A-F). Furthermore, we explored the relationships between the module gene and molecular traits, including CD8+ T cell infiltrations, cytotoxicity lymphocytes infiltrations, and CD274 expression. The heatmap of module trait correlations showed that the green-yellow module was highly positively correlated with CD8+ T cells and CD274 levels in LUAD BrMs, which meant that the green-yellow module (149 genes) may play an essential role in estimating the tumor immune microenvironment classification and the clinical response of ICIs (Fig. 2B). The scatter plots in Fig. 2C, D showed similar results. Thus, the green-yellow module was recognized as the key module.

Fig. 2

Module-trait relationships and green-yellow-related enrichment analysis. A Cluster dendrogram of DEGs based on different metrics. Each color depicts a module that possesses weighted co-expressed genes. B Correlation of module eigengenes to CD8+ T cells and PD-L1 in LUAD BrMs. The values in the cells of each column are presented as p-value. Red represents a positive correlation, and blue represents a negative correlation. *: P < 0.05, **: P < 0.01, ***: P < 0.001, ****: P < 0.0001. C, D Scatter plot show the correlation between Gene significance for CD8+ T cells or CD274 and module eigengenes in greenyellow module. E, F Gene Ontology and KEGG enrichment analysis for 149 genes in greenyellow module

Furthermore, the GO and KEGG analysis demonstrated that the genes of the green-yellow module were enriched in various immune processes and pathways. GO analysis indicated that lymphocyte activation, regulation of T cell activation, T cell activation and positive regulation of leukocyte activation were enriched in the green-yellow gene set (Fig. 2E). KEGG analysis showed that hematopoietic cell lineage, natural killer cell mediated cytotoxicity, cytokine-cytokine receptor interaction, T cell receptor signaling pathway were enriched in green-yellow module (Fig. 2F). These pathways had been demonstrated to be linked to immune process and immunotherapy [37, 38]. This shown the importance of the green-yellow module genes for the immune microenvironment and immunotherapy in LUAD BrMs, and it is required to further explore its characteristics and roles.

Construction and immune prediction of GYMS

Considering the individual differences in the immune microenvironment and fully exploiting the predictive immune value of green-yellow module genes, we constructed an immune scoring system (green-yellow module score, GYMS) by applying the GSVA algorithm based on 149 genes of green-yellow module in Xiangya BrMs cohort (n = 70). These patients were divided into low-GYMS (n = 35) and high-GYMS (n = 35) group based on the median value 0.0051 as a constant threshold across all cohort. To explore the performance of GYMS to predict the outcome in LUAD BrMs patients, we conducted a survival analysis for patients with survival information (n = 55). The survival analysis indicated that the high-GYMS group (n = 29) had a better OS benefit than low-GYMS group (n = 26) in our Xiangya BrMs cohort (P = 0.025, Fig. 3A). Moreover, the result of survival analysis had been verified in an immunotherapy cohort in lung cancer from Samsung Medical Center (high-GYMS vs low-GYMS; P = 0.06, Figure S2A) and the TCGA-LUAD cohort (high-GYMS vs low-GYMS; P = 0.031, Figure S3A). Furthermore, we investigated the correlation between the GYMS and clinical molecular information. Our analysis revealed no significant difference in the proportion of sex, age, smoking status, and molecular mutations between the high- and low-GYMS groups (Figure S4). Nevertheless, we observed a higher proportion of patients with KRAS mutations in the high-GYMS group (high-GYMS vs low-GYMS, 23.5% vs 6.1%), although statistical significance was not reached (Figure S4).

Fig. 3

Immune characteristics of the green-yellow module score (GYMS). A Kaplan–Meier curve shows the effect of GYMS on OS with the constant threshold of 0.0051. B Heatmap shows the spearman’s correlation between GYMS and co-inhibitors. C Heatmap shows the expression in immunomodulators (MHC I, MHC II, co-stimulators, chemokines, and receptors) between low- and high-GYMS. D Correlation between GYMS and immune cells. E Correlation between GYMS and the infiltration levels of seven types of immune cells (CD8+ T cells, NK cells, CTL, macrophages, dendritic cells, MDSC, and Treg), which were estimated using six independent algorithms. F Heatmap shows the spearman's correlation between GYMS and ESTIMATE score. G Comparison of ESTIMATE score, immune score, and stromal score between low- and high-GYMS

Apart from the prognostic value, the GSEA was performed to explore the transcriptomic features produced by high-GYMS BrMs. By using GO, KEGG, and Hallmark gene sets as references, GSEA was performed between the high-GYMS and low-GYMS groups in patients from the Xiangya cohort. The level of GYMS in LUAD BrMs patients was positively correlated with the immune-related biological process and pathways, such as adaptive immune response, antigen processing and presentation, chemokine signaling pathway, natural killer cell mediated cytotoxicity, IFNG response, and inflammatory response (Figure S5; Table S5). Therefore, high-GYMS BrMs were correlated with high antigen presentation capacity and immune cell infiltrations.

Then, we evaluated the correlation between immunological characteristics and GYMS in LUAD BrMs. The GYMS positively correlated with the multiple immunomodulators, including MHC, Co-stimulators, and chemokines, which was consistent with the GSEA results (Fig. 3C). The high-GYMS group had higher expression of MHC I/II molecules than the low-GYMS group, showing a solid antigen presentation capacity (Figure S6D). A majority of chemokines and acceptor, including CXCL9, CXCL10 and CXCR3, were significantly upregulated in the high-GYMS group (P < 0.01; Figure S6A), promoting the recruitment of effector immune cells. Generally, the high expression of immune checkpoint molecules such as PD-L1/PD-1/CTLA4 was reported in inflamed TIME [39]. We found that the GYMS was positively correlated with many immune checkpoint molecules, including PD-L1, PD-1, CTLA-4, HAVCR2, LAG3, IDO1, PDCD1LG2, TIGIT, and BTLA (P < 0.01; Fig. 3B).

The immune score, stromal score and ESTIMATE score were also positively correlated with the GYMS (P < 0.001, Fig. 3F). High-GYMS had higher immune score, stromal score, and ESTIMATE score than the low-GYMS (P < 0.0001, Fig. 3G). Next, the correlation between GYMS and immune cells infiltration was evaluated using six different algorithms. We found that GYMS was positively correlated with overall level of immune cells (Fig. 3D, E). In addition to showing superior immune-activated cells (e.g. cytotoxic lymphocytes, CD8+ T cells, NK cells), high-GYMS BrMs also indicated a more increased abundance of immunosuppressive cells (e.g. M2-type TAMs, MDSC, Treg, and fibroblast). Meanwhile, the consistent results were observed in primary lung cancer in two external cohorts (the SMC immunotherapy cohort and the TCGA-LUAD cohort; Figure S2D-G, Figure S3B-E). Then, resident cells (astrocyte and microglia) in brain also were explored in this study. High-GYMS BrMs had higher infiltrations of the microglia and reactive astrocyte than low-GYMS BrMs (P < 0.01, Figure S7A,B). Therefore, these results demonstrated that high-GYMS BrMs were correlated with the inflamed tumor microenvironment that may well be susceptible to immunotherapy. However, it also indicates significant immune escape in high-GYMS BrMs.

The above results were validated in protein level. 66 FFPE samples from BrMs of LUAD was obtained to conduct immunohistochemical staining for immune cell markers (Table S3). High-GYMS BrMs had a higher infiltrating density of CD8+ T cells, CD68+ macrophages, and CD163+ M2 macrophages compared to low-GYMS BrMs (P < 0.05, Fig. 4B, D, F). Representative images are shown in Fig. 4A, C, E.

Fig. 4

Immunohistochemical analysis shows the immune infiltration characteristics. Comparison of infiltrating density of CD8+ T cells (A, B), CD68+ macrophages (C, D), and CD163+ M2 macrophages (E, F) between high- and low-GYMS group. All tumor slides were stained with antibodies against CD8 (Cat. ab209775, 1:2000, Abcam), CD68(ZM-0060 Ready-to-Use, ZSGB-BIO), and CD163(ZA-0428 Ready-to-Use, ZSGB-BIO). Scale bars: 100 μm in 200x and 50 μm in 400x

To further explore new immune-related targets, the Cytohubba functions of the Cytoscape were used to discover ten most essential hub genes in the green-yellow module, including CD8A, LAG3, TMEM176B, TNFRSF13B, SH2D2A, KLHL6, PARP15, CXCL10, IGLL5, and SLAMF6 (Figure S1G). The CD8A is the marker gene of CD8+ T cells and the LAG3 is the important immune checkpoint molecule, which further demonstrates the accuracy of the above results. Given the clinical utility of 10 genes would be easier than 149 genes, utilizing the expression profiles of these 10 genes, we developed a 10-gene scoring system using the GSVA method. Nonetheless, our analysis revealed no significant difference in prognosis between the high and low score groups, although there was a trend toward longer median survival in the high score group (P > 0.05; Figure S1H).

The potential role of GYMS in ICIs responsiveness

To investigate whether GYMS can predict the clinical response to ICIs, we calculated the TIDE score in our Xiangya BrMs cohort. We found that the GYMS was positively correlated with CD274 (r = 0.34; P = 0.0042), IFNG (r = 0.71; P < 0.0001), CD8 (r = 0.64; P < 0.0001), and Dysfunction (r = 0.59; P < 0.0001) and negatively correlated with MDSC (r = − 0.32; P = 0.0061), TAM M2 (r = − 0.54; P < 0.0001), and Exclusion (r = − 0.37; P = 0.0018) (Fig. 5B). Besides, a strong negatively correlation (r = − 0.37; P = 0.0015) between GYMS and TIDE score had been found in the Xiangya BrMs cohort (Fig. 5B). The high-GYMS group had a borderline significant lower TIDE score than the low-GYMS group (P = 0.076, Fig. 5D). It has been reported that the low-TIDE score is a surrogate biomarker to predict a good response to cancer immunotherapy [28]. As expected, higher proportion of responder was found in the high-GYMS group compared to low-GYMS group (Fig. 5C). Meanwhile, we also used the SubMap algorithm to evaluate the clinical response to ICIs (PD-1 and CTLA-4 inhibitors) in high- and low-GYMS patients with LUAD BrMs. In this study, the high-GYMS BrMs displayed more promise in response to anti-PD-1 inhibitors (Bonferroni-corrected P = 0.008, Fig. 5E). In addition, we found that GYMS was positively correlated with T-cell inflamed score (r = 0.85; P < 0.0001) that can predict the clinical response of ICIs (Fig. 5A). Meanwhile, the expression of several therapeutic immune checkpoints, including PD-L1, CTLA-4, TIM-3, and LAG-3, was significantly higher in the high-GYMS group (P < 0.05, Figure S6B). Consistent with above results, in the anti-PD-1/PD-L1 NSCLC cohort, the high-GYMS group had a higher proportion of responders with durable clinical benefits (DCB) than the low-GYMS group (Figure S2B,C). Therefore, these results indicated that the GYMS is a potential biomarker with applications in predicting immunotherapy response in both LUAD patients with and without BrMs.

Fig. 5

Immunotherapeutic responsiveness of the green-yellow module score (GYMS). A Scatter plot shows the correlation between GYMS and T-cell inflamed score. B Bar plot shows the correlation between GYMS and TIDE score. C Comparison of the proportion of the responder between low- and high-GYMS. D Comparison of the TIDE score between low- and high-GYMS. E Immunotherapeutic responses to anti-CTLA-4 and -PD-1 treatments in low- and high-GYMS patients

GYMS predicts immune phenotypes in external validation cohort of BrMs

In a NSCLC BrMs cohort (n = 43) from Sun Yat-Sen University Cancer Center, we observed that high-GYMS (n = 25) had higher expression of Co-stimulators and chemokine than low-GYMS (n = 18), such as CXCR3, CXCL9, CXCL10, CCL4, and CCL3 (Fig. 6C). The GYMS also was positively correlated with immune score (P < 0.001; Fig. 6A). Meanwhile, GYMS was uncovered to be positively correlated with T cells, CD8+ T cells, Cytotoxic lymphocytes, NK cells, Myeloid dendritic cells, monocytes, Treg, and MDSC (P < 0.05, Fig. 6D, E). GYMS was positively correlated with the T-cell inflamed score, showing high-GYMS BrMs had an inflamed tumor microenvironment (r = 0.58, P < 0.0001, Fig. 6F). Moreover, GYMS was negatively correlated with TIDE score; Reduced TIDE score was found in the high-GYMS group (P < 0.05, Fig. 6G). As expected, borderline significant higher proportion of responder was found in the high-GYMS group compared to low-GYMS group (P = 0.064; Fig. 6H). Then, the SubMap algorithm was used to evaluate the clinical response to ICIs, showing the high-GYMS BrMs might have better response to anti-PD-1 inhibitors (Bonferroni-corrected P = 0.024, Fig. 6I). Furthermore, GYMS was also positively correlated with expression of several immune checkpoints, including CD274, PDCD1, CTLA4, TIM3, LAG3, TIGIT (P < 0.05, Fig. 6B). Therefore, the GYMS can distinguish immune phenotypes and immunotherapy responsiveness in the external validation cohort of NSCLC BrMs.

Fig. 6

Immune characteristics and immunotherapeutic responsiveness of GYMS in OMIX575 cohort (43 BrMs). A Heatmap shows the spearman’s correlation between GYMS and ESTIMATE score. B Heatmap shows the spearman’s correlation between GYMS and co-inhibitors. C Heatmap shows the expression in immunomodulators (MHC I, MHC II, co-stimulators, chemokines, and receptors) between low- and high-GYMS. D Correlation between GYMS and immune cells. E Correlation between GYMS and the infiltration levels of seven types of immune cells (CD8 + T cells, NK cells, CTL, macrophages, dendritic cells, MDSC, and Treg), which were estimated using six independent algorithms. F Scatter plot shows the correlation between GYMS and T-cell inflamed score. G Bar plot shows the correlation between GYMS and TIDE score. H Comparison of the proportion of the responder between low- and high-GYMS. I Immunotherapeutic responses to anti-CTLA-4 and -PD-1 treatments in low- and high-GYMS patients

To broaden the application of GYMS, we explored the performance in Pan-cancer BrMs cohorts. In melanoma BrMs and breast cancer BrMs, the T-cell inflamed score was highly positively correlated with the GYMS (P < 0.0001, Figure S8A,E). Besides, the higher CD274 expression and lower TIDE score was found in the high-GYMS group compared to low-GYMS group (Figure S8 B,C,F,G). As expected, the high-GYMS BrMs may benefit more from anti-PD1 therapy in Pan-cancer BrMs by using SubMap algorithm (P < 0.01, Figure S8D,H).

Identification of potential therapeutic drugs for low-GYMS group

To improve the poor prognosis of low-GYMS BrMs patients, we employed two methods to identify potential drugs based on the drug response data from the CTRP and PRISM database. On the one hand, differential drug response between low-GYMS (bottom quintile) and high-GYMS (top quintile) groups was performed to determine drugs of lower estimated AUC values in the low-GYMS group with log2FoldChange > 0.10 and Wilcox rank sum test P < 0.05 (Fig. 7A). On the other hand, correlation analysis between the GYMS and estimated AUC values was utilized to search drugs with positive correlation coefficient with Spearman’s correlation coefficient > 0.25 (Fig. 7A). Then, 2 CTRP-derived compounds (ABT-737 and navitoclax) and 1 PRISM-derived compound (romidepsin) had lower estimated AUC values in low-GYMS group and positively correlated with GYMS (Fig. 7B,C). Therefore, the three compounds were identified as potential candidate drugs for low-GYMS BrMs patients.

Fig. 7

Identification of potential therapeutic drugs for low-GYMS BrMs patients. A The workflow for identifying drugs with higher drug sensitivity in low-GYMS BrMs patients. B Spearman’s correlation analysis and differential drug response analysis of 2 CTRP-derived compounds. C Spearman's correlation analysis and differential drug response analysis of 1 PRISM-derived compounds

Development and validation of a GYMS-related risk signature

To identify a biomarker to predict prognostic value, we constructed a GYMS-related risk signature in the Xiangya BrMs cohort. Differential expression analysis was used between the high-GYMS (n = 35) and low-GYMS (n = 35). A total of 354 differential expression genes (DEGs) with log2FC > 1 and adjust P < 0.05, including 283 up-regulated genes and 71 down-regulated genes (Fig. 8A; Table S6). Then, in 55 BrMs samples with survival information, the GYMS-related risk signature composed of 4 genes (MAB21L1, C14orf132, IL12RB1, and CCR2) was developed based on the 354 DEGs by applying Univariate Cox and LASSO Cox analysis (Fig. 8B, C; Table S7). Risk scores were calculated for each sample (risk score = 0.084* MAB21L1 + 0.029* C14orf132 + (− 0.005)* IL12RB1 + (− 0.077)*CCR2). BrMs patients of LUAD were divided into low-risk group (n = 28) and high-risk group (n = 27) using the median value as threshold (Fig. 8D). As shown in Fig. 8E, BrMs patients with high-risk had significantly shorter overall survival than those patients with low-risk (log-rank P < 0.0001), and consistent results were acquired in two external validation sets (the immunotherapy cohort and the TCGA-LUAD cohort; P < 0.05, Figure S9A, S10A). The time-dependent AUC of the risk signature was 0.66, 0.74, and 0.77 at 1-, 2-, and 3-year, respectively, which showed robust predictive value of the GYMS-related risk signature (Fig. 8F). Simultaneously, the AUC values at 1, 2, and 3 years for the GYMS-related risk signature outperformed those markers derived from immune score, CD8+ T cells, and CTL (Fig. 8F). The prognostic value of the risk signature was further validated in two external cohorts (Figure S9B,S10B). Then, the effect of the risk signature on the prognosis of BrMs patients was further validated. We found that the higher risk score correlated with poor prognosis in the univariate Cox regression analysis, which was statistical significance (P < 0.001; HR = 67.604, 95% CI 8.819–518.264, Fig. 8G). Next, three factors (Postoperative radiotherapy, Age, and GYMS-related risk signature) with statistical significance in univariate Cox analysis were imported into Multivariate Cox regression analysis (Fig. 8GH), and indicated that the GYMS-related risk signature was an independent prognostic biomarker in LUAD BrMs and primary LUAD by using the Xiangya cohort, the SMC immunotherapy cohort in NSCLC (Figure S9C,D), and the TCGA-LUAD cohort (Figure S10C,D) (P < 0.05). In addition to GYMS-related risk signature, Age and Postoperative radiotherapy are also independent prognostic factors in LUAD patients with BrMs (P < 0.05; Fig. 8H). Therefore, our findings confirmed that the GYMS-related risk signature could be used as a valuable prognostic indicator in LUAD patients with or without BrMs.

Fig. 8

Developing a GYMS-related prognostic signature using LASSO Cox regression. A Volcano plot of differently expression genes between the high- and low-GYMS group. The x-axis range is from − 6 to 10. B, C LASSO Cox analysis identified genes most correlated to overall survival in the Xiangya cohort. D Visualization of the association of the risk scores with survival status and gene expression profiles in LUAD BrMs. E Kaplan–Meier survival analysis of low- and high-risk group in LUAD BrMs. F Comparison of the 1-,2-,3-year AUC value between the GYMS-related signature and other similar score, including Stromal score, Immune score, ESTIMATE score, CD8+ T cell infiltration, and CTL infiltration in LUAD BrMs. G, H Forest plot of the GYMS-related risk signature in univariate and multivariate cox analysis. CTL: Cytotoxic T Lymphocyte. KPS: Karnofsky Performance Status