Screening M2-like TAMs-related genes
To elucidate the relationship between macrophages and prognosis of BRCA, the contents of M1 and M2 macrophages in TCGA-BRCA samples were calculated by CIBERSORT algorithm. Next, BRCA patients were divided into high-content and low-content M1 or M2 macrophages group. Kaplan-Meier survival analyses showed some difference between high-content and low-content M1 macrophage groups (Fig. 1A, P-value = 0.032). Furthermore, significant difference in survival has shown between high-content and low-content M2 macrophages group, and the patients in the high-content M2 macrophages had worse prognosis than that of the low-content M2 macrophages (Fig. 1B, P-value < 0.001). This result indicated that the effects of M1 and M2 macrophages on prognosis were quite opposite, and M2 macrophages had a more pivotal role in BRCA than M1 macrophages. GSE58812 datasets confirmed the above result (Fig. 1C, D). Based on this observation, WGCNA was performed to identify the key module related to M2 macrophages. Figure 1E demonstrates that the green module was positively and negatively correlated with M2 macrophages and M1 macrophages, respectively. And 3,127 genes in the green module were selected for downstream analyses. Then, 2,098 DEGs (935 up-regulated and 1,163 down-regulated genes) were screened out from the breast cancer tissues compared to the normal tissues in GSE42568. Meanwhile, 2,753 DEGs (1,078 up-regulated and 1,675 down-regulated genes) were identified from TCGA-BRCA. Lastly, we searched for co-expressed genes between WGCNA-derived M2-like TAMs-related genes and up-regulated DEGs associated with BRCA, and eventually screened out 85 overlapping genes, which may play an important role in the development and progression of breast cancer (Fig. 1F).
Go enrichment and KEGG pathway and PPI network analysis
GO and KEGG analyses were conducted to further explore the potential roles of intersecting hub genes. GO analysis showed that the overlapping 85 genes mainly affect the cellular component of bicellular tight junction, tight junction and apical junction complex (Fig. 2A). KEGG pathway analysis showed that these 85 genes mainly affect tight junction, Rap1 signaling pathway and PI3K-Akt signaling pathway (Fig. 2B). It is well known that the Rap1 signaling pathway and PI3K-Akt signaling pathway are closely associated with proliferation, invasion and metastasis of cancer, while the relationship between tight junction and tumor has been rarely reported. Breast cancer is a complex and heterogeneous disease, and proper adhesion between adjacent epithelial cells of the breast is essential for the normal structure and function of the epithelial tissue. There is growing evidence that dysregulation of cell-cell adhesion is associated with many cancers. There is a study that suggests that tight junction may be involved in the development or progression of breast cancer [20]. The PPI network of candidate hub genes was constructed by using the STRING online tool. Subsequently, the top 10 up-regulated genes in breast cancer were then visually analyzed by using "CytoNCA" toolkit and degree algorithm. Briefly, FOXA1, AGR2, MUC1, ERBB3, TFF1, EZR, GATA3, FGFR3, SPDEF and AGR3 were sorted out (Fig. 2C). The deeper the red color, the higher the relevance score, and the more significant the importance in the BRCA. It is clear that FOXA1, AGR2, MUC1 and ERBB3 are the most important genes.
Construction of nomogram and ROC diagnostic curve
We constructed a nomogram model to predict the risk of BRCA. Therefore, our nomogram performs well in BRCA prediction, the higher the scores, the higher the risk of disease. As shown in Fig. 2D, the patient has 100 total points of these four hub genes, the risk of breast cancer is up to more than 90%. Figure 2E verifies the reliability of nomogram. Subsequently, we calculated ROC curves for the four hub genes (FOXA1, AGR2, MUC1 and ERBB3) to assess reliability of prediction, and the AUC values of FOXA1, MUC1, ERBB3 and AGR2 are 0.877, 0.907, 0.930 and 0.865, respectively. It demonstrated that these four hub genes could differentiate breast cancer from normal tissues (Fig. 2F).
Functional analysis of M2-like TAMs-related hub genes
According to the correlation analysis between four hub genes and M2-like macrophages, three of the four hub genes were closely related with M2-like macrophages. Compared with other genes, FOXA1 and ERBB3 were more closely associated with M2-like macrophages, among which FOXA1 and ERBB3 had the highest correlation (Fig. 3A). In TCGA datasets, high FOXA1 expression correlated with shorter overall survival (OS) (Fig. 3B). But the high expression of AGR2 group showed a better survival rate than the low expression of AGR2 group (Fig. 3C). The expressed differences of four hub genes between normal tissues and breast tumor tissues were analyzed based on TCGA-BRCA, and the four hub genes were all highly expressed in breast tumor tissues compared with normal tissues (Fig. 3D). GSEA analysis showed that FOXA1 and ERBB3 were mainly related to tumor-related pathways, such as complement and coagulation cascades [21], colorectal cancer, mtor signaling pathway and prostate cancer (Table 1). Complement and coagulation cascade are important components of the human immune system. Platelets in TME can promote tumor cell invasion and metastasis by activating the coagulation cascade. And then, the relationship between the expression levels of four hub genes and the tumor immune microenvironment was analyzed. CIBERSORT algorithm analysis showed that four hub genes were positively correlated with the infiltration of immune cells such as M2 macrophages, resting Mast cells, resting CD4+ T cells, monocytes and native B cells (Fig. 3E-H). In addition, FOXA1, ERBB3, AGR2 had a higher positive correlation with M2 macrophages, which was consistent with the correlation analysis results in Fig. 3A. The expression of these four hub genes was positively correlated with M2-like macrophages infiltration in BRCA, but it was completely opposite to M1-like macrophages infiltration in BRCA.
Table 1
Gene set enrichment analysis (GSEA) of four hub genes
Gene | KEGG Pathway | Enrichment Score | P-value |
FOXA1 | GLYCOSYLPHOSPHATIDYLINOSITOL_GPI_ANCHOR_BIOSYNTHESIS | 0.64890 | 0.00064 |
COMPLEMENT_AND_COAGULATION_CASCADES | 0.46633 | 0.00027 |
VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION | 0.51269 | 0.00147 |
NUCLEOTIDE_EXCISION_REPAIR | 0.49271 | 0.00267 |
PROTEIN_EXPORT | 0.57175 | 0.00729 |
AGR2 | COMPLEMENT_AND_COAGULATION_CASCADES | 0.44045 | 0.00459 |
GLYCOSYLPHOSPHATIDYLINOSITOL_GPI_ANCHOR_BIOSYNTHESIS | 0.59598 | 0.02040 |
PRIMARY_IMMUNODEFICIENCY | -0.60826 | 0.02514 |
SPLICEOSOME | 0.32993 | 0.02020 |
ERBB3 | PROSTATE_CANCER | 0.42079 | 0.00012 |
MTOR_SIGNALING_PATHWAY | 0.50957 | 0.00023 |
GLYCOSYLPHOSPHATIDYLINOSITOL_GPI_ANCHOR_BIOSYNTHESIS | 0.67080 | 0.00039 |
COLORECTAL_CANCER | 0.42077 | 0.00204 |
RNA_DEGRADATION | 0.44905 | 0.00061 |
MUC1 | VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION | 0.54677 | 0.00018 |
FATTY_ACID_METABOLISM | 0.51296 | 0.00144 |
CIRCADIAN_RHYTHM_MAMMAL | 0.69052 | 0.00805 |
GLYCOSYLPHOSPHATIDYLINOSITOL_GPI_ANCHOR_BIOSYNTHESIS | 0.51112 | 0.01853 |
GALACTOSE_METABOLISM | 0.48660 | 0.03024 |
In pan-cancer analysis, four hub genes were found to be highly expressed in most kinds of tumor tissues, and FOXA1, ERBB3 had higher expression in breast tumors than other tumors particularly (Fig. 4A-D). Moreover, it was found that the age of patients over 65 years old were positively associated with the expression level of four hub genes from the clinical analysis of BRCA (Fig. 4E-H).
Hub genes causally associated with the risk of BRCA
We assessed the causal association between four hub genes and breast cancer. The SNPs characteristics of four hub genes were strong instrumental variables (P value < 5*10− 6). Using the fixed effects IVW method, we found that four hub genes were all associated with the risk of breast cancer (Table 2). β-value is greater than zero indicates that exposure factor is a risk factor, and the risk of morbidity increases with the increase of exposure factor. Except for AGR2, the expressions of FOXA1, ERBB3 and MUC1 were positively correlated with the occurrence of breast tumors (β-value > 0), this result was consistent with the prognostic analysis of AGR2 (Figs. 3C). Obvious causal association between FOXA1 and breast cancer was observed in Figs. 5A. The combined effect value of all SNPs is greater than zero, indicating that FOXA1 is a risk factor, and the risk of breast cancer increases with the increase of FOXA1 expression (Figs. 5B). The weighted median and weighted mode method also confirmed the result of fixed effects IVW method. The causal effect of FOXA1 funnel plot was roughly symmetrical along the symmetry axis (Figs. 5C), and no horizontal pleiotropy is observed at the intercept of MR Egger regression, which further indicates that pleiotropy does not bias the causal effect. In the leave-one-out plot, the MR analysis was performed again on the remaining SNPs after removing each SNP, the comprehensive effect value was consistent with the main effect value, indicating that the calculation results of all SNPs made causality significant (Figs. 5D).
Table 2
The results of mendelian randomization (MR) analysis
Gene | Exposure | Outcome | Method | nSNPs | β-value | P-value |
FOXA1 | prot-a-2611 | Breast cancer id:ebi-a-GCST007236 | MR Egger | 26 | 0.1462 | 0.9402 |
Weighted median | 26 | 0.9526 | 5.54E-11 |
Inverse variance weighted (fixed effects) | 26 | 0.6474 | 5.44E-245 |
Simple mode | 26 | 0.4266 | 0.1621 |
Weighted mode | 26 | 0.5861 | 9.18E-07 |
MUC1 | prot-a-1967 | Breast cancer id:ebi-a-GCST007236 | MR Egger | 16 | 4.6470 | 0.0002 |
Weighted median | 16 | 1.8530 | 4.25E-18 |
Inverse variance weighted (fixed effects) | 16 | 0.7446 | 0 |
Simple mode | 16 | 1.6230 | 0.0040 |
Weighted mode | 16 | 1.6600 | 6.72E-10 |
AGR2 | prot-a-58 | Breast cancer id:ebi-a-GCST007236 | MR Egger | 10 | -1.3710 | 0.6640 |
Weighted median | 10 | -0.6746 | 7.32E-05 |
Inverse variance weighted (fixed effects) | 10 | -0.4915 | 6.48E-53 |
Simple mode | 10 | -1.6350 | 0.0004 |
Weighted mode | 10 | -0.6215 | 0.0061 |
ERBB3 | prot-c-2617_56_35 | Breast cancer id:ebi-a-GCST007236 | MR Egger | 3 | 1.4770 | 0.9365 |
Weighted median | 3 | 0.3960 | 2.55E-06 |
Inverse variance weighted (fixed effects) | 3 | 0.03712 | 1.93E-78 |
Simple mode | 3 | 0.3527 | 0.0267 |
Weighted mode | 3 | 0.3025 | 0.0229 |
Drug Sensitivity Analysis
To investigate the sensitization or resistance of four hub genes to drugs during chemotherapy, we integrated drug IC50 and gene expression profile data from GDSC BRCA cell lines (Table 3). Through this method, the drugs with the highest degree of correlation with hub genes can be screened out. The results showed that the IC50 of A-443654, Lapatinib, KIN001-102, CP724714 and ZSTK474 were negatively correlated with FOXA1 expression, and the IC50 of Pyrimethamine and BX-795 were positively correlated with FOXA1 expression according by Pearson’s correlation analysis (Fig. 6A). The IC50 of Lapatinib, A-770041, Sorafenib, Roscovitine, A-443654, MS-275 and Paclitaxel was negatively correlated with the expression of ERBB3, indicating that high ERBB3 expression has the sensitizing effect on these drugs (Fig. 6B). The IC50 of Saracatinib, Imatinib, BMS-509744 and Lapatinib was negatively correlated with MUC1 expression, and Olaparib, YK 4-279 was positively correlated with MUC1 expression (Fig. 6C). Interestingly, the expression of FOXA1, ERBB3 and MUC1 all increased the sensitivity of Lapatinib to BRCA patients. In 2007, FDA approved Lapatinib as a common treatment for HER2-positive advanced breast cancer. Whether the combined expression of FOXA1, ERBB3 and MUC1 could be considered as a potential opportunity for the using of Lapatinb in clinical breast cancer patients. Moreover, the expression of AGR2 is closely related to the resistance of TAE684, XMD8-85, AZ628 and Pyrimethamine (Fig. 6D). In other words, the expression of AGR2 is more likely to lead to drug resistance.
Table 3
The results of drug sensitivity analysis
Gene | Drug | Pearson Correlation Coefficient | P-value |
FOXA1 | A-443654 | -0.502 | 0.00736 |
Lapatinib | -0.488 | 0.00883 |
KIN001-102 | -0.419 | 0.02051 |
CP724714 | -0.414 | 0.02173 |
ZSTK474 | -0.391 | 0.02821 |
Pyrimethamine | 0.634 | 0.02205 |
BX-795 | 0.616 | 0.02705 |
ERBB3 | Lapatinib | -0.633 | 0.00420 |
A-770041 | -0.589 | 0.00724 |
Sorafenib | -0.530 | 0.01437 |
Roscovitine | -0.493 | 0.02152 |
A-443654 | -0.490 | 0.02222 |
MS-275 | -0.474 | 0.02627 |
Paclitaxel | -0.473 | 0.02655 |
MUC1 | Saracatinib | -0.708 | 2.11E-05 |
Imatinib | -0.402 | 0.01027 |
BMS-509744 | -0.347 | 0.02361 |
Lapatinib | -0.334 | 0.02839 |
Olaparib | 0.491 | 0.01859 |
YK 4-279 | 0.481 | 0.02153 |
AGR2 | TAE684 | 0.709 | 0.00735 |
XMD8-85 | 0.68 | 0.01109 |
KIN001-102 | -0.355 | 0.01299 |
CP724714 | -0.314 | 0.02230 |
AZ628 | 0.625 | 0.02303 |
Pyrimethamine | 0.613 | 0.02678 |
ZSTK474 | -0.291 | 0.02974 |
Afatinib | -0.291 | 0.02974 |