We obtained somatic mutation data of breast invasive carcinoma (BRCA) samples from the TCGA GDC database by choosing the “Masked Somatic Mutation” (https://portal.gdc.cancer.gov/) (Zhang et al., 2021). The preprocessing was employed with VarScan software and the somatic mutations were visualized using the MAFtools R package (Mayakonda et al., 2018). Then, we downloaded the RNA sequencing data (FPKM values) of the BC patients and subsequently converted FPKM values to TPM values. Moreover, the data were divided into the lncRNA and mRNA expression profiles. Further, we download the clinicopathologic features and outcomes in the same population from the UCSC Xena platform (http://xena.ucsc.edu/), such as sex, age, stage, and microsatellite instability (MSI) status (Speir et al., 2016). In addition, two datasets including somatic mutation and clinical data in BC patients were downloaded from the International Cancer Genome Consortium (ICGC) database (https://daco.icgc.org/), which were Breast Cancer-FR (BRCA-FR) and Breast Cancer-KR (BRCA-KR) (Zhang et al., 2019).

To analyze the copy number variations (CNVs) of TP53 in TCGA-BRCA patients, we obtained the data of Masked Copy Number Segment using the TCGAbiolinks package in R language (Colaprico et al., 2016). The CNV data was processed using GISTIC 2.0 by performing the GenePattern5 function (Reich et al., 2006). During the analytical process, GISTIC 2.0 with default settings was used except for several parameter (i.e., the confidence was 0.99 and X chromosome was included before the analysis). Finally, the results of GISTIC 2.0 were visualized with the MAFtools R package.

To predict the response to ICI therapy caused by the TP53 mutation in BC patients, we computed the TMB, MSI, and tumor immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu) for each BC sample. The total number of the somatic mutations per megabase of the genome detected in the tumor was defined as the TMB (Yarchoan et al., 2017); the insertion or deletion of repeat units results in a change in the microsatellite length, which is referred to as MSI (Vilar and Gruber, 2010); TIDE is a computational framework that can evaluate the response to immunotherapy and predict tumor immune escape by analyzing the gene expression profiles of cancer cases (Jiang et al., 2018). We calculated all TMB, MSI, and TIDE values for each sample, and compared their differences between patients with wild-type TP53 and those with mutant TP53 using a Wilcoxon rank-sum test.

To explore the significance of mutant TP53 in BC progression, we classified the TCGA patients into “TP53-mutant” and “TP53-wt” groups. The holistic analysis was employed by principal component analysis (PCA), which is a multivariate statistical technique under the broad title of factor analysis, that focus on pattern recognition and signal processing (Ringnér, 2008). PCA was conducted with the R packages factoextra and FactoMineR. DEGs were determined using the Bioconductor R package DESeq2 (Love et al., 2014), and the threshold for DEGs was p < 0.01 and |logFC| > 1.5. The results were presented in heatmap and volcano plots.

Gene ontology (GO) analysis is a common bioinformatics tool applied in large-scale functional enrichment studies that can annotate genes and analyze the biological process, cellular component, and molecular function of these genes (Yu et al., 2012). Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) is a database to explore the comprehensive biological systems and functions generated by experimental techniques in high-throughput biology from massive molecular datasets (Yu et al., 2012). The GO annotation and KEGG pathway enrichment analyses of signature genes was implemented using the ClusterProfiler package and the DAVID online database (Yu et al., 2012). Results with a false discovery rate (FDR) less than 0.05 were considered statistically significant.

