182 core immune escape encoding genes were extracted from a study by Keith A. Lawson. et al. IMMPORT database (https://www.immport.org/) provided 2483 immune encoding genes (21). Venn diagram was drawn by Hiplot mapping tool (https://hiplot.com.cn/) to identify two sets of intersection genes. Downloaded data from TCGA database (https://tcga-Data.Nci.nih.gov/TCGA/), GTEX database (https://www.Gtex-portal.org/home/) and UCSC database (http://xena.Ucsc.edu/), which provided us with clinical data from 1095 BC cases and 292 normal controls, and their matched transcriptome RNA second-generation sequencing data.

The “t-test” of “GraphPad Prism 8” were used for differential expression analysis of normal and tumor samples, and the “survival” R package for survival analysis. The cBioportal database (http://www.cbioportal.org/) was used for further differentially expressed proteins analysis (22). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of differential expression of protein-coding genes were performed using “clusterprofiler” package. Category, Count, and -log10pvalue were used by GO analysis to obtain important metabolic pathways.

First, the gene sequencing data and patient survival data were combined, and univariate Cox analysis was performed using the “surviving” and “forestplot” R packages, and critical genes associated with prognosis were further explored using the “glmnet” R package. Multivariate cox analysis of candidate genes was performed with the “surv” package, and it constructed a prognostic risk assessment model.

Then, risk prediction and risk score calculation were performed using a prediction model. Based on the risk scores, patients were divided into high- and low- risk groups.

The “Hiplot” plotting tool was used to show the distribution of survival situation of low- and high- risk patients. The “pheatmap” R package was used to observe the difference in expression levels between high- and low- risk groups. Kaplan-Meier Survival analysis was used to assess the survival rates of BC patients in both groups. The validity of candidate genes in predicting the prognosis of BC patients was measured by plotting ROC curves using the “pROC” R package, which based on survival status, survival time and candidate gene expression in BC patients. The TCIA database (https://tcia.at/home) was used to analysis the enrichment of immune cells in the high- and low- risk groups (23). The prognostic model was externally validated by GSE20685 and GES42568 datasets.

The “T test” of “GraphPad Prism 8” was used for differential expression analysis of normal and tumor samples, and “survival” package was used to obtain survive-related hub genes, and P<0.05 was regarded statistically significant. The nomogram was used to predict the effect of each differentially expressed protein coding genes on 1-, 3-, and 5- year overall survival. The nomogram is composed of central prognostic factors, whose point scale is assigned to each variable. We used a horizontal line to determine the points for each variable and calculated the total points for each patient by adding points for all variables by normalizing the distribution from 0 to 100. Then we established the performance of the calibration curve for the visual nomogram. We compared the predictions in the calibration curves with the observed results. The best prediction was when the slope was close to 1.

The “XY” of “GraphPad Prism 8” performed correlation analysis of central protein-encoding genes with immune checkpoints PD-1, PD-L1 and CTLA4 and showed the relationship between them using Hiplot’s chord charts. The Gene Set Cancer Analysis (GSCA) database (http://bioinfo.life.hust.edu.cn/GSCA/) was used to analyze the expression of central protein-encoding genes and associated immune infiltration (24). The TIMER database (https://cistrome.shinyapps.io/timer/) was used to capture the immune cells and the correlation of tumor purity and quantity of gene expression, and represented by the heat map of the Hiplot (25).

Used the CELLMINER database (https://discover.nci.nih.gov/cellminer/) download the processed data “Processed Data Set” (26), and the data was processed by the “readxl” package. There were some “NA” missing values in the drug sensitivity data, and the “impute. KNN ()” function was used to evaluate and complement the missing values. Pearson correlation coefficients between each gene expression and different drugs were calculated and the two groups of drugs with the maximum and minimum correlation for each protein-coding gene were found based on the correlation coefficients. Correlation analysis was performed using “XY” of “GraphPad Prism 8” for visualization and validation purposes. The GES97681 dataset was used to verify whether Vemurafenib affects the expression of B-cell lymphoma 2-associated protein A1 (BCL2A1) and to find the pathway of drug effects through the literature, integrating the genes and pathways in the KEGG data (https://www.genome.jp/kegg/) to discover the effects of drugs on genes.

Human BC lines (MCF-7, MDA-MB-231, SK-BR-3) were obtained from the School of Biomedical Sciences, Chinese University of Hong Kong. MCF-7 maintained in RPMI-1640 medium and MDA-MB-231 were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco; Thermo Fisher Scientific, Inc.). Both media contained 10% fetal bovine serum (FBS; Thermo Fisher Scientific, Inc.) and 100 μg/ml streptomycin/100 U/ml penicillin (Gibco; Thermo Fisher Scientific, Inc.). SK-BR-3 were cultured in DMEM complete medium (iCell-h189-001b). Human breast epithelial cells (MCF-10A) purchased from iCell Bioscience Inc, China. MCF-10A was cultured in special medium (iCell-h131-001b). Vemurafenib (MedChemExpress, HY-12057) was dissolved in dimethyl sulfoxide (DMSO, final concentration is 0.1%) to prepare required concentrations. All cell lines were kept at 37°C in a humidified incubator with 5% CO2.

For the MTT assay, MCF-7, SK-BR-3 and MDA-MB-231 cells were detached by trypsinization and counted in a haemocytometer. Cells were seeded in 96-well plates at a density of 4000 cells/well in 100 μl medium per well. Twenty-four hours after incubation in the CO2 incubator, adherent cells were treated with increasing concentrations of drugs: vemurafenib (0, 20, 40, and 60 μM) in fresh 1640 and DMEM medium. For measuring the concentration of vemurafenib that resulted in 50% control growth inhibition (IC50), at 48 hours following drug treatment, 10 μl MTT (3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) (Sigma-Aldrich) solution (5 mg/ml in PBS) was added to each well after 48h of drug treatment. In addition, cell viability was measured by adding MTT at 0, 24, 48, 72h. Then incubation was continued for additional 4h in the CO2 incubator. The blue formazan product, formed by reduction in live attached cells, was dissolved by adding 100 μl of 100% DMSO per well. The plates were gently swirled at room temperature for 10 minutes to dissolve the precipitate. Absorbance was monitored at 490 nm using a microplate reader (CYTATION 3; Agilent Technologies, Inc.). The effect of vemurafenib on cell viability was assessed as the percent of cell viability compared with vehicle-treated control cells, which were arbitrarily assigned 100% viability. Dose-response curves and IC50 values were obtained using GraphPad Prism software. At least 3 dose-response experiments were performed for each compound, and the mean IC50 ± SD was calculated.

Cells were seeded at the density described above and incubated for 24 hours before the addition of vemurafenib. The compounds were added to the final concentrations of 0, 20, 40, and 60 μM. After 48 hours of treatment, Cell pellets were collected for experiments. Total RNA was isolated from MCF-7, SK-BR-3 and MDA-MB-231 by using Trizol reagent (Thermo Fisher Scientific, Inc.) according to the manufacturer’s instructions. Total RNA (1 μg) from each group of treated cells was converted to cDNA using a FastKing RT reagent kit (Tiangen, Inc.). The RT reaction was performed at 42°C for 15 min and 95°C for 3 min. qPCR was performed with a SYBR Green Real Time PCR kit (Thermo Fisher Scientific, Inc.) on CFX96 Touch Real Time PCR System (BioRad Laboratories, Inc.) under the following conditions: 95°C for 1 min, then 40 cycles of 95°C for 5 sec and 60°C for 15 sec. The gene expression level was calculated as 2(^−△△Ct) method, and the △Ct means the Ct of target gene minus the Ct of reference gene. The primers used for real-time PCR were BCL2A1 (forward) 5’-AAATTGCCCCGGATGTGGAT-3’ and (reverse) 5’-ACAAAGCCATTTTCCCAGCCT-3’; GAPDH (forward) 5’-CTGGGCTACACTGAGCACC-3’ and (reverse) 5’-AAGTGGTCGTTGAGGGCAATG-3’.

Western blotting was used to test the PI3K/AKT signaling pathway and BCL2A1 protein. RIPA lysis buffer (Beijing Solarbio Science & Technology Co., Ltd.) containing protease inhibitor PhosSTOP EASYpack (Roche Diagnostics) was used for total protein extraction according to the manufacturer’s protocol. Protein concentration was measured by the BCA kit (Beyotime Biotechnology, China). Then, protein samples were separated by 10% SDS−PAGE gel and transformed into PVDF membranes (Millipore, USA). Afterwards, membranes were incubated using 5% bovine serum albumin (BSA, Sigma−Aldrich; Merck KGaA) at room temperature for 1h and incubated with primary antibodies under 4°C overnight. The antibodies are as follows: anti-PI3K (1: 1,000, 4249, CST), anti-AKT (1: 1,000, 9272S, CST), anti-phosphorylated (p)−AKT (1: 1,000, 9271S, CST), anti-BCL2A1 (1: 1,000, 14093, CST), and anti-GAPDH (1: 2,000, GTX100118, GENE TEX) with GAPDH being the endogenous control. Afterwards, membranes were incubated with HRP−conjugated secondary antibodies at room temperature for 1h using a secondary antibody (1: 3,000, A0208, Beyotime). Finally, ECL blotting detection reagents (Clarity; Bio−Rad Laboratories, Inc.) was utilized to observe protein blots and the signals were detected by the ChemiDoc™ Imaging System (Bio−Rad Laboratories, Inc.).

Databases we used have: IMMPORT database, TCGA database, GTEX database, UCSC database, cBioportal database, TCIA database, GEO database (include GSE20685, GES42568 and GES97681 datasets), GSCA database, TIMER database, CELLMINER database. The softwares we used have: Hiplot mapping tool and GraphPad Prism.

Our research mainly used “GraphPad Prism 8” and “R language” for data difference analysis, visualization, etc. From the “GraphPad Prism 8”, “T test” was used to compare the differences between two data sets, “one-way ANOVA” was used to compare differences between three or more groups of data, and “XY” was used for correlation analysis. P<0.05 was considered statistically significant, indicating significant differences between the data were used. R is the language and operating environment for statistical analysis and mapping. we used the “Pheatmap”, “Survival”, “GGPUBR” and “clusterprofiler “ package fort expression analysis, survival analysis and enrichment analysis for differentially expressed genes and proteins. All experiment data are presented as the mean ± standard deviation (SD), further analyzing using GraphPad Prism 9.0. Differences in the results of two groups were evaluated using either two-tailed Student’s t test or one-way ANOVA followed by post-hoc Dunnett’s test. The differences with P< 0.05 were considered statistically significant.

Venn plot was used to identify the overlap of 2483 human immune associated genes from Immport database with 182 immune escape associated genes from Keith A. Lawson’s study ( Figure 1A ), yielded 31 crossover genes ( Figure 1B ). We compared the expression levels of immune escape encoding genes between BC patients and healthy individuals, found that 18 genes were significantly up-regulated and 12 genes were significantly down-regulated in BC patients, and ERAP1 was not significantly different in BC patients and healthy individuals ( Figure 1C ). BC patients were divided into high and low groups according to “res.cut” expression levels, and Kaplan-Meier survival analysis was performed for BC immune escape gene expression levels. Among them, 15 genes showed significant differences in survival and expression levels in BC patients ( Figure 1D ), and the survival curves of the other 16 genes were also shown ( Figure S1A ).

We identified 181 differentially expressed proteins regulated by BC immune escape encoding genes through the cBioPortal database ( Table S1 ). The volcanic plots ( Figure 2A ) show proteins with significant changes (P<0.05) and were used for further analysis. To determine which pathway differential proteins were mainly enriched, we performed GO and KEGG enrichment analysis using “clusterprofiler” R package, which showed that these differentially expressed proteins play important roles in protein serine/threonine kinase activity, ubiquitin-like protein ligase binding, membrane raft, cell substrate connection, gland development andreproductive structure development. The PI3K-Akt signaling pathway is considered to be the most prominent downstream signaling pathway for immune escape related genes in BC ( Figures 2B, C ).

Univariate Cox analysis was performed on 181 differentially expressed protein coding genes associated with immune escape in BC, and 31 of them were identified to be associated with survival ( Figure 3A ). To evaluate differentially expressed protein-coding genes meaningfully associated with BC prognosis and to obtain a better-fitting model, we used LASSO to downscale the high-dimensional information by adding a constraint to the absolute value of the coefficients and further screening by 10-fold cross-validation. The optimal λ value (λmin=0.025) was obtained from the minimum local likelihood deviation, and 16 genes were found to be significantly associated with prognosis ( Figures 3B, C ). The previous one looks at each feature individually to see if it is related to survival, while multivariate cox regression is to test whether multiple features are related to survival at the same time. Multivariate Cox analysis recognized 11 excellent prognostic gene models, namely EIF4EBP1, BCL2A1, NDRG1, MRE11A, ERRFI1, CLDN7, PIK3CA, CASP9, BRD4, PDCD4, and G6PD ( Figure 3D ). The risk score of patients was reckoned based on the scoring formula: Risk score = (EIF4EBP1 expression level * 0.162425485) - (BCL2A1 expression level * 0.234522698) + (NDRG1 expression level * 0.158199675) + (MRE11A expression level * 0.3105604) - (ERRFI1 expression level * 0.268609023)+ (CLDN7 expression level * 0.173133629) + (PIK3CA expression level * 0.442654646) - (CASP9 expression level * 0.337535821) - (BRD4 expression level * 0.497664344) - (PDCD4 expression level * 0.140816004) + (G6PD expression level * 0.190558316). There were more deaths in the high-risk group than the low-risk group, while fewer patients survived longer than 5 years than in the low-risk group ( Figure 4A ). The heatmap shows the risk of prognostic genes with a multivariate outcome P value less than 0.05. As the risk score increased, the expression levels of MRE11A, PIK3CA, EIF4EBP1 and NDRG1 increased, while the expression levels of BCL2A1 and BRD4 decreased ( Figure 4C ). Kaplan-Meier survival curve assessed the difference in survival between high- and low-risk patients in the risk model, which showed that survival decreased more significantly over time in the high-risk group than in the low-risk group, and the mean prognosis was worse in the high-risk group than in the low-risk group (P<0.0001) ( Figure 4E ). The prognostic value of Kaplan-Meier survival curves was determined by the P-ROC curve, and the area under the curve (AUC) was 0.672, indicating the median reliability of the kaplan-Meier survival curve ( Figure 4G ). The AUC of the survival assessment model was 0.68 for 2 years, 0.71 for 5 years, 0.75 for 8 years, and 0.77 for 10 years ( Figure 4I ). The applicability of the prognostic model was externally validated using the GSE20685 dataset ( Figures 4B, D, F, H, J ) and GES42568 dataset ( Figure S2 ), and the results further confirmed the reliability of the risk model. We imported the high- and low- risk groups into the TCIA database and found that the high-risk groups had a lower percentage of immune cells clustering ( Figure S2 ).

Seven of the 11 prognostic indicators were statistically significant, when Kaplan-Meier survival analysis was performed on them according to expression level, six of them were significantly different ( Figure 5A ). Analysis of their expression levels in tumor patients and normal subjects identified five core prognostic markers, which were EIF4EBP1, BCL2A1, NDRG1, ERRFI1 and BRD4 ( Figure 5B ). We plotted the expression levels of five central prognostic factors based on 1069 TCGA-BC patients. By testing these parameters to obtain the corresponding assigned scores for their expression, the individual gene scores were summed to find the total score, which can predict the survival of BC patients for 1, 3, and 5 years ( Figure 5C ). A calibration curve was plotted to verify the accuracy of the prediction of the line graph. In the calibration diagram, the predicted result (blue line) was very close to the actual result (black line), with a c-index of 0.695, indicating that the prediction quality of Nomogram was very high ( Figure 5D ). In conclusion, five core prognostic markers can accurately predict the prognosis of BC patients.

Correlation analysis of 5 core prognostic factors and immune checkpoints (PD-1, PD-L1 and CTLA4) showed that BCL2A1 had the highest correlation with immune checkpoints, while other factors with no significant difference ( Figure 5E ). Correlation analysis of the 5 central protein-encoding genes with 26 immune cells showed that BCL2A1 had the highest correlation with immune cells ( Figure 5F ). The correlation of 5 central protein-coding genes with tumor purity and 6 kinds of immune cells (B cells, CD4+ T cells, CD8+ T cells, Neutrphils, Macrophages and Dendritic cells) was further analyzed. BCL2A1 was negatively correlated with tumor purity and positively correlated with immune cells ( Figure 5G ).

We analyzed the expression levels of five core prognostic factors in two sets of BC classification data from public databases TCGA ( Figure 6A ) and GSE96058 ( Figure 6B ). The expression levels of EIF4EBP1, BCL2A1, NDRG1, ERRFI1 and BRD4 in all BC patients are consistent. The expression levels of EIF4EBP1 and BCL2A1 in BC patients are significantly higher than those in normal breast tissues, and with the increase of tumor malignancy. Its expression level increased.

Analysis of the potential correlation between the expression of 5 core protein coding genes and drug sensitivity in different human cancer cell lines from CellMiner database, which showed that the expression of BCL2A1 was positively correlated with the drug sensitivity of Vemurafenib ( Figure 7A ). The results of the correlation between other drugs and gene expression are placed in the Supplementary Table ( Table S2 ). We used RNA sequencing data from 78 BC patients in GSE97681, 24 without and 54 with Vemurafenib, the expression levels of BCL2A1 were analyzed and were found to be significantly reduced with the drug ( Figure 7B ). A review of the literature and the KEGG database revealed a map of the regulatory pathways of Vemurafenib on BCL2A1-sensitive targets (19) ( Figure 7C ).

We first detected the expression levels of PI3K, AKT, p-Akt and BCL2A1 in MCF7, MDA-MB-231, SKBR3 and MCF10A cells by Western Blot analysis. It was confirmed that the expression of these proteins in BC cells (MCF7, MDA-MB-231, SKBR3) was higher than that in normal breast cells (MCF10A)( Figure 7D ).

To explore the effect of vemurafenib on BC cells, the concentration of vemurafenib that resulted in 50% control growth inhibition (IC50) of MCF7, MDA-MB-231 and SKBR3 cells was assessed by MTT ( Figure 7E ). The IC50 of MCF-7 is 42 μM, MDA-MB-231 is 34 μM and SKBR3 is 51 μM. Next, the inhibitory activity of different concentrations of vemurafenib (0, 20, 40, 60 μM) on MCF7, SKBR3 and MDA-MB-231 cells at 0, 24, 48, and 72 h were further studied by MTT ( Figure 7F ). The results showed that vemurafenib inhibited cell proliferation in a dose-dependent manner, and the inhibitory effect increased with time.

To further explore the mechanism by which vemurafenib inhibit the growth of MDA-MB-231, MCF7 and SKBR3 cells, the expression level of PI3K/AKT signaling pathway and BCL2A1 in vemurafenib treated cells was evaluated by western blot ( Figures 8A–C ). The expression level of PI3K, AKT, p-AKT and BCL2A1 were reduced in vemurafenib treated cells compared to without vemurafenib. This suggested that the drug inhibits the PI3K/AKT signaling pathway and BCL2A1. In addition, apoptosis-associated genes expression was detected. BCL2A1 downregulation was confirmed using RT-qPCR ( Figures 8D–F ), and the results indicated that vemurafenib promoted cell apoptosis of MDA-MB-231 and MCF-7.