Development of a prognostic signature based on immune-related genes and the correlation with immune microenvironment in breast cancer

Data acquisition and preprocessing

The RNA-sequencing datasets and the corresponding clinical characteristics data of BC patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). Patients with less than 30 days of follow-up and male BC patients were excluded. Finally, 1039 patients were included in this study and randomized into training group (n = 520) and testing group (n = 519) by using the R package “caret” for subsequent analysis [11]. The training set was used to construct the prognostic immune gene signature, while the entire set and testing set were used to validate the predictive ability of the established prognostic immune model. The list of IRGs was downloaded from the Immunology Database and the Analysis Portal (ImmPort) database (https://www.immport.org/) [12], an open platform of human immunology database for clinical and translational research. These genes could be identified and extensively participate in the important process of immunology, thus providing a good foundation for immunology research (Figure 1).

Figure 1. The flowchart of this study.

Differentially expressed genes (DEGs) and IRGs analysis

Then, we performed the DEG analysis on the transcriptional data from the TCGA database with cutoff values of false discovery rate (FDR) < 0.05 and log2 |fold change| > 1. After that, we extracted the differentially expressed IRGs from the above DEGs and used the limma package of R software (http://bioconductor.org/packages/limma/) to acquire the differentially expressed IRGs associated with BC between cancer and adjacent non-tumor samples obtained differentially expressed IRGs associated with BC using the limma package of R software (http://bioconductor.org/packages/limma/). Finally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway analysis and the Gene Ontology (GO) annotation were utilized to explore the underlying molecular mechanism of these differentially expressed IRGs.

Extract Transcription factors (TFs) and construct the regulatory network

The Cistrome Cancer Database summarizes the TCGA genome data accompanied with more than 23,000 chip-SeQ and chromatin accessibility maps to demonstrate regulatory relationships between TFs and genes [13]. The OS time of patients was obtained by downloading clinical information from the TCGA database. The univariate Cox analysis was performed by using the R survival package to screen out the IRGs that were closely related to survival prognosis. Accordingly, regulatory network was constructed for seeking the interrelation between these genes and TFs. In this study, TFs were compared with previously obtained IRGs to screen differentially expressed TFs related to survival-related IRGs. The regulatory network diagram of TFs and survival-related IRGs was drawn using Cytoscape software (version 3.8.0) [14].

Construction and verification of the immune-related prognostic signature

To construct an IRG-based immune-related prognostic model, genes significantly associated with prognosis were screened by univariate Cox regression analysis, and risk coefficients were obtained by Lasso-Penalized Cox regression analysis. For predicting OS, BC patients were divided into a high-risk group and low-risk group based on the calculation results of the R packages “survival” and “survminer”, using the median risk score as the best cut-off value. The time-dependent prognostic capability of the gene signature was measured by calculating the area under the curve (AUC) [15], using the R package “survivalROC” [16]. In addition, univariate and multivariate Cox regression analyses were performed to evaluate the prognostic significance and routine clinicopathological features, including age, grade, clinical stage, and TNM stage. These results were presented through the R package “Ggpubr” [17]. The principal components analysis (PCA) is used as a statistical method to find key variables in multidimensional datasets. PCA could predict the analysis and visualize of multidimensional data sets [18] through limma [19] and scatterplot 3d [20] packages. Subsequently, PCA was performed to verify the grouping ability of signatures by identifying a small set of synthetic variables through a dimensionality reduction process. Finally, the m6Ascore was determined in a manner similar to the Genomic Grade Index (GGI) [21]: riskScore = ∑ (PC1_i + PC2_i), where i was the expression of overlapping IRGs with a significant difference in prognosis of BC patients.

Gene set enrichment analysis (GSEA)

GSEA is a computational method for identifying whether defined gene sets are statistically significantly and consistently different across biological states [22]. We performed GESA using the JAVA program (https://www.broadinstitute.org/gsea), to explore the potential KEGG pathway enrichment terms of the IRG signature in the BC cohort. A total of 1000 random sample permutations were included, and enriched gene sets with nominal p < 0.05 and FDR < 0.25 were considered statistically significant.

Comparison of TME cell infiltration among high- and low-risk group

To understand the degree of immune cell infiltration in the two subgroups, the relative abundance of each cell infiltration in the TME of BC samples was quantified using single-sample gene-set enrichment analysis (ssGSEA). From the research of Charoentong, we obtained the gene sets for each type of TME infiltrating immune cells, and stored a relatively comprehensive subset of human immune cells, including activated CD8 T cells, natural killer T cells, activated DCs, macrophages, and regulatory T cells.

BC sample validation experiments

Here, we adopted the immunofluorescence (IF) and immunohistochemistry (IHC) to verify one of the model-associated IRGs (SDC1) in the collected BC samples. The BC tumor samples were Surgically excised in our department of thyroid and breast surgery (Tongji Hospital). Specifically, for immunohistochemistry (IHC), all the BC tissues were deparaffinized and were heated in citrate buffer in sequential. Then, the obtained sections were immune-stained with a primary anti-SDC1 antibody (CD138/Syndecan-1 Rabbit pAb, ABclonal, China Catalog: A1235) overnight at 4°C. After washing, these incubated sections were then incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies. DAB peroxidase substrate Kit (Maxin, China) was used to observe peroxidase activity and sections were restained with hematoxylin. Digital images of sections were collected by SOPTOP CX40 microscope (China). The IF experiment was conducted to verify the distribution of SDC1 and its correlation with M2 macrophages in the BC samples. For IF, the obtained sections were incubated with the above anti-SDC1 antibody, anti-CD2062 antibody, and nuclear 4,6-diamidino-2-phenylindole (DAPI, Sigma, USA) for counterstaining. Digital image acquisition was realized by using a fluorescence microscope (Olympus, Japan).

Statistical analysis

All the statistical analyses were performed by the R (v.3.6.3) software. The Fisher’s exact test or Pearson χ² test was used to analyze qualitative variables as appropriate. P < 0.05 was considered statistically significant.

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Data sources and identification of differentially expressed IRGs

A total of 4575 DEGs were identified from 1109 BC samples and 113 normal samples, including 2698 up-regulated and 1877 down-regulated genes (Figure 2A, 2B). Among them, 366 differentially expressed IRGs (193 up-regulated and 173 down-regulated) were obtained from these DEGs by immune gene co-expression analysis using the ImmPort database (Figure 2C, 2D). The result of GO analysis showed that the differentially expressed IRGs were enriched in several biological processes, including cellular response to chemokine, leukocyte migration, and chemokine-mediated signaling pathway. Cellular component analysis demonstrated that IRGs were mostly enriched in T cell receptor complex, immunoglobulin complex, and circulating, cytoplasmic vesicle lumen. And in the molecular function, these genes were mostly enriched in receptor-ligand activity and cytokine activity (Figure 2G). The KEGG analysis indicated that IRGs were mainly involved in Cytokine-cytokine receptor interaction (Figure 2H).

Figure 2. Differentially expressed genes, immune-related genes, and TFs. The differentially expressed genes between BC and non-tumor tissues were shown in the heatmap (A) and the volcano plot (B). Heatmap (C) and the volcano plot (D) indicated the differentially expressed immune-related genes (IRGs). Heatmap (E) and volcano plot (F) referred to the differential TFs between breast cancer and non-tumor tissues. Red dots represented the upregulated genes or TFs, green dots represent downregulated genes or TFs, and black dots represented not differentially expressed genes or TFs. N, normal tissue. T, tumor. The GO (G) and KEGG (H) functional enrichment analysis of immune-related genes (IRGs).

Identification of survival-related IRGs and construction of the TF regulatory network

By integrating the mRNA expression and clinical information of BC patients in the TCGA database, we finally screened 39 IRGs with significant differences in the prognosis of BC patients (P < 0.05) (Supplementary Table 1). Then, to further explore the potential molecular mechanism of these survival-related IRGs (Supplementary Figure 1A), the expression patterns of 318 TFs were examined and 80 TFs were found to be differentially expressed in BC and non-tumor breast tissues in the TCGA database (Figure 2E, 2F). Subsequently, we analyzed the relationship between prognostic immune genes and differentially expressed TFs, and screened out 22 TFs and 23 prognostic immune genes, thus constructing a regulatory network with a correlation score > 0.4 and p-value < 0.001. In particular, the visualization of TF-based regulatory network clearly showed the interactions between these genes (Supplementary Figure 1B).

Construction of the immune-related prognostic model

BC patients were randomly divided into a training set and a test set, including 520 and 519 BC patients, respectively. By performing univariate Cox regression on the training data, it found that 40 IRGs were retained in the training data with p < 0.05 (Table 1, Supplementary Table 2). The Cox regression analysis and least absolute shrinkage selection operator (LASSO) regression were used to analyze the expression profiles of differentially expressed IRGs in the training group, to define the candidate genes (Figure 3A, 3B). Finally, a total of 6 IRGs were obtained, including PSME2, ULBP2, IGHE, SCG2, SDC1, and SSTR1. The specific calculation formula was as follows: [Expression level of PSME2^* (−0.0143)] + [Expression level of ULBP2^* (0.1070)] + [Expression level of IGHE^* (0.0619)] + [Expression level of SCG2^* (0.0044)] + [Expression level of SDC1^* (0.0021)] + [Expression level of SSTR1^* (0.0451) (Table 2). Based on these prognostic genes, we successfully established the immune-related prognosis model for OS.

Figure 3. The LASSO coefficient profiles about 6 immune-related genes were shown in (A, B). The lower X-axis indicated log (λ), the upper X-axis indicated the average number of OS-related genes, and the Y-axis showed the partial likelihood deviance error. Red dots indicated the average partial likelihood deviances about the model with a given λ, the vertical bars represented the range of the partial likelihood deviance errors, and the vertical black dotted lines meant the best fit with the optimal λ values. The training set (C), the testing set (D) and the entire set (E) showed the survival curves of high-risk group and low-risk group respectively. The survival-dependent receiver operating characteristic (ROC) curves prognostic value in 1 year of the three sets were shown in the (F–H).

Validation of the immune-related prognostic model

Then, we examined the comparisons of survival differences between the high- and low-risk group in the training set (Figure 3C), the testing set (Figure 3D), and the entire set (Figure 3E). The Kaplan-Meier log-rank analysis revealed significant differences in OS between the two risk groups. In addition, we found that the AUCs of OS in training set, test set and entire set were 0.782, 0.888, and 0.775, respectively (Figure 3F–3H). Our IRG model had the maximum AUC value compared to other clinical features, demonstrating good predictive power. Subsequent risk curve, scatter diagram, and heatmap were used to analyze the risk score distribution (Figure 4G–4I), survival status (Figure 4J–4L), and the expression of the candidate genes (Figure 4M–4O) for each BC patient in the training set, testing set, and the entire set, respectively. The univariate and multivariate Cox regression analysis verified that the IRGs in the prognosis model could serve as independent predictors of prognosis (Figure 4A–4F).

Figure 4. The univariate (A–C) and multiple (D–F) regression analysis of BC, including the relationships of the age, stage, T stage, distant metastasis, lymph node metastasis, and riskScore in the three sets. The training set (A, D), the testing set (B, E), and the entire set (C, F). The green squares indicated that the median value of hazard ratio (HR) was less than 1, while the red squares indicated that the median value of HR was greater than 1. Analysis of risk score, OS, and the expression of the six genes in the training set (G–I), testing set (J–L), and entire set (M–O). The risk score, OS, and heat map were listed from top to bottom.

Next, by analyzing the model and different clinicopathological factors, the prognosis of low-risk group was significantly better than that of high-risk group in terms of age (≤ 65/ > 65) (Supplementary Figure 2A and 2B), T stage (T1 + T2/T3 + T4) (Supplementary Figure 2C and 2D), and stage (I + II/III + IV) (Supplementary Figure 2E and 2F). Similarly, the patients without lymph node metastasis (Supplementary Figure 2G) and distant metastasis (Supplementary Figure 2H and 2I) had a similar outcome. The PCA was further to detect the difference between the low-risk group and the high-risk group according to the immune-related signals (Figure 5A), immune genes (Figure 5B), differential expression genes (Figure 5C), and the entire gene expression profiles (Figure 5D). As seen in Figure 5, compared with the other three groups, the high-risk group and the low-risk group in Supplementary Figure 2A were distributed in different directions, proving that our prognosis model could effectively distinguish the high-risk group from the low-risk group.

Figure 5. The PCA based on the immune-related signature (A), immune-related genes (B), differently expressed genes (C), and the entire gene expression profiles (D) between the high-risk group and low-risk group. The correlation of the immune-related signature with clinicopathological characteristics. SDC1 was associated with age (E) and lymph node metastasis (F). IGHE were associated with age (G). The risk score of our prognostic model was significantly associated with a higher tumor stage (H).

Assessment of the correlation between candidate genes and clinicopathological characteristics

We then analyzed the relationship between the expression of 6 candidate IRGs and different clinical features, including age, gender, tumor grade, T stage, clinical stage, lymph node metastasis, and distant metastasis. Among 6 candidate IRGs, IGHE and SDC1 were correlated with age, but IGHE expression was increased in patients over 65 years old, while SDC1 expression was increased in patients under 65 years old (Figure 5E, 5G). Besides, the expression difference of SDC1 between lymph node metastasis groups was statistically significant (Figure 5F). Finally, the risk score based on our prognostic model was significantly associated with the higher tumor grade (Figure 5H).

Different states of functional enrichment analysis between high-risk and low-risk groups

GSEA was used to investigate the differences between the high- and low-risk groups. The results revealed that the GO biological process “Regulation of cell substrate junction organization” (Figure 6A) and “Regulation of chondrocyte differentiation” (Figure 6B), molecular function “Extracellular matrix structural constituent” (Figure 6C) were differentially enriched in low-risk groups (P < 0.01). Besides, the biological process “Glycosyl compound catabolic process” (Figure 6D), molecular function “Oxidoreductase activity acting on NADPH” (Figure 6E) and cellular component “Oxidoreductase complex” (Figure 6F) were associated with the high-risk groups (P < 0.01).

Figure 6. Enrichment plots of Gene Ontology annotation from gene set enrichment analysis (GSEA). GSEA results showed that the regulation of cell substrate junction organization (A), Regulation of chondrocyte differentiation (B), and Extracellular matrix structural constituent (C) were differentially enriched in low-risk phenotype, while Glycosyl compound catabolic process (D), Oxidoreductase activity on NADPH (E) and Oxidoreductase complex (F) were enriched in the high-risk phenotype. The expression abundance of different TME infiltrating cells in the high- and low-risk group (G). The upper and lower end of the box represented the quartile range of the value, the middle line represented the median value, and the asterisk represented the statistical p-value (^*P < 0.05, ^**P < 0.01, ^***P < 0.001).

Besides, the KEGG pathway analysis showed that the genes of low-risk group were mainly enriched in the “TGF beta signaling pathway” (Supplementary Figure 3A), “Hedgehog signaling pathway” (Supplementary Figure 3B), and “Adherens junction” (Supplementary Figure 3C) (P < 0.01) while the “Antigen processing and presentation” (Supplementary Figure 3D), “Cytosolic DNA sensing pathway” (Supplementary Figure 3E) and “Oxidative phosphorylation” (Supplementary Figure 3F) were enriched in high-risk group (P < 0.01). The function enrichment of IRGs between high- and low-risk groups was explored, showing that the IRGs in our model were mainly involved in immune-related signaling pathways (Supplementary Figure 1C, 1D).

The tumor-infiltrating immune cells in risk signature

Subsequently, to explore the relationship between the IRG-based prognostic risk model and TME, the differences in tumor-infiltrating immune cells between the high- and low-risk groups defined in our prognostic model were analyzed. It was found that the abundance of activated B cells, activated CD4 T cells, immature DCs and eosinophils were significantly more abundant in the low-risk group of the OS model (Figure 6G). These results potentially shed light on the regulatory mechanisms of BC TME.

The immune-related risk signature and mutation profile

Using somatic mutation data available from the TCGA database, the relationship between mutation characteristics and this model was evaluated in BC patients. The frequently mutated genes in the high- and low-risk groups were presented in Figure 7A, 7B. It was intriguing that tumor mutation burden (TMB) was significantly higher in low-risk groups and associated with longer OS (Figure 7G).

Figure 7. The waterfall diagram of tumor somatic mutation in patients with high- (A) and low- (B) risk groups. (C–F) The association between IPS and the risk model, the IPS, IPS-PD1, IPS-CTLA4, and IPS-PD1/CTLA4 scores were significantly increased in the low-risk group. Kaplan-Meier curve was used to analyze the survival of the high and low TMB load (G). The results showed that the survival of the two cohorts with high and low TMB load was significantly different. (H, I) Wilcoxon test was used to analyze the difference of PD1 and CTLA4 expression between high- and low-risk groups.

The immune-related risk signature and response to immune checkpoints-inhibitors (ICIs)

At last, we further explored the relationship between IPSips and this prognostic model. The IPS, IPS-PD1, IPS-CTLA4 and IPS-PD1/CTLA4 scores were designed to evaluate the feasibility of ICIs applying for BC patients. This result showed that the IPS was significantly increased in the low-risk group compared with the high-risk group (Figure 7C–7F). Moreover, the expression of PD1 and CTLA4 was higher in the low-risk group (Figure 7H, 7I). These results collectively suggested that IPS levels were higher in the low-risk group, and that these patients exhibited more immunophenotypes and were peculiarly prone to benefit from immune checkpoint therapy.

The validation of immune-related SDC1 expression in BC samples

In the end, the IHC and IF were performed to verify the SDC1 expression characteristics, which was one of the 6 IRGs in our risk model. The IHC results showed that SDC expressions were significantly higher in BC samples of high-risk patients, compared with those of low-risk patients (Figure 8A). Additionally, IF results showed the SDC1, as well as M2 macrophage biomarker CD206, were both of higher abundance in BC samples of high-risk patients, suggesting that SDC1 was indeed a tumor-promoting factor as previously reported and was positively linked to the expression of M2 macrophages (Figure 8B). The histological verification of BC samples successfully confirmed that the expression feature of IRGs in our model was in line with expectations.

Figure 8. The SDC1 expression features in BC samples. (A) The IHC assay showed the SDC1 expression level in the high-risk and low-risk patients. (B) The IF assay verified the SDC1 and CD206 expression levels, as well as their co-expression in the high-risk and low-risk patients. CD206, green; SDC1, red; DAPI, nucleus.