To investigate the differences in biological processes between TP53-mutant and TP53-wt groups, we performed GSEA, based on the gene expression profile of the TCGA-BRCA dataset. GSEA is also a functional enrichment analysis, based on a predefined set of genes between two groups, which can determine whether there is a statistical difference (Subramanian et al., 2005). It is used frequently in analyzing the enrichment of signaling pathways and biological processes. The geneset of c2. cp.kegg.v6.2.-symbols was downloaded from the Molecular Signature Database (MsigDB, http://software.broadinstitute.org/gsea/msigdb/). GSEA was performed, and adjusted p-values less than 0.05 were regarded as statistically significant.

To quantify the immune cell proportions in the TCGA-BRCA samples, we used the CIBERSORT algorithm (https://cibersort.stanford.edu/) and the LM22 gene signature matrix (Newman et al., 2015). Highly sensitive and specific discrimination was performed for the phenotypes of 22 immunocytes (T cells, B cells, natural killer cells, and macrophages) in the TME (Hinshaw and Shevde, 2019). CIBERSORT was run to deconvolute samples, and used the expression values of a set of reference genes (547 genes), which were considered the minimal representative values for each type of cells. Based on these values, we deduced the cell type proportioning from the data of samples with mixed cells. Thus, we analyzed the effect of TP53 gene mutations on immune cell infiltration in TCGA-BRCA patients.

Meanwhile, the ESTIMATE algorithm was applied to assess the immune infiltration levels of BC patients according to the interpretation of gene expression profiles (Yoshihara et al., 2013). The ImmuneScore and StromalScore were calculated for each sample using the using the ESTIMATE package in R (https://www.r-project.org/). We performed Mann-Whitney U tests to compare the differences in immune cell infiltration in tumor tissues between TP53-mutant and TP53-wt groups.

In this study, we implemented the STRING (https://string-db.org) (Szklarczyk et al., 2019) to infer the protein-protein interaction (PPI) network. STRING is an online tool that can predict protein-protein interactions and construct the PPI network of selected genes. Interactions with a confidence score greater than 0.7 were included to construct the PPI network in Cytoscape software (Version 3.7.2). We defined the high-density areas as hub-genes based on the vertex-weighting scheme by using the MCODE plugin (Shannon et al., 2003).

The miRNA-mRNA interaction data was downloaded from the mirTarBase (http://mirtarbase.mbc.nctu.edu.tw/index.php) (Hsu et al., 2011). Then, we predicted the target miRNAs of the hub-genes from the PPI network, and carried out using the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw) (Hsu et al., 2011). Moreover, the regulatory relationships between miRNA and lncRNA were further established. Based on these hub gene-miRNA pairs and miRNA-lncRNA pairs, a ceRNA network for mRNA-miRNA-lncRNA was illustrated using Cytoscape software (version 3.7.2) (Shannon et al., 2003).

To validate the association between the expression of CXCL1, CXCL10, and CCL20 and TP53 status, we collected 10 cancer tissues with TP53 mutation and 10 tissues without mutation from BC patients. We performed IHC to compare the level of CXCL1, CXCL10, and CCL20 between two groups. IHC was performed as previously described (Wang et al., 2021), with antibodies specific for TP53 (Affinity, 1:100), CXCL1 (Affinity, 1:100), CXCL10 (Affinity, 1:100), or CCL20 (Affinity, 1:100). Pictures were taken with a microscope (Nikon DS-Ri2, Tokyo, Japan). Pathological samples were evaluated and scored separately by two qualified pathologists. The IHC scoring is as follows: 0 for no staining, 1+, 2+, 3+, and 4 + for 1–24, 25–49%, 50–74%, and over 75% staining intensity, respectively.

Genomics of Drugs Sensitivity in Cancer (GDSC) is a public online database (http://www.cancerrxgene.org/downloads/) and is used to determine anticancer drug response and somatic mutations in cancer (Yang et al., 2013). We identified the association between TP53 mutations and anticancer drug sensitivity in BC patients, based on the data of the gene mutation status in cancer cell lines and IC50 values of anticancer drugs.

To identify the prognostic role of TP53 status based on clinicopathological features, we analyzed the OS rate by conducting univariate and multivariate Cox regression analyses to test the risk score. The potential prognostic parameters were included to construct a nomogram using the TCGA-BRCA datasets. We constructed the nomogram using the rms R package. To analyze the performance of models, a calibration plot was graphically mapped by the nomogram predicted vs. observed probability. Moreover, the concordance index (C-index) was commonly obtained to quantitatively examine the discrimination ability of the nomogram.

In this study, all data processing and analysis were carried out using R software (Version 4.0.2). For continuous variables, a Student’s t-test was used to compare the means between the normally distributed variables, whereas a Mann-Whitney test was used for the variables that were not normally distributed. Moreover, a Chi-square test or Fisher exact test was used for discontinuous variables. The correlations among genes were determined by Pearson correlation analysis. Prognostic analysis was performed using the R package survival. The Kaplan-Meier curves were plotted to show the survival time of BC patients, and the log-rank test was used for the survival comparisons between the two groups. The independent prognostic factors in BC were identified using univariate Cox regression and multivariate Cox regression analyses. We plotted the receiver operating characteristic (ROC) curves using the pROC R package (Robin et al., 2011). The area under the ROC curve (AUC) was calculated to assess the prognostic risk scores (Robin et al., 2011). A two-sided p value less than 0.05 was considered as statistically significant.

To analyze the effects of TP53 mutations on the genomic mutations in BC patients, we downloaded three BRCA datasets from TCGA and ICGC databases (n = 943). First, we evaluated the mutation profile in BRCA patients as shown in Figure 1A. The results indicated that the missense mutations accounted for a major portion, single nucleotide polymorphisms (SNPs) were more often observed than insertion-deletion (indel) mutations, and the C > T single nucleotide variants were the most common variant in BC patients. The frequency of TP53 mutations was the second most in all the BRCA patients. Subsequently, we subdivided all patients into two groups, TP53-mutant and TP53-wt groups, according to the TP53 status. The somatic mutations of BRCA samples were calculated and visualized by the “Maftools” R package, and were presented in Figure 1B. The waterfall plots presented the mutation profile of associated genes (Figure 1B for TCGA-BRCA; Figure 1C for BRCA-FR; and Figure 1D for BRCA-KR). Moreover, the amino acid substitutions in the TP53 gene were evaluated and shown in Figure 1E. The location of each amino acid variant was corresponding to the coordinate axis below. The mutation type was distinguished by different colors, and the tag indicates the meaning of each color. The results showed that the main mutant form of TP53 amino acid was missense mutation in all three datasets. We also separated TCGA-BRCA patients into TP53-mutant and TP53-wt groupsTP53, and analyzed the CNV status. The data were analyzed via GISTIC 2.0 to obtain gene-level estimates of CNV, with the default settings except for several parameters (e.g., confidence: 0.99; X chromosome was not excluded from the analysis). Finally, the GISTIC 2.0 output was visualized using the MAFtools package, and shown in Figure 1F. This indicated that significant alterations in CNV levels located in related genes were observed in the TP53-mutant group.

Further, we explored the biological effect of TP53 mutations based on the mutational signature analysis. According to the biological characteristics, somatic mutational processes could be characterized by the mutation patterns, and 96 mutation patterns were translated into 30 different mutational signatures (Alexandrov et al., 2020). The results indicated that significant changes in Signature 1, 3, and 13 were observed compared to those in the TP53-wt group (Figures 2A,B). In addition, compared with that in the TP53-wt group, the level of the TP53 gene was substantially increased in the TP53-mutant group (p = 0.037; Figure 2C), and the TMB value (p < 0.001; Figure 2D), MSI value (p = 0.077; Figure 2D), and TIDE score for immunotherapy (p < 0.001; Figure 2E) were also elevated in the TP53-mutant group.

To detect the effect of TP53 mutations on drug sensitivity in BC patients, we assessed the correlation between TP53 mutations and IC50 values of molecules from the GDSC database. The result showed that multiple drugs related to the frequency of TP53 mutation (Figure 3A). The pathway analysis revealed that the TP53 pathway was significantly enriched (Figure 3B), and the high mutation rates of 6 genes in this pathway were also prevalent in BC patients (Figure 3C). Moreover, the TP53 mutation had some effect on BC sensitivity to multiple chemotherapy agents and small molecule substances (Figure 3D), especially to Nutlin-3a (Figure 3E).

To assess the effect of the TP53 mutation on BC tumorigenesis, the TCAG-BRCA patients were separated into TP53-mutant and TP53-wt groups. As shown in Table 1, TP53 mutation status was significantly correlated with a younger age (<60 vs. ≥ 60, p = 0.007) and earlier M stage (M0 vs. MX, p = 0.007). As evaluated by PCA analysis, significant differences were shown (p < 0.05) between TP53-mutant and TP53-wt groups (Figure 4A). Moreover, DEGs analysis identified that 845 upregulated DEGs and 237 downregulated DEGs were associated with the TP53 mutation (|log2 fold change|> 1.5 and (adjust) p-value < 0.01; Figures 4B,C).

Subsequently, to analyze the cellular functions of 1082 DEGs, we conducted GO and KEGG enrichment analyses using the R package clusterProfiler. The results of the GO analysis demonstrated that DEGs were involved in the biological processes of cornification, keratinization, skin development, intermediate filament cytoskeleton, and peptidase inhibitor activity (Table 2; Figures 4D–F). KEGG pathway analysis suggested that the immune-related DEGs were significantly enriched in neuroactive ligand-receptor interaction, nicotine addiction, salivary secretion, and the IL-17 signaling pathway (Table 3; Figure 4G). Besides, the GSEA results of the TP53-mutant group revealed that the significant pathways (p < 0.05 and FDR q-value < 0.25) were enriched in focal adhesion, ribosome, antigen processing and presentation, and ECM receptor interaction, and the details are shown in Table 4 and Supplementary Figure S1A, S1B.

The PPI network of DEGs was constructed using the STRING online database (Supplementary Figure S2A), and the results were imported into Cytoscape software for further analysis (Supplementary Figure S2B); the red color represented up-regulated gene expression and the blue color represented down-regulated gene expression. Then, we used the plugin MCODE in Cytoscape to analyze the important modules. In the regions of high density, the central nodes were identified as hub-genes (Supplementary Figure S2C). Based on the information of miRNA-mRNA interactions in the miRTarBase, we predicted the miRNAs associated with the hub-genes, and lncRNA associated with the miRNAs. Thus, the mRNA-miRNA-lncRNA ceRNA network was constructed based on the predicted relationship shown in Supplementary Figure S2D. The results above indicated that the chemokines CXCL1, CXCL10, and CCL20 was significantly upregulated in the TP53-mutant group (Supplementary Figure S2C, S2D). Further, the CXCL10 and CCL20 expression level was lower in BC tissues (TCGA-BRCA patients) compared with normal tissues (Figures 5A,B). We also examined the expression level of them in paired tissue samples. The results indicated that the level of CXCL10 and CCL20 in BC tissues was also significantly lower than those in paired samples (Figures 5C,D).

We determined the effects of TP53 mutation on the expression of hub-genes by IHC in BC tissues, including CXCL1, CXCL10, and CCL20. As shown in Figure 5E, the upregulated expression of CXCL1, was significantly associated with TP53 mutation (p < 0.05). The similar results were also found for the expression of CXCL10, CCL20, and TP53 (p < 0.05; Figure 5E).

To determine how the TP53 mutation influences BC immunogenicity, we compared the expression differences in immune-related genes and stromal-related genes between the TP53-mutant group and TP53-wt group. The results indicated that in the mutation group, the levels of the ImmuneScore and StromalScore were both significantly increased (ImmuneScore: p < 0.001; StromalScore: p = 0.003; Figure 6A). Moreover, the expression of multiple HLA gene families was significantly upregulated in the mutation group (Figure 6B). Next, we used CIBERSORT to evaluate the composition ratio of 22 immune cell types in each BC sample and the result showed individual differences (Figure 6C). We also compared the levels of 22 immune cells between the TP53-mutant group and TP53-wt group. The results demonstrated that the proportions of Tregs, T helper cells, and M0 type macrophages were significantly upregulated in the TP53-mutant group (Figure 6D, p < 0.05), whereas the proportion of resting CD4⁺ T cell and M2-type macrophages was lower (Figure 6D, p < 0.05).

We also performed the Kaplan-Meier analysis to assess the prognostic significance of TSPOAP1-AS1 expression. In the TCGA-BRCA patients, the TP53 mutation was significantly associated with a shorter OS (p = 0.038; Figure 7A), whereas there was no significance for BRCA-FR (p = 0.819; Figure 7B) and BRCA-KR (p = 0.301; Figure 7C) patients. Then, to further confirm the prognostic value of the TP53 mutation, we conducted univariate and multivariate Cox regression analyses for OS. The results revealed that TP53 mutation (p = 0.0298), age (p < 0.001), tumor stage (p < 0.001), T stage (p = 0.01), N stage (p < 0.001), and M stage (p < 0.001) were correlated with BC prognosis (Table 5). Then, these variables were included to build the multivariable Cox models of OS (Table 5). The TP53 mutation remained independently associated with OS [HR: 1.76 (1.24–2.50), p = 0.002], which was also true for age [HR: 1.94 (1.37–2.76), p < 0.001], tumor stage [HR: 2.46 (1.26–4.80), p = 0.009], and M stage [HR: 1.67 (1.05–2.66), p = 0.03]. These results revealed that the TP53 mutation is an independent predictive factor of OS in BC patients. Further, to develop a clinical quantitative tool to predict the OS for BC patients, a nomogram was constructed based on the results of multivariable cox regression. In this nomogram, the significant variables including the TP53 mutation, age, stage, and TNM status were used to assign points (Figure 7D). The C-index of this nomogram was 0.772, and the calibration plots suggested that there was good consistency between the nomogram and observed OS probabilities in BC (Figure 7E).