Association between the TME, the TMB, and prognosis in BRCA patients.
The immune scores were divided into high-immune and low-immune score groups according to the optimal cut-off value. The K-M survival analysis and log-rank test were performed to identify the relationship between either OS or the PFI with the immune score. The survival probability of the high-immune score group was significantly higher than the low-immune score group, with a p-value of 0.0028 [Fig. 2A]. The probability of a progression-free interval was also significantly higher for the high-immune score group than the low-immune score group (p-value = 0.0084) [Fig. 2D]. Similarly, the stromal score and the TMB were classified into high-stromal and low-stromal score groups and high-TMB and low-TMB groups, respectively. The K-M survival analysis and log-rank test demonstrated that patients with a high-stromal score had no statistical difference with either OS [Fig. 2B] or the PFI [Fig. 2E] compared to patients with a low-stromal score. While patients with high-TMB have a lower survival probability than patients with low-TMB (P-value = 0.0055) [Fig. 2C], patients with high-TMB had no statistical difference with PFI compared to patients with low-TMB (P-value = 0.86) [Fig. 2F]. Patients in the high-immune score group had a significantly higher stromal score and TMB than patients in the low-immune score group [Fig. 3].
Construction of the WGCNA and identification of corresponding modules.
Setting the criterion of protein-coding genes expressed in at least half of BRCA tissues, a total of 13721 mRNA sequences were screened for the WGCNA. The screen identified 12 modules [Fig. 4A]. The association between modules and traits was constructed and identified correlations between the green module and immune-score (0.94) and the black module and stromal-score (0.77). Both p-values were less than 0.001 [Fig. 4B]. In the present study, eight was defined as a set point for the soft-threshold power and the related scale-free topology index was 0.9 [Fig. 4C]. Therefore, we further analyzed the green and black modules. The relationships between GS and MM were identified for both the green and black modules, the correlation was 0.99, 0.85, respectively (p-values < 0.001) [Fig. 4D, E]. In the green and black modules, 753 and 180 mRNA sequences were identified, respectively. We defined the cut-off as a MM > 0.8 and a GS > 0.6, narrowing our results for further study to 191 and 40 mRNA sequences in the green and black modules, respectively.
Calculation of the risk score by random forest, univariate, and multivariate Cox analysis.
Random forest was performed for a survival analysis of the 191 mRNA sequences from the green module and 40 mRNA sequences from the black module, respectively. A univariate Cox analysis was used to examine the relationships between the expression of the 191 mRNA sequences in the green module, and the expression of the 40 mRNA sequences in the black module, with OS. Based on these results, 37 mRNA sequences [Fig. 5A] were selected for a stepwise multivariate Cox analysis. In addition, a four-mRNA risk score formula was established [Fig. 5B]. Risk score= (-0.259*expression level of “MYO1G”)+ (-0.27*expression level of “TBC1D10C”)+(0.292*expression level of “SELPLG”)+ (-0.11*expression level of “LRRC15”). Based on the median value, the risk score was classified into high-risk and low-risk groups. A K-M survival analysis and log-rank test demonstrated that the high-risk group had a lower probability of survival than the low-risk group (p-value < 0.001) [Fig. 5C].
Development of the nomogram for predicting survival probability.
A multivariate Cox regression was performed to construct the prediction model [Fig. 6A]. The nomogram illustrates that the M-stage made a robust contribution to the prediction of survival probability, followed by age, the tumor size, and the risk score. The score for each variable subtype was presented on a point scale. By calculating the total score and identifying it on the scale for the total possible points, we can easily obtain the survival probability of a patient. DCA (Decision Curve Analysis) is an emerging method to evaluate the discrimination of a predictive model(21), it is widely used in many top journals, such as BMJ, JAMA, NATURE, and others. In our study, DCA to evaluate the predictability of our model demonstrated that the net benefit of our nomogram was significantly higher than others [Fig. 6B]. The C-index for the nomogram was 0.713, and the calibration curve demonstrated good performance by the nomogram [Fig. 6C].
Relationship between four TME genes and BRCA endpoints.
We divided four genes (MYO1G, TBC1D10C, SELPLG, and LRRC15) into high and low groups based on the median expression value. The K-M analysis and log-rank test were applied. SELPLG and MYO1G expression had statistically significant correlations with OS in both the high and low expression groups, but not with the PFI [Fig. 7A, B, E, F]. The differential expression of LRRC15 was not correlated with either OS or the PFI [Fig. 7C, G]. The TBC1D10C high-expression group had a higher positive association with OS and the PFI than the low-expression group [Fig. 7D, H]. These data demonstrated that TBC1D10C is a protective marker. We next calculated the differential expression of all four genes in tumor tissue from the TCGA database and normal tissue from the GTEx database. These data indicated that all four genes were more highly expressed in tumor tissues than in normal tissue [Fig. 8].
Enrichment pathway analysis of TBC1D10C.
Considering the differential prognosis for differential expression of TBC1D10C, we made a gene list to rank the correlation between TBC1D10C and other genes. GSEA indicated that it primarily correlated with, and thereby might be involved in, the allograft rejection, interferon gamma response, inflammatory response, K-Ras signaling, TNFA signaling via the NFKB, and complement pathways [Fig. 9A]. The top 6 of high enrichment score pathway was displayed in Fig. 9B.
Tumor-infiltrating immune cells and immune checkpoints correlation with TBC1D10C in BRCA patients.
We next evaluated the relationship of TBC1D10C expression levels with tumor-infiltrating immune cells in the breast cancer microenvironment using ssGSEA [Fig. 10A]. Only 3 of 28 immune cell types (central memory CD8 T cells, memory B cells, and neutrophils) that were differentially expressed had no statistical significance. The tumor-infiltrating immune cells, e.g., activated CD8 T cells and activated B cells, had a robust correlation with TBC1D10C (p-value < 0.05) [Fig. 10B].
Immune checkpoint blockade therapies are emerging and effective strategies for treating cancer(24). We next explored the associations between TBC1D10C and immune checkpoints including PD1, PD-L1, TIGIT, CTLA-4, TIM-3, and LAG-3. The chord chart indicated that TBC1D10C had a robust positive correlation with PD1, TIGIT, and CTLA-4 [Fig. 10C].
We further analyzed the relationship between the clinical characteristics of BRCA patients and TBC1D10C expression levels. These data indicated that the differential expression of TBC1D10C had a statistically significant difference according to the PAM50 subtype, M stage, N stage, stromal score, and immune score [Table 1].
Table 1
clinical characteristic of patients between high and low expression of TBC1D10C
| | TBC1D10C | |
| Overall | High-expression | Low-expression | P value |
n | 898 | 449 | 449 | |
ER_Status (%) | | | | 0.065 |
Negative | 196 (21.8) | 112 (24.9) | 84 (18.7) | |
Positive | 658 (73.3) | 314 (69.9) | 344 (76.6) | |
unknow | 44 (4.9) | 23 (5.1) | 21 (4.7) | |
PR_Status (%) | | | | 0.954 |
Negative | 280 (31.2) | 142 (31.6) | 138 (30.7) | |
Positive | 571 (63.6) | 284 (63.3) | 287 (63.9) | |
unknow | 47 (5.2) | 23 (5.1) | 24 (5.3) | |
HER2_Status (%) | | | | 0.007 |
Negative | 674 (75.1) | 323 (71.9) | 351 (78.2) | |
Positive | 122 (13.6) | 60 (13.4) | 62 (13.8) | |
unknow | 102 (11.4) | 66 (14.7) | 36 (8.0) | |
PAM50Subtype (%) | | | | < 0.001 |
Basal | 118 (13.1) | 61 (13.6) | 57 (12.7) | |
Her2 | 62 (6.9) | 30 (6.7) | 32 (7.1) | |
LumA | 362 (40.3) | 160 (35.6) | 202 (45.0) | |
LumB | 166 (18.5) | 59 (13.1) | 107 (23.8) | |
unknow | 190 (21.2) | 139 (31.0) | 51 (11.4) | |
M stage (%) | | | | 0.004 |
M0 | 882(98.2) | 447(99.6) | 435 (96.9) | |
M1 | 16 (1.8) | 2 (0.4) | 14 (3.1) | |
N stage (%) | | | | 0.013 |
N0 | 413 (46.0) | 212 (47.2) | 201 (44.8) | |
N1 | 304 (33.9) | 142 (31.6) | 162 (36.1) | |
N2 | 119 (13.3) | 53 (11.8) | 66 (14.7) | |
N3 | 62 (6.9) | 42 (9.4) | 20 (4.5) | |
T stage (%) | | | | 0.09 |
T1 | 236 (26.3) | 128 (28.5) | 108 (24.1) | |
T2 | 532 (59.2) | 259 (57.7) | 273 (60.8) | |
T3 | 101 (11.2) | 53 (11.8) | 48 (10.7) | |
T4 | 29 (3.2) | 9 (2.0) | 20 (4.5) | |
TNM stage (%) | | | | 0.051 |
I | 157 (17.5) | 80 (17.8) | 77 (17.1) | |
II | 526 (58.6) | 269 (59.9) | 257 (57.2) | |
III | 201 (22.4) | 98 (21.8) | 103 (22.9) | |
IV | 14 (1.6) | 2 (0.4) | 12 (2.7) | |
age [IQR] | 58.00 [48.00, 67.00] | 58.00 [49.00, 66.00] | 58.00 [48.00, 67.00] | 0.92 |
TMB [IQR] | 0.97 [0.63, 1.78] | 0.93 [0.58, 1.89] | 1.01 [0.66, 1.71] | 0.412 |
stromal_score [IQR] | 414.52 [-163.90, 838.41] | 554.52 [90.80, 967.12] | 204.14 [-339.70, 700.80] | < 0.001 |
immune_score [IQR] | 121.41 [-362.02, 691.89] | 638.86 [170.02, 1224.44] | -287.32 [-614.74, 72.88] | < 0.001 |
ER, estrogen receptor; PR, progesterone receptor; HER2, human epithelial growth factor receptor 2. IQR, interquartile range. TMB, tumor mutational burden. n, sample numbles. |