CENPL, ISG20L2, LSM4, MRPL3 are four novel hub genes and may serve as diagnostic and prognostic markers in breast cancer

As a highly prevalent disease among women worldwide, breast cancer remains in urgent need of further elucidation its molecular mechanisms to improve the patient outcomes. Identifying hub genes involved in the pathogenesis and progression of breast cancer can potentially help to unveil mechanism and also provide novel diagnostic and prognostic markers. In this study, we integrated multiple bioinformatic methods and RNA in situ detection technology to identify and validate hub genes. EZH2 was recognized as a key gene by PPI network analysis. CENPL, ISG20L2, LSM4, MRPL3 were identified as four novel hub genes through the WGCNA analysis and literate search. Among these, many studies on EZH2 in breast cancer have been reported, but no studies are related to the roles of CENPL, ISG20L2, MRPL3 and LSM4 in breast cancer. These four novel hub genes were up-regulated in tumor tissues and associated with cancer progression. The receiver operating characteristic analysis and Kaplan–Meier survival analysis indicated that these four hub genes are promising candidate genes that can serve as diagnostic and prognostic biomarkers for breast cancer. Moreover, these four newly identified hub genes as aberrant molecules in the maintenance of breast cancer development, their exact functional mechanisms deserve further in-depth study.

Breast cancer is one of the most common malignant tumors that present serious and major threats to female life and health. Although current breast cancer therapeutic methods have been well developed and improved, latest data showed that breast cancer still has a high mortality rate among women worldwide. Thus, there is still an urgent need to explore the potential molecular mechanisms for improving the patient outcomes 1,2 . In the past decade, applications of high-throughput chip and sequencing technologies have resulted in accumulation of a wealth of novel research data resource that can be analyzed by a series of bioinformatic methods, providing a novel approach to explore the molecular mechanism of tumorigenesis and tumor development 3 . Among a wide range of different bioinformatics tools, weighted gene co-expression network analysis (WGCNA) algorithm is the most commonly used method for gene co-expression network research. By constructing coexpression gene modules and associating external information, the key gene modules and potential hub genes can be identified [4][5][6] . In general, hub genes show high connectivity in the gene co-expression network, which often located in the upstream of the gene regulatory network and play a predominant role in the gene network coordination 7,8 . Therefore, identification of potential novel hub genes is of great significance for exploring the mechanism of tumor initiation and progression. Selection of appropriate dataset is an important prerequisite for screening hub genes, and multiple types of datasets that were generated from different platforms are now available in the public database 9 . To explore the best potential of these datasets, it would be of great advantage to integrate them for downstream analysis. To achieve this goal, we used Robust Rank Aggregation (RRA) analysis algorithm for the process of breast cancer datasets. RRA is a reliable bioinformatic method that can remove substantial inter-study variations and statistical analysis difficulties existed in individual studies via integrating the gene expression profiles of different cross-platform datasets 10 . It has been used in various malignant tumor studies, such as in hepatocellular carcinoma, colorectal cancer, lung cancer and thyroid cancer [11][12][13][14] . Heterogeneity is one of the characteristics of tumor cells, which is reflected by different expression patterns of genes at the transcription level 15 . Analyzing and identifying the temporal and spatial heterogeneity information of RNA expression can be of great value to reveal the structural relationship between tissues and cells, as well as to uncover the potential functions of genes in disease state. RNA in situ detection technology can be used for studying the heterogeneity of RNA expression, and under the condition of maintaining tissue and cell morphology integrity, it can obtain the spatial localization and abundance of intracellular RNA at the single cell level 16 .
In this study, we first integrated DEGs from multiple breast cancer datasets based on RRA algorithms, and then identify the key gene of known functions using PPI network analysis. At the same time, WGCNA algorithm was applied to construct a weighted gene co-expression network and screen for potential novel hub genes related to breast cancer. The diagnostic performance and prognostic value of these novel hub genes were evaluated and their possible molecular mechanisms in breast cancer were explored by bioinformatic methods. Finally, we also made full use of RNA in situ detection technology to detect the expression abundance and spatial localization of each hub gene at single cell level, and further analyzing the expression differences and correlations, such that to validate those results from the bioinformatic analysis mentioned above.

