3.1 BRCA-related DEG identification
By pooling data from three GEO microarray datasets, we were collectively able to analyze 304 BRCA tumor tissue samples and 32 normal breast tissue control samples. With the aid of the GEO2R tool, we were able to extract 2384, 1176, and 260 BRCA-related DEGs from the GSE45827, GSE42568, and GSE61304 datasets, respectively. Using a Venn diagram approach, we then identified 85 DEGs that were common to all three datasets. Of these DEGs, 54 were downregulated (logFC<0) and 31 were upregulated(logFC>0) in BRCA tumor tissues (Supplementary Table 2 & Figure 1).
3.2 Functional enrichment analyses of identified DEGs
In order to understand the functional roles of these 85 BRCA-associated DEGs, we next conducted GO and KEGG analyses thereof. GO analyses examined the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) with which these DEGs were affiliated. The GO BP analysis revealed the upregulated DEGs to be primarily enriched for processes including mitotic nuclear division, cell division, G2/M transition of mitotic cell cycle, collagen catabolic process, and endodermal cell differentiation, whereas downregulated DEGs were enriched lipid metabolic processes, lipid transport, regulation of cardiac muscle contraction by regulation of the release of sequestered calcium ions, positive regulation of cell proliferation, and positive regulation of cell-matrix adhesion. The GO CC analysis revealed upregulated DEGs to be primarily enriched in the proteinaceous extracellular matrix, midbody, condensed chromosome kinetochore, centriole, and microtubule cytoskeleton, whereas downregulated DEGs were enriched in the extracellular space, extracellular region, lipid particles, proteinaceous extracellular matrix, and the apical plasma membrane. With respect to the GO MF analysis, upregulated DEGs were enriched for functions including collagen binding, heparin binding, ATP binding, and metalloendopeptidase activity, with downregulated DEGs being enriched lipid binding, heparin binding, lipoprotein particle binding, transporter activity, zinc-dependent alcohol dehydrogenase activity (Table 1).
We further conducted a KEGG pathway analysis of these DEGs. This approach revealed the upregulated DEGs to primarily be associated with the p53 signaling pathway, ECM-receptor interactions, and oocyte meiosis, whereas downregulated DEGs were primarily associated with the PPAR signaling pathway, tyrosine metabolism, cytochrome P450 drug metabolism, and protein digestion and absorption. KEGG pathway enrichment analysis results for upregulated and downregulated DEGs are shown in Figure 2.
Table 1 GO enrichment analysis of BRCA-related DEGs (The terms with the top 5 P-values are shown)
Expression
|
Category
|
Term
|
Count
|
P-value
|
Fold enrichment
|
FDR
|
Up-regulated
|
GOTERM_BP_DIRECT
|
GO:0007067~mitotic nuclear division
|
11
|
9.55E-12
|
24.02601457
|
1.31E-08
|
|
GOTERM_BP_DIRECT
|
GO:0051301~cell division
|
9
|
1.29E-07
|
13.92884793
|
1.77E-04
|
|
GOTERM_BP_DIRECT
|
GO:0000086~G2/M transition of mitotic cell cycle
|
6
|
4.06E-06
|
23.72309866
|
0.005568755
|
|
GOTERM_BP_DIRECT
|
GO:0030574~collagen catabolic process
|
5
|
4.88E-06
|
42.31854839
|
0.006689601
|
|
GOTERM_BP_DIRECT
|
GO:0035987~endodermal cell differentiation
|
4
|
1.46E-05
|
80.24850657
|
0.020038146
|
|
GOTERM_CC_DIRECT
|
GO:0005578~proteinaceous extracellular matrix
|
8
|
2.10E-07
|
17.5483871
|
2.23E-04
|
|
GOTERM_CC_DIRECT
|
GO:0030496~midbody
|
6
|
2.03E-06
|
27.34283571
|
0.002166272
|
|
GOTERM_CC_DIRECT
|
GO:0000777~condensed chromosome kinetochore
|
4
|
3.89E-04
|
27.02855024
|
0.413450196
|
|
GOTERM_CC_DIRECT
|
GO:0005814~centriole
|
4
|
8.34E-04
|
20.80959178
|
0.885431205
|
|
GOTERM_CC_DIRECT
|
GO:0015630~microtubule cytoskeleton
|
4
|
0.001454401
|
17.16411585
|
1.539162614
|
|
GOTERM_MF_DIRECT
|
GO:0005518~collagen binding
|
4
|
1.62E-04
|
36.30322581
|
0.177777537
|
|
GOTERM_MF_DIRECT
|
GO:0008201~heparin binding
|
4
|
0.002811565
|
13.61370968
|
3.047438273
|
|
GOTERM_MF_DIRECT
|
GO:0005524~ATP binding
|
8
|
0.013957193
|
2.913971302
|
14.3155518
|
|
GOTERM_MF_DIRECT
|
GO:0004222~metalloendopeptidase activity
|
3
|
0.017098388
|
14.45703683
|
17.26865537
|
Down-regulated
|
GOTERM_BP_DIRECT
|
GO:0006629~lipid metabolic process
|
6
|
5.42E-05
|
14.26072187
|
0.080640059
|
|
GOTERM_BP_DIRECT
|
GO:0006869~lipid transport
|
4
|
0.00103261
|
19.63976608
|
1.52497516
|
|
GOTERM_BP_DIRECT
GOTERM_BP_DIRECT
|
GO:0010881~regulation of cardiac muscle contraction by regulation of the release of sequestered calcium ion
GO:0008284~positive regulation of cell proliferation
|
3
7
|
0.001115403
0.001276576
|
58.91929825
5.605340963
|
1.646301724
1.882087273
|
|
GOTERM_BP_DIRECT
|
GO:0001954~positive regulation of cell-matrix adhesion
|
3
|
0.001499269
|
50.88484848
|
2.207004361
|
|
GOTERM_CC_DIRECT
|
GO:0005615~extracellular space
|
14
|
4.96E-05
|
3.713932195
|
0.055590194
|
|
GOTERM_CC_DIRECT
|
GO:0005576~extracellular region
|
14
|
3.01E-04
|
3.107246377
|
0.33641381
|
|
GOTERM_CC_DIRECT
|
GO:0005811~lipid particle
|
4
|
7.87E-04
|
21.65656566
|
0.878449732
|
|
GOTERM_CC_DIRECT
|
GO:0005578~proteinaceous extracellular matrix
|
6
|
8.18E-04
|
8
|
0.912134466
|
|
GOTERM_CC_DIRECT
|
GO:0016324~apical plasma membrane
|
6
|
0.001180873
|
7.367697595
|
1.314786411
|
|
GOTERM_MF_DIRECT
|
GO:0008289~lipid binding
|
4
|
0.005831972
|
10.64711447
|
6.708154954
|
|
GOTERM_MF_DIRECT
|
GO:0008201~heparin binding
|
4
|
0.00684227
|
10.04821429
|
7.827449644
|
|
GOTERM_MF_DIRECT
|
GO:0071813~lipoprotein particle binding
|
2
|
0.012086409
|
160.7714286
|
13.44231459
|
|
GOTERM_MF_DIRECT
|
GO:0005215~transporter activity
|
4
|
0.012887447
|
7.958981612
|
14.27185096
|
|
GOTERM_MF_DIRECT
|
GO:0004024~alcohol dehydrogenase activity, zinc-dependent
|
2
|
0.014486531
|
133.9761905
|
15.90609305
|
Abbreviations: BP, biological process; CC, cellular component; Count, the number of enriched genes in each term; FDR, false discovery rate; MF, molecular function; BRCA, breast cancer.
3.3 PPI network construction and analysis
We next generated a PPI network in an effort to better understand the relationships between the 85 BRCA-related DEGs identified above. To that end, we leveraged the STRING database to construct a network that successfully incorporated 63 of these 85 DEGs (28 and 35 up- and down-regulated, respectively). The final network contained 63 nodes and 215 edges (Figure 3A). Using the MCODE plugin, we then identified 16 nodes corresponding to hub genes within this network using filtering degree≥14 as a criterion for hub identification (Figure 3B).
3.4 Assessment of the prognostic relevance of core genes
Next, we tested the prognostic relevance of these 16 core genes using the Kaplan-Meier Plotter database. This analysis revealed that the upregulation of 14 of these genes was significantly associated with poorer patient survival (P <0.05, Supplementary Table 3 & Figure 4). We then further validated the expression of these 14 prognostic hub genes using the GEPIA tool, confirming that all 14 of these genes were expressed at higher levels in BRCA tumor tissue samples relative to normal breast tissue (P<0.01, Figure 5). Importantly, patients with higher expression of these genes had poorer outcomes than did patients with lower expression thereof. The gene symbols, full names, and functions of these prognosis-related genes are shown in Table 2.
Tab 2 Functions of key prognosis-related genes
Gene symbol
|
Full name
|
Function
|
TOP2A
|
Topoisomerase (DNA) II α
|
TOP2A can be therapeutic targets for anticancer, and mutations in this gene is associated with the development of drug resistance[34].
|
ASPM
|
Abnormal spindle microcephaly associated
|
ASPM plays a role in the regulation of mitotic spindle, with a priority role in regulating neurogenesis. Overexpression of ASPM links with tumor progression and poor prognosis in various tumors, including breast cancer[35].Furthermore, it has been reported that ASPM promotes tumor development by regulating Wnt-β-catenin signaling[36, 37].
|
NEK2
|
Never in Mitosis (NIMA) Related Kinase 2
|
NEK2 plays a critical role in multiple aspects of mitotic processes, such as centrosome duplication and separation, microtubule stabilization, kinetochore attachment, and spindle assembly checkpoint[38].NKK2 is frequently overexpressed in a variety of human malignancies, has been implicated in tumorigenesis, tumor progression and drug resistance[39, 40].
|
RRM2
|
ribonucleotide reductase M2
|
RRM2, a rate-limiting enzyme for DNA synthesis and repair, overexpression exhibited enhanced cellular invasiveness, serving as a potential prognostic marker and a possible therapeutic target for cancer treatment.
|
UBE2T
|
Ubiquitin-conjugating enzyme E2T
|
UBE2T plays a key role in cellular processes, for example, signal transduction, cell cycle control and tumorigenesis[41].
|
CEP55
|
Centrosomal protein 55
|
CEP55 localizes to the centrosome throughout the cell cycle, and to regulate both the mitotic spindle and microtubule organization[42].
|
ANLN
|
Anillin, actin binding protein
|
ANLN plays an important role in cytokinesis, highly express in numerous cancer types[43, 44].
|
DLGAP5
|
discs large associated protein 5
|
DLGAP5 is a regulator of cell cycle involved in carcinogenesis. Elevated DLGAP5 expression correlates with poor prognosis in breast cancer patients[45].
|
CCNB2
|
Cyclin B2
|
Cyclin B2 may regulate the cell cycle with transforming growth factor beta-mediated. Aberrant expression of CCNB2 is correlated with tumor invasion, metastasis, poor prognosis of various human cancers[46].
|
BIRC5
|
Survivin, Baculoviral IAP repeat containing 5
|
BIRC5 may inhibit apoptosis and is usually overexpressed in most tumors[47].
|
TPX2
|
Targeting protein for Xenopus kinesin-like protein 2
|
TPX2 is a microtubule-associated protein required
for mitosis and spindle assembly.
|
AURKA
|
Aurora kinase A
|
AURKA, a serine/threonine kinase, has important functions in regulating cell cycle and mitosis.
|
KIF2C
|
Kinesin family member 2C
|
KIF2C, a crucial mitotic regulator, is highly expressed in many cancer types[48].
|
UBE2C
|
Ubiquitin-conjugating enzyme E2C
|
UBE2C is overexpressed in many human tumor types[49]
|
3.5 Integrated mRNA-miRNA interaction network analysis
Subsequently, Candidate target miRNAs associated with fourteen prognosis-related genes were predicted using miRWalk, TargetScan, and miRDB. Overlapping miRNAs predicted by all three databases were then selected, with a specific focus on miRNAs capable of binding to 3’UTR regulatory regions. An interaction network incorporating these miRNAs and core genes was then constructed using Cytoscape (Figure 6). The network contained 217 nodes and 214 edges. By comparing the targets of core genes, RRM2, AURKA, NEK2, and BIRC5 were found to be potential targets of 53 miRNAs, 41 miRNAs, 30 miRNAs and 28 miRNAs, respectively. Those miRNAs that simultaneously regulated a high number of gene cross-links (≥2) were selected (Supplementary Table 4).
3.6 Construction of a mRNA-target miRNA-TF regulatory network
The selected miRNAs and corresponding genes may play a crucial role in the pathogenesis of BRCA. Hence, RRM2, CEP55, ANLN, UBE2C, NEK2, TOP2A, AURKA, KIF2C, and BIRC5 were considered to be key genes of interest in this context. In order to further understand the function of these genes, a miRNA-target gene-TF regulatory network for these key genes was constructed containing 9 genes, 10 miRNAs, and 3 TFs (Figure 7).
3.7 Verification of candidate molecules expression by RT-PCR and Western Blot
For the purpose of enhancing the credibility of results, the mRNA levels of the three hub genes (AURKA, RRM2 and BIRC5) were measured by RT-PCR and E2F1protein expression was determined with western blot in MDA-MB-231 and MCF-10A cells. As illustrated in Figure 8A, all three genes were significantly upregulated in breast cancer cells (P < 0.01), as predicted by the bioinformatics analysis. E2F1 protein was highly expressed in breast cancer cells as shown in Figure 8B.