1. Analysis process of the study
We first downloaded the transcriptome profiling and clinical data of 414 BC samples and 19 normal samples from TCGA database and further analyzed the ratio of immune and stromal component in TME of each tumor sample through ESTIMATE algorithm. Difference analysis was conducted to find out the association between immune/stromal component in TME and the clinic-pathological characteristics of BC patients. A total of 284 DEGs were further acquired based on the immune and stromal component in TME of BC. GO and KEGG enrichment analysis of these 284 DEGs were performed to learn their biological functions. Finally, CXCL12 was found out through the intersection analysis of PPI network and univariate COX regression. Then, we focused on the expression profile of CXCL12 in BC samples and normal samples, and the association between the expression of CXCL12 and survival and clinic–pathological characteristics of BC patients. GSEA of different gene collections was also carried out to learn the functional profiles of CXCL12. TICs profile in TIME of BC samples was calculated by CIBERSORT algorithm. Difference analysis and correlation analysis were applied to find out the correlation between the expression of CXCL12 and TICs.
2. Immune and stromal component in TME were associated with the clinic-pathological characteristics of BC
To learn the association between the proportion of immune/stromal component in TME and the clinic-pathological characteristics of BC patients, we first set up the immune/stromal component evaluating system of TME in BC samples. Immune, Stromal and ESTIMATE Score represented the proportion of immune, stromal and the total component in TME of each tumor sample, respectively. The higher Score suggested the larger proportion of immune or stromal component in TME. The clinic-pathological characteristics of BC patients, including age, gender, grade, stage and TNM classification, were concluded in Supplement Table 1. Difference analysis revealed that Stromal Score was associated with patients’ age and gender, especially the grade (Figure 2a-c, p=0.018; 0.032; <0.001). The higher Stromal Score indicated the higher grade of BC (Figure 2c, p<0.001). It was also found that Stromal Score was related with T, N classification and clinical stage of BC (Figure 2d-g). Regarding Immune Score, it was significantly up-regulated in female patients (Figure 2i, p=0.037) and the higher Immune Score also suggested the higher grade of BC (Figure 2j, p<0.001). However, neither TNM classification nor clinical stage of BC was associated with Immune Score (Figure 2k-n, p>0.05). As for ESTIMATE Score, the results similarly showed that it was up-regulated in female patients (Figure 2p, p=0.028) and indicated the higher grade of BC (Figure 2q, p<0.001). In addition, ESTIMATE Score was connected with T classification and clinical stage of BC (Figure 2r-u). From above, we might conclude that the stromal and immune component in TME were significantly associated with the grade of BC, and partly related with TNM classification and clinical stage. This provided us a new sight to explore the underlying mechanisms of the development and progression of BC.
3. DEGs were primarily enriched in immune-associated activities
In order to further understand the underlying mechanisms of BC TME, we conducted difference analysis to acquire the DEGs profile in TME. BC samples were first divided into two groups, including high Stromal Score group and low Stromal Score group (or high Immune Score group and low Immune Score group). Heat map displayed the DEGs profile in TME of BC (Figure 3a-b). To be specific, there were 537 DEGs obtained from the comparison between high Stromal Score and low Stromal Score, among which 496 DEGs were up-regulated and 41 DEGs were down-regulated (Figure 3c, Supplement Table 2). Besides, there were 531 DEGs obtained from the comparison between high Immune Score and low Immune Score, among which 466 DEGs were up-regulated and 65 DEGs were down-regulated (Figure 3c, Supplement Table 3). In conclusion, a total of 284 DEGs were commonly up-regulated or down-regulated in Stromal or Immune Score group (Figure 3c, Supplement Table 4). Next, we carried out GO and KEGG enrichment analysis to assess the biological functions of these 284 DEGs. Go enrichment analysis revealed that 284 DEGs were primarily enriched in immune-associated activities, such as leukocyte proliferation, migration and chemotaxis, lymphocyte and mononuclear cell proliferation and so on (Figure 3d-f). Similarly, KEGG enrichment analysis suggested that these DEGs were mainly enriched in immune-associated activities, including complement and coagulation cascades, cytokine-cytokine receptor interaction, B cell receptor signaling pathway, chemokine signaling pathway and so on. Therefore, we considered that DEGs acquired from Stromal and Immune Score group of BC were significantly associated with immune-associated activities in TME.
4. PPI network and univariate COX regression analysis of 284 DEGs
To learn the detailed reciprocity among these 284 DEGs, we further ploted PPI network through STRING database and cytoscape software. The potential interactions among DEGs were displayed in Figure 4a. The top ten DEGs, which shared the leading nodes, were ranked in Figure 4b, including ITGAM, CXCL12, CXCL13 and so on (Supplement Table 5). Univariate COX regression of DEGs was conducted at the same time to find out the specific DEGs, which were significantly associated with BC patients’ survival (Figure 4c). MMP9, COMP, F13A1 and CXCL12 came to our eyes (Figure 4c, p<0.05; Supplement Table 6). Intersection analysis was finally carried out to search the DEGs who shared the leading nodes in PPI network and were significantly related with BC patients’ survival (Figure 4d). Fortunately, CXCL12 emerged.
5. CXCL12 was down-regulated in BC tissues and associated with the clinic-pathological characteristics
The CXC chemokine CXCL12 greatly participated in multiple physiological and pathological processes through interacting with its receptors CXC chemokine receptor 4 (CXCR4) and atypical chemokine receptor 3 (ACKR3) [29]. In our study, CXCL12 was significantly down-regulated in BC tissues compared with normal tissues (Figure 5a, p<0.001). Besides, the paired analysis also indicated that CXCL12 was obviously down-regulated in BC tissues compared with the paired normal tissues (Figure 5b, p<0.001). However, survival analysis suggested that the expression of CXCL12 was not significantly associated with the survival of BC patients (Figure 5c, p=0.050). As for clinic-pathological characteristics of BC patients, the results showed that the expression of CXCL12 changed with ages and higher expression of CXCL12 indicated the higher grade of BC (Figure 5d, f; p<0.001, p<0.001). In addition, the expression of CXCL12 was partly correlated with T, N classification and clinical stages of BC (Figure 5h-k). Hence, we might conclude that CXCL12 had the potential to participate in the progression of BC.
6. CXCL12 greatly participated in immune-associated activities
On account of the potential role of CXCL12 in the progression of BC, we further explored the underlying mechanisms of CXCL12. BC cases were first divided into two groups, including high CXCL12 expression group and low CXCL12 expression group. GSEA suggested that for C5 collection, the gene ontology sets, genes in high CXCL12 expression group were enriched in cytokine binding, leukocyte migration and negative regulation of interleukin 6 production (Figure 6a, p<0.05). For C2 collection, the KEGG gene sets database, genes in high CXCL12 expression group were enriched in allograft rejection, antigen processing and presentation, B cell receptor signaling pathway, complement and coagulation cascades, and so on (Figure 6b, p<0.05). For hallmark gene sets, genes in high CXCL12 expression group were primarily enriched in allograft rejection, complement, IL2 stats signaling, IL6 jak stat3 signaling and so on. Finally, for C7 collection defined by MSigDB, the immunologic gene sets, genes were found to be enriched in multiple immune-associated activities, which were related with CD4 T cell, CD8 T cell, B cell and so on (Figure 6d, p<0.05). As a result, according to GSEA, CXCL12 was significantly related with immune-associated activities and might play an important role in regulating TIME of BC.
7. CXCL12 communicated with TICs in TIME
Since GSEA suggested that CXCL12 was significantly correlated with immune-associated activities, we speculated that there might be underlying connections between CXCL12 and TIME of BC. In consequence, we calculated the proportion of TICs in BC samples through CIBERSORT algorithm. The specific ratio of 22 TICs in tumor samples was shown in Figure 7a. And the association among these TICs was exhibited in Figure 7b. Difference analysis was carried out to learn the association between the expression of CXCL12 and specific TICs. The results showed that 10 kinds of TICs, such as Macrophages M2, B cells naïve, T cells gamma delta, T cells CD4 naïve, T cells follicular helper and so on, were significantly associated with the expression of CXCL12 (Figure 7c, Supplement Table 7). The correlation analysis suggested that the expression of CXCL12 were positively or negatively related with 8 kinds of TICs, among which B cells naïve, macrophages M2 and mast cells resting were positively associated, and dendritic cells activated, dendritic cells resting, NK cells resting, T cells CD4 naïve and T cells follicular helper were negatively associated (Figure 7e-l, p<0.05; Supplement Table 7). Finally, intersection analysis between difference analysis and correlation analysis suggested that 6 kinds of TICs were significantly related with the expression of CXCL12 (Figure 7d). In conclusion, CXCL12 was obviously related with various TICs in TIME of BC and might greatly participated in regulating the TIME of BC via closely communicating with multiple TICs.