Results
Identification of robust DEGs in breast cancer by the RRA analysis. Differentially expressed genes (DEGs) of eight datasets from the GEO database were integrated to perform RRA analysis, and the characteristics for each dataset are shown in Table 1. We used |log2FC|> 1 and p value < 0.05 as screening criteria to obtain the robust DEGs between breast cancer tissues and normal tissues. A total of 512 robust DEGs were identified, containing 202 up-regulated genes and 310 down-regulated genes (Supplementary Table S1). Supplementary  Fig. S1 shows the top 20 most significant up-regulated and down-regulated robust DEGs obtained by RRA methods from these eight different datasets. Of those, COL11A1 (P = 2.47E−19, adjusted P = 6.53E−15, logFC = 2.86) and S100P (P = 1.24E−17, adjusted P = 3.28E−13, logFC = 3.50) were the two most significant up-regulated genes. Meanwhile, LEP (P = 2.68E−14, adjusted P = 7.06E−10, logFC = − 3.12) and FGF2 (P = 2.84E-14, adjusted P = 7.48E−10, logFC = − 1.87) were the two most significant down-regulated gene in breast cancer tissues.
GO functional enrichment analysis and KEGG pathway enrichment analysis of robust DEGs. To gain insight into the known biological processes and pathways involved in breast cancer, GO functional enrichment analysis and KEGG pathways analysis of 512 robust DEGs were performed. The results showed that those robust DEGs were significantly enriched in 720 GO terms and 9 KEGG pathways, respectively (Supplementary Table S2). Figure 1a-c show the top 20 GO terms, GO terms related to biological process including extracellular structure organization, extracellular matrix organization, ossification, mitotic nuclear division, and regulation of lipid metabolic process (Fig. 1a); Cellular component GO terms were mainly distributed in collagen − containing extracellular matrix, extracellular matrix component, lipid drople and fibrillar collagen trimer (Fig. 1b). The molecular function GO terms consisting of extracellular matrix structural constituent, glycosaminoglycan binding, heparin binding, sulfur compound binding growth factor binding those DGEs were significantly enriched (Fig. 1c). What′s more, KEGG pathway enrichment analysis revealed that PI3K-AKT signaling pathway, PPAR signaling pathway, ECM − receptor interaction, Relaxin signaling pathway, IL − 17 signaling pathway, AMPK signaling pathway were significantly associated with these robust DEGs identified. (Fig. 1d).
EZH2 as a key gene by PPI network analysis. The PPI network of the 512 DEGs, including 493 nodes and 2993 edges was constructed via STRING database (minimum required interaction score: 0.4). By ranking the PPI network nodes using 9 topological analysis methods including both local-and global-based algorithms from cytoHubba plugin of Cytoscape software, we found the EZH2 score ranked in the top 10 by 9 algorithms (Table 2). Furthermore, by performing gene module analysis using MCODE plugin in Cytoscape software, EZH2 gene was also found in Module 1 ( Supplementary Fig. S2), which is the most important module (MCODE score = 31.942) among all modules. In addition, GEPIA database analysis showed that the mRNA expression levels of EZH2 were significantly higher in breast cancer tissues than normal breast tissues ( Supplementary Fig.  S3) and its expression level was significantly associated with the poorer prognosis of patients in breast cancer ( Supplementary Fig. S4). GSEA demonstrated that EZH2 high expression level group was significantly enriched in "Cell cycle" and "DNA replication", which are known to be tumor cell proliferation related pathways (Supplementary Fig. S5).   www.nature.com/scientificreports/ www.nature.com/scientificreports/ module prove our our judgement (Supplementary Table S3-4 showed the GO functional enrichment analysis and KEGG pathways analysis results of genes in blue module and brown module separately). The blue module contained 920 genes and the brown module contained 730 genes. Next, we set the filter standard of hub gene associated with breast cancer: module membership (MM) value > 0.6 and gene significance (GS) value > 0.3, and found that 64 hub genes from the blue module and 38 hub genes from the brown module meet the eligibility criteria (Table 3). By combining with literature searches, four hub genes (CENPL, ISG20L2, MRPL3, and LSM4) were obtained for further analysis. None of these four selected genes had been reported in breast cancer molecular mechanisms studies.

