Differences in gene expression between high and low stromal scores
The stromal score was calculated to predict the non-tumor cell infiltration [21]. Based on the stromal score, the samples were split into high and low stromal score groups, and a total of 478 genes were significantly differentially expressed in the two groups, including 104 upregulated and 374 downregulated genes (Figure 2A). These genes were considered to be associated with the stroma in tumor tissue.
Differences in gene expression between high and low immune scores
The immune score refers to the infiltration of immune cells in tumor tissue, which is considered an available indicator of prognosis in tumors [30]. Based on the immune score, the samples were split into high and low immune score groups; 796 genes were found to be significantly differentially expressed in the 2 groups (Figure 2B). Of these, 503 genes were upregulated whereas 293 genes were downregulated, suggesting that these genes might be implicated in the immune status in the tumor microenvironment.
Overlapping genes between the stroma score-related and immune score-related genes
The overlapping genes between the stroma score-related and immune score-related genes were identified by VENN analysis; 167 genes were obtained, including 58 upregulated and 109 downregulated genes (Figure 2C). The functions of these overlapping genes were further investigated, and the downregulated genes were found to be significantly implicated in one KEGG pathway (hsa05146, amoebiasis) and five biological process terms, for example: GO:0006885~regulation of pH (Table 1).
Immune cell infiltration in BRCA
The infiltration abundance of six immune cell types in BRCA was analyzed using the Cibersort algorithm. From the bar charts of immune cell subset proportions (Figure 3), CD8+ T cells and activated memory CD4+ T cells accounted for a large proportion of immune cell infiltration in BRCA. Therefore, CD8+ and CD4 T+ cells were selected in the following analysis.
Identification of CD4+ and CD8+ T cell-related genes
The infiltrating abundance of CD4+ and CD8+ T cells is related to prognosis in BRCA [13]. The infiltrating abundance of CD4+ and CD8+ T cells was calculated using the Cibersort algorithm. Then, the immune cell-related genes were identified using the Spearman correlation test between overlapping genes and the infiltrating abundance of immune cells. A total of 39 CD4+ T cell-related genes and 78 CD8+ T cell-related genes were identified. These genes were regarded as crucial genes involved in CD4+ and CD8+ T cell infiltration.
PPI network of CD4+ and CD8+ T cell-related genes
Proteins and their functional interactions form the backbone of cellular machinery, and connectivity networks are conducive to fully understand biological phenomena [31]. For CD4+ T cell-related genes, the PPI network contained 13 genes and 11 interactions (Figure 4A). Of the 13 genes, 5 genes were upregulated whereas 8 genes were downregulated. For CD8+ T cell-related genes, the PPI network consisted of 56 genes (16 upregulated and 40 downregulated) and 76 interactions (Figure 4B).
ceRNA network of CD4+ and CD8+ T cell-related genes
For CD4+ T cell-related genes, 17 miRNA-mRNA interactions (Supplemental Figure 1A), and 13 lncRNA-mRNA interactions (Supplemental Table 1) were finally predicted. The ceRNA network is shown in Figure 5A, and contains 2 lncRNAs, 12 miRNAs, and 5 mRNAs, consisting of 25 ceRNA regulatory axes. For example: the chr22-38_28785274-29006793.1–miR-34a/c-5p–CAPN6 axis was identified.
For CD8+ T cell-related genes, a total of 51 miRNA-mRNA pairs (Supplemental Figure 1B), and 43 lncRNA-mRNA pairs were predicted (Supplemental Table 2). The ceRNA network contained 10 lncRNAs, 32 miRNAs, and 14 mRNAs, consisting of 81 ceRNA regulatory axes (Figure 5B). For example: the chr22-38_28785274-29006793.1–miR-494-3p–SLC9A7 axis was identified.
Breast neoplasm-related chemicals target mRNAs in the ceRNA network
The etiology of many diseases involves the interactions between environmental chemicals and genes that regulate physiological processes [32]. The CTD provides information about chemical–gene/protein-disease relationships [32]. In this study, the chemical–gene interactions involving breast neoplasms were identified from the CTD, and were filtered using the mRNAs in the ceRNA network. A total of 31 chemicals were found to interact with the five genes in the CD4+ T cell-related ceRNA network (Figure 6A, Supplemental Table 3). For example, atrazine (CAS, 1912-24-9) was predicted to result in the increased expression of CAPN6 mRNA.
Similarly, 57 chemicals were found to interact with the 12 genes in the CD8+ T cell-related ceRNA network (Figure 6B, Supplemental Table 4). For example, arsenic (CAS, 7440-38-2) was predicted to affect the methylation of the SLC9A7 gene.
CD4+ and CD8+ T cell-related genes associated with prognosis of BRCA
To explore the prognostic value of CD4+ and CD8+ T cell-related genes, survival analysis was performed. In total, 14 genes were found to be significantly related to prognosis of BRCA patients; of these, eight genes were commonly related to both CD4+ and CD8+ T cells (Figure 7, Table 2), for example, MUC2. Among the 14 genes, 6 genes were found to be specific to CD4+ or CD8+ T cells; METTL5 and MSTN were CD4+ T cell-related genes, whereas NMNAT3, GNAL, ZBED4, and RDH12 were CD8+ T cell-related genes.