3.1 The identification and validation of mrDEGs in BC
The flow diagram was shown in Figure 1. GSE128610 and GSE72319 datasets were analyzed online by GEO2R. The mrDEGs in BC were identified based on the cut-off criteria of |logFC|>=1, P<0.05 and P. adjust<0.05. We found out 1756 up-regulated and 3225 down-regulated mrDEGs in GSE128610. Subsequently, 251 up-regulated and 1162 down-regulated mrDEGs in GSE128610 were validated by GSE72319. After initial validation, they were further verified by TCGA data, and then 23 up-regulated and 71 down-regulated mrDEGs were finally identified (Table 1, Figure 2).
3.2 GO function enrichment and KEGG pathway map for mrDEGs
GO enrichment analysis of 94 validated mrDEGs was performed by STRING. The results were presented in Table 2 and Table S1. The mrDEGs in BC were shown to play roles in cancer-related biological processes, such as neural crest cell migration involved in autonomic nervous system development, regulation of cell migration, cell surface receptor signaling pathway and cell differentiation. However, no significant molecular function or cellular component was observed.
The mrDEGs were mapped in KEGG pathway. They were suggested to participate in the following cancer-related pathways: PI3K-ALT pathway, TGF-beta pathway, evading apoptosis, and resistance to chemotherapy (Figure 3).
3.3 mrDEGs-related PPI network construction and Modeling analysis
PPI network of 94 mrDEGs in BC was constructed by STRING database, with a total of 94 nodes and 60 edges (Figure 4a). Three models were found to meet the criteria of score>=3 and nodes>=3 by MOCODE plug-in of Cytoscape software (Figure 4b). GO and KEGG enrichment analyses were carried out for the three models respectively (Table 3 and Table S2). The results showed that model 1 was mainly associated with the structure and function of nerve cells, model 2 was involved in extracellular matrix, and model 3 took parts in tissue development and gene expression.
3.4 Selection and tissues validation of hub mrDEGs
We screened out 9 hub mrDEGs according to MCC>=6 by cytoHubba plug-in of Cytoscape. Among them, up-regulated hub mrDEGs contained FN1, BGN, EFNA3, COL5A2 and SEMA3F, and down-regulated hub mrDEGs comprised RHOQ, SEMA3A, NRP1 and DDR2 (Table 4). Consistent results were obtained after we validated and visualized the expression of 9 hub mrDEGs by Ualcan with tissue samples (Figure 5).
3.5 Analysis of TF-hub mrDEGs network
In order to explore the potential regulatory relationships of hub mrDEGs, we predicted the TF targeting hub mrDEGs with ENCODE database. The result demonstrated that 8 hub mrDEGs were matched excepting for COL5A2. The Cytoscape software was applied to visualize the TF-hub mrDEGs network, and 167 associations between 8 hub mrDEGs and 121 TF were predicted (Figure 6). As is shown, MAZ could regulate 5 hub mrDEGs (eg. EFNA3, SEMA3A and SEMA3F), and SP2, MXD4, KLF9, KLF16, HDGF and ARID4B regulated 3 hub mrDEGs.
3.6 Analysis of miRNA-hub mrDEGs network
Subsequently, miRNA-hub mrDEGs pairs of 9 hub mrDEGs were performed with TarBase and miRTarBase. Finally, only 6 hub mrDEGs were mapped excepting for BGN, SEMA3F and SEMA3A, next 31 associations between 25 miRNA and 6 hub mrDEGs were obtained by Cytoscape software (Figure 7). We found that hsa-mir-21-5p could interact with 3 hub mrDEGs (eg. COL5A2 and DDR2).
3.7 Construction of TF-miRNA-hub mrDEGs network analysis
The hub mrDEGs (FN1, EFNA3, NRP1, RHOQ and DDR2) which coregulated by TF and miRNA were selected and their interactive regulators were extracted, and then building a TF-miRNA-hub mrDEGs network by Cytoscape (Figure 8). A total of 5 hub mrDEGs, 21 miRNA and 117 TF were included in the TF-miRNA-hub mrDEGs network. Next, we analyzed the interactive results of TF-mrDEGs and miRNA-hub network respectively (Table 5). We found that MAZ, HDGF and SP2 could regulate 3 hub mrDEGs. Simultaneously, hsa-mir-21-5p, hsa-mir-1-3p, hsa-mir-218-5p, hsa-mir-26a-5p, and hsa-mir-335-5p regulated 2 hub mrDEGs. In addition, FN1, EFNA3, and NRP1 were the highest degree score in interaction network.
3.8 Survival analysis of coregulated mrDEGs by TF and miRNA
Kaplan Meier plotter was adopted to analyze overall survival curve of with the expression of 5 coregulated hub mrDEGs by TF and miRNA in BC. Two hub mrDEGs with statistical significance were identified in the survival analysis (P<0.05, n=1402), including FN1 (HR = 1.28 (1.03 − 1.59), P = 0.023) and DDR2 (HR = 0.77 (0.62 − 0.96), P = 0.017). Therefore, the up-regulated FN1 and down-regulated DDR2 might confer to poor BC prognosis (Figure 9).