Correlation analysis of the four novel hub genes with clinicopathological variables in breast cancer.
In view of the bc-GenExMiner v.4.6 contains a relatively large number of samples and rich clinical information, we explored the relationships between the expression levels of the novel four hub genes and the clinicopathological variable in this platform. First, we found the expression levels of four novel hub genes (CENPL, ISG20L2, MRPL3, and LSM4) were significantly higher in the subjects aged ≤ 51 years, and high expression of CENPL, ISG20L2, MRPL3, and LSM4 were associated with lymph node metastasis and higher SBR grade (P < 0.05, Fig. 3a−c). Moreover, The TNBC and basal-like BC patients both displayed significantly increased expression of CENPL, ISG20L2, MRPL3, and LSM4 than the non-TNBC and non-basal-like patients, and the expression of CENPL, ISG20L2, MRPL3, and LSM4 were also significantly positively related to Ki67 status ( Supplementary Fig. S6). Thus, we can conclude that CENPL, ISG20L2, MRPL3, and LSM4 are closely related to clinicopathological variables of BC.
Validation of the expression differences of the four novel hub genes. Based on the TCGA_BRCA and match TCGA normal and GTEx data of GEPIA database, the mRNA expression levels of these genes (CENPL, ISG20L2, MRPL3 and LSM4) were also significantly higher in breast cancer tissues than normal tissues ( Supplementary Fig. S7). Immunohistochemistry analysis from The Human Protein Atlas database also showed that these four hub genes were up-regulated in protein expression level in breast lobular carcinoma ( Supplementary Fig. S8). To elucidate the underlying mechanisms of abnormal up-regulation of these four hub genes in breast cancer, we first investigated the association between gene expression and their methylation levels. DiseaseMeth version 2.0 analysis displayed that the mean methylation levels of CENPL, MRPL3 and LSM4 were all significantly reduced in breast cancer compared to normal breast tissues (P < 0.05) (Fig. 4a,c,d). While the mean methylation levels of ISG20L2 significantly increased in breast cancer compared to normal breast tissues (P < 0.05) (Fig. 4b). Additionally, genetic alterations of CENPL, ISG20L2, MRPL3, and LSM4 were further examined in cBioPortal database, showing these four hub genes were altered in 570 (26%) of 2173 breast cancer patients (Fig. 4e). CENPL and ISG20L2 showed the highest alteration levels (20%) with gene amplification as the main alteration type.
Identifying the diagnostic performance of each hub gene in breast cancer. To (Fig. 5b), meaning that the four hub genes all possess good diagnostic performance.

Module
RNA in situ detection. We measured the expression abundance and spatial localization of the five hub genes by RNA in situ detection technology. The results indicate that the expression of these five hub genes were mainly distributed in the cytoplasm and nucleus, and the amount of signal originates from hybridization to each probe varied greatly. Among these five genes, ISG20L2 showed the highest expression level, while CENPL and LSM4 have fewer signal (Fig. 8a-d). Compared to normal mammary epithelial cell (MCF10A), the RCPs of each hub gene in both breast cancer cell lines (MCF7 and MDA-MB-231) showed a significant increase (p < 0.05), but the RCPs of ISG20L2 were reduced and the LSM4 were similar as observed in SKBR3 cell (Fig. 8e-i). Considering the five hub genes, which were identified based on PPI networks and WGCNA, share the same signaling pathways during breast cancer progression, we conducted correlation analysis between the four novel hub genes (CENPL, ISG20L2, LSM4 and MRPL3) and EZH2. As shown in Table 4, the expression levels of each hub gene (CENPL, ISG20L2, MRPL3, and LSM4) was correlated with EZH2 in three different breast cancer cell lines (p < 0.01). Furthermore, we also performed the correlations analysis in GEPIA database to assess the four hub genes correlation with EZH2 in breast cancer and the results remained satisfactory (Supplementary Fig. S10, P < 0.0001). www.nature.com/scientificreports/

Discussion
As a highly prevalent tumor disease worldwide, the complex mechanism involved in the development of breast cancer has not been fully elucidated so far. Thus, identifying potential hub genes involved in breast cancer is not only helpful to elucidating the molecular mechanisms, but also possessing great potentials for the searching of an effective diagnostic biomarkers and prognosis predictors. Previously, hub genes were mainly produced by using a small-scale dataset in most research and showed distressing inconsistent results. In this paper, we applied 10 breast cancer datasets to identify and validate potential novel hub genes so as to guarantee the credibility of the results. First, 512 robust DEGs between breast tumor tissues and normal breast tissues were identified using RRA method. Functional annotation analysis revealed that the robust DEGs significantly enriched in GO terms were associated with proliferation and energy metabolism, such as extracellular matrix organization, extracellular structure organization, mitotic nuclear division, regulation of lipid metabolic process, and glycosaminoglycan binding, which are implicated in the progression of tumor cell [17][18][19] . Meanwhile, the KEGG pathway enrichment analysis showed that 22 genes from those robust DEGs were most associated with PI3K-AKT signaling pathway, which serves as a pivotal intracellular signaling path that plays a crucial role in cell cycle regulation and thus involved in breast cancer development 20 . In addition, we also found that those robust DGEs were significantly enriched in PPAR signaling pathway, EC-receptor interaction, AMPK signaling pathway, and multiple studies have shown these pathways activation are participated in the development and progression of breast cancer and affect the final outcomes [21][22][23] . Based on the PPI networks analysis of robust DEGs, EZH2 was found to be a key gene in the development of breast cancer, and this has been shown in various breast cancer-related studies 24,25 .
Subsequently, a total of 104 hub genes associated with breast cancer were found using the "WGCNA" approach. Among those 104 hub genes, we chose four significantly up-regulated genes, including CENPL,   www.nature.com/scientificreports/ ISG20L2, MRPL3, and LSM4 as hub genes of interest due to the reason that they have rarely been studied in breast cancer and closely related to clinicopathological variables of BC patients. Centromere Protein L (CENPL) is a component of the CENPA-CAD (nucleosome distal) complex, which participated in the assembly process of kinetochore proteins, mitotic progression and chromosome segregation 26 . Interferon stimulated exonuclease gene 20 Like 2 (ISG20L2) encodes a 3'-5' exoribonuclease that involved in the 12S pre-rRNA processing, as a target gene of miR-139-3p, has been reported to take part in the pathogenesis of hepatocellular carcinoma 27,28 . Mitochondrial ribosomal protein L3 (MRPL3), which belongs to the L3P ribosomal protein family, encodes a 39S subunit protein, and plays a regulatory role in the process of Combined Oxidative Phosphorylation 29 . Small nuclear ribonucleoprotein Sm-like4 (LSM4) encodes a member of the LSM family of RNA-binding proteins, has an important role in pre-mRNA splicing by mediating U4/U6 snRNP formation, and this gene has been reported involved in the pathogenesis of pancreatic cancer 30,31 . Interestingly, Joseph S. Baxter et al 32 have also identified that LSM4 is one of 110 target genes at 33 breast cancer risk loci based on Capture Hi-C technology, which supports the conclusions of the present study. In this study, our data not only showed the expression level of CENPL, ISG20L2, LSM4 and MRPL3 were strongly associated with age, lymph node metastasis and higher SBR grade, but also significantly upregulated in triple negative breast cancer patients compared to non-triple negative breast cancer. Furthermore, they are also correlated with the expression of the tumor cell proliferative marker MKI67 (Ki-67), indicating that these four genes may play vital roles in the pathogenesis and progression of breast cancer. Next, ROC analysis revealed that the mRNA expression levels of these four hub genes had excellent diagnostic performance for breast cancer. Prognosis analysis showed that these hub genes were high risk genes, and the higher expression levels were related to poorer prognosis for breast cancer patients. Early diagnosis and accurate evaluation of prognosis play an important role in improving the prognosis of breast cancer patients. thus, these four hub genes also have potentials to serve as promising candidate diagnostic biomarkers and prognosis predictors for breast cancer.
Alternatively, we have also undertaken a preliminary analysis of up-regulation mechanism of these four hub genes refer to DiseaseMeth 2.0 and cBioPortal Database. In DiseaseMeth 2.0 database, the correlation analysis of DNA methylation patterns of hub genes with mRNA expression revealed that CENPL, MRPL3 and LSM4 were hypomethylated in breast cancer samples compared with adjacent normal ones, while ISG20L2 was hypermethylated in breast cancer sample. Generally speaking, the DNA methylation patterns negatively correlates with mRNA expression, which is consistent with the observed up-regulation of CENPL, MRPL3 and LSM4 in breast cancer, while the latest research suggested DNA hypermethylation can lead to mRNA upregulation 33 . Moreover, in cBioPortal database, we found that the abnormal expression of the above hub genes in breast cancer were associated with genetic alterations. In Summary, we think that the gene regulation complexity and up-regulation mechanism of the four hub genes should be further studied.
In addition, we used GSEA to further explore the biological functions of the four hub genes in breast cancer and the results showed that the high-expression groups of CENPL, ISG20L2, MRPL3 and LSM4, were significantly enriched in pathways related to cell proliferation such as "cell cycle" pathway. The results of GSVA were in accordance with the GSEA results. Cell-cycle dysregulation is one of the hallmarks of cancer and several researches have reported that cell cycle disturbance is the most important mechanism for cancer occurrence and progression 34 . Given that the specific functions of those novel hub genes remain unclear, additional research are required to investigate the underlying molecular mechanisms in breast cancer.
Finally, RNA in situ detection technology was applied to detect the five hub genes obtained based on bioinformatics methods, and the results showed that the expression of those five hub genes are different but the spatial localization are similar. While verifying the expression differences of each gene, the correlation analysis showed that the CENPL, ISG20L2, MRPL3 and LSM4 was correlated with EZH2 expression. Since EZH2 are currently more studied in breast cancer, and the four novel hub genes (CENPL, ISG20L2, MRPL3, LSM4) and EZH2 are involved in similar signaling pathways derive from GSEA results, we speculate that each novel hub gene has www.nature.com/scientificreports/ correlation with EZH2 expression, which were verified by the RNA in situ detection results. As an oncogene, growing evidence has identified EZH2 was closely related to breast cancer development and progression through multiple molecular mechanisms 35,36 , but no studies are related to the roles of CENPL, ISG20L2, MRPL3 and LSM4 in breast cancer. Thus, these four novel hub genes as aberrant molecules in the maintenance of breast cancer progression, their exact functional mechanisms deserve further in-depth study.

Collection of breast cancer-related gene expression profile datasets.
Thirteen different datasets comprising eleven datasets from GEO database, one dataset from TCGA database (TCGA_BRCA dataset) and one METABRIC dataset, thousands of breast cancer samples and 280 normal breast tissues samples, were included in our study. The eight series matrix files and corresponding platform annotation information files in RRA analyses were downloaded from Gene Expression Omnibus (GEO) database (http:// www. ncbi. nlm. nih. gov/ geo/), and processed using R package "GEOquery" 37 . The RNA sequencing data normalized by FPKM method, which contains 1066 breast cancer samples and 112 adjacent normal breast tissues samples, were downloaded from The Cancer Genome Atlas (TCGA) data portal (https:// portal. gdc. cancer. gov/) (up to May 01, 2020). At the same time, survival time and vital status of each breast cancer sample in TCGA_BRCA dataset were also extracted and used for subsequent overall survival (OS) analysis. The mRNA expression data and clinicopathological characteristics of 1,904 breast cancer samples in METABRIC dataset were acquired from the cBioPortal website (https:// www. cbiop ortal. org/), of which the mRNA expression levels were determined by Illumina Human v3 microarray and normalized by logarithm 38 . In addition, the raw data of three breast cancer related datasets (GSE21422, GSE42568, and GSE65194) derived from the same microarray platform (GPL570 Affymetrix human genome U133A U133 Plus 2.0 array) were collected separately, then merged and preliminarily cleaned using the "GEOquery" package. The SVA function and Combat function were used to standard and remove the batch effect of three different datasets 39,40 . The merged dataset (GEO_BRCA dataset) was used to validate the diagnostic performance of single hub genes in breast cancer. Furthermore, Breast cancer samples microarray gene expression data from three datasets (GSE20685, GSE20711 and GSE58812) and concomitant follow-up information were also acquired from the Gene Expression Omnibus (GEO) database using "GEOquiry" package and merged with TCGA_BRCA dataset using the "ComBat" function, which served as TCGA_ GEO BRCA dataset (including 1596 breast cancer patients) to explore the prognostic value of novel hub genes.

RRA analysis and identification of robust DEGs.
To discern the DEGs between breast cancer and normal breast tissue in each dataset from GEO database. The "limma" package in R was adopted to normalize the gene expression data and conduct differential gene expression analysis 41 . The differentially expressed genes (DEGs) in each dataset were sorted by their fold change value. subsequently, R package "RobustRankAggreg" 10 was applied to integrate the ranked DEGs of 8 datasets from GEO database so that to find the most important and robust DEGs. Finally, those robust DEGs were determined according to the thresholds: |log2fold change|≥ 1 and P < 0.05.

Pathways and GO function enrichment analyses.
To identify the biological functions and pathways of those robust DEGs, Gene Ontology (GO) Function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using the "clusterProfiler" R package 42 . The GO terms or KEGG pathways with Adjusted P values less than 0.01 indicated statistical significance. Plus, bubble plots were used for visualizing the top 20 enrichment results of GO terms and KEGG pathways.

PPI network construction and analysis of modules.
To identify the key gene of known functions in breast cancer, 512 robust DEGs were mapped to STRING database (STRING, https:// string-db. org/, database version 11.0) to construct a PPI network 43 . Nine topological algorithms in Plug-in CytoHubba 44 , consisting of "MCC", "MNC", "Degree", "BottleNeck", "EcCentricity", "Closeness", "Radiality", "Betweenness" and "Stress" were selected to identify the hub genes in PPI, and the top 10 genes in each topological algorithm were viewed as most stable key gene in PPI analysis. Moreover, the plug-in Molecular Complex Detection (MCODE) 45 in Cytoscape software was also applied to analyze and recognize the modules in the PPI network. All parameters of the above analysis procedure used were set at default values.
WGCNA and potential hub genes identification. To screen potential novel hub genes related to breast cancer, the WGCNA algorithm 46 was used to construct weighted gene co-expression network and identify gene modules that are highly associated with breast cancer. First, gene expression data of the top 3776 up-regulated DEGs obtained by RRA analysis (according to P < 0.05) was extracted from TCGA breast cancer dataset and associated with sample information to construct a sample clustering tree. Second, appropriate soft threshold value (5, scale free R 2 = 0.97) was selected to convert the correlation matrix into adjacency matrix. Subsequently, the resulting adjacency matrix was further converted to topological overlap matrix (TOM) by the TOM similarity algorithm. Referring to the TOM-based dissimilarity calculation formula, these 3776 genes were classified into different gene modules marked by different colors. Third, the minimal module size was set as 50 genes and the height cut-off as 0.25 to merge the highly similar gene modules. Meanwhile, the correlation value between each module's module eigengene (ME) and samples information were calculated using Pearson correlation coefficient. The candidate gene modules were identified based on the degree of correlation between the module's ME values and samples traits. Genes with gene significance (GS) value greater than 0.3 and module membership (MM) value greater than 0.6 in candidate modules were defined as hub genes for breast cancer. These genes may www.nature.com/scientificreports/ have stronger association with the progression and development of breast cancer. Finally, these hub genes were further filtered out based on bioinformatics analyses and literature searches.
Correlation analysis of each hub gene with Clinicopathological Parameters using Breast Cancer Gene-Expression Miner v4.6 (bc-GenExMiner v4.6) online tool. Bc-GenExMiner v4.6 (http:// bcgen ex. ico. unica ncer. fr) 47is a widely used statistical mining online tool for exploring the "correlation", "expression" and "prognostic" analyses of genes of interest in breast cancer by incorporating published annotated breast cancer transcriptomic data (DNA microarrays [n = 11 359] and RNA-seq [n = 4 712]). In this part of the study, the association of each hub genes expression with clinicopathological features in breast cancer was assessed using the "expression" analysis functionalities in bc-GenExMiner v4.6, and the clinicopathological parameters used in this study mainly contained nodal status (N), Age status, Pathological tumor stage, Scarff Bloom & Richardson grade status (SBR), Ki67 status and Basal-like (PAM50) & triple-negative breast cancer status, etc. A p value less than 0.05 was considered statistically significant.
Validation of differential expression of novel hub gene. The Gene Expression Profiling Interactive Analysis (GEPIA, http:// gepia. cancer-pku. cn/) database and The Human Protein Atlas (HPA; http:// www. prote inatl as. org/) database were used to validate the differential expression of each hub gene between breast cancer tissue and normal breast tissue from gene expression and protein levels separately.

Diagnostic performance analyses.
With the aid of R package "pROC" 48 , the receiver operating characteristic (ROC) curves analysis was used to evaluate the diagnostic value of each hub gene using in TCGA_BRCA dataset and GEO_BRCA dataset respectively.

Prognostic value analyses.
To assess the prognostic value of each hub gene, both the samples in META-BRIC dataset and TCGA_GEO_BRCA dataset were divided into high-expression group and low-expression group based on each hub gene's best separation cut-off values. Using built-in "survminer" package and "survival" package in R software, the overall survival (OS) rates were calculated via the Kaplan-Meier method, and the difference in the OS rates between high expression group and low expression group of each hub gene was compared by the log-rank test, P < 0.05 was considered as difference significant. In parallel, hazard ratio (HR) value at 95% confidence interval (95% CI) of each hub gene was also calculated. HR greater than 1 suggested that the gene increase the risk of breast cancer, and HR less than 1 indicated that the gene was a beneficial factor for breast cancer. Moreover, the prognostic value of each novel hub gene was further assessed using the "prognostic" analysis functionalities in bc-GenExMiner v4.6 Platform.
Correlation analysis of methylation level and gene expression of hub genes. The human disease methylation database (DiseaseMeth, version 2.0, http:// bioin fo. hrbmu. edu. cn/ disea se meth/) is a database that integrates massive methylation data from microarray and sequencing results, providing the methylation status annotation information of human diseases 49 . This web database was used to compare the difference of methylation levels of each hub gene between breast cancer and normal breast tissues.
Association analysis of genetic alteration and gene expression of hub genes. The genetic alteration data for each hub gene in the METABRIC dataset samples at the cBioPortal website (http:// www. cbiop ortal. org/) was used to investigate the correlation of genetic alteration and gene expression in breast cancer.
Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) for single hub genes. To find the potential biological functions of single hub genes in breast cancer, R package "clusterprofiler" was chosen to conduct GSEA using METABRAC dataset. Refer to the split-group approach of OS analysis, the samples were divided into "high-expression group" and "low-expression group" based on each hub gene's best expression separation cut-off value. Gene differential expression analysis between each hub gene's "high-expression group" and "low-expression group" was carried out using the "limma" R package. Subsequently, based on the ordered list of all genes according to the logFC value, we performed GSEA using the "clusterProfiler" R package, p. adjust < 0.05 was regarded as statistically significant. Moreover, the GSVA was implemented to verify the differential KEGG pathways of high-expression group and low-expression group via R package "GSVA" 50 . The reference gene set "c2.cp.kegg. v7.0. symbols" were obtained from the Molecular Signature Database (MSigDB, http:// softw are. broad insti tute. org/ gsea/ msigdb/ index. jsp). Cutoff value of differential KEGG pathways was set |logFC|˃0.2, and P < 0.01 was regarded as statistically significant.

RNA in situ detection technology and image analysis. RNA in situ detection technology was used
to determine each hub gene expression at the cellular level. Specific operations were performed referring to literature reported by Ruijie Deng et al 51 and the main steps include: Design of padlock probe complementary to the target RNA, after the padlock probe hybridize to its target and the padlock probe is connected into a ring through specific Splint R DNA ligase, the rolling-circle amplification (RCA) is initiated under the action of primer and Phi29 DNA polymerase, and Finally, fluorescently labeled probes were added to achieve signal detection. For high detection efficiency, three padlock probes were designed for each hub gene in this study (Supplementary Materials: The padlock probe designed for five hub genes). The cell lines used in this study include: MCF10A, MCF7, MD-MB-231 and SKBR3. Image analysis and quantification of signal intensity from each probe was performed in CellProfiler software. A minimum of 1000 cells was counted for each cell line probe www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.