CENPL, ISG20L2, LSM4 and MRPL3 Are Four Novel Hub Genes Involved in Breast Cancer Development and Progression

As a highly prevalent tumor disease worldwide, Further elucidation of the molecular mechanisms of the occurrence, development and prognosis of breast cancer remain an urgent need. Identifying hub genes involved in these pathogenesis and progression can potentially help to unveil its mechanism and provide novel diagnostic and prognostic markers for breast cancer. In this study, we systematically integrated multiple bioinformatic methods, including robust rank aggregation (RRA), functional enrichment analysis, protein-protein interaction (PPI) networks construction and analysis, weighted gene co-expression network analysis (WGCNA), ROC and Kaplan-Meier analyses, DNA methylation analyses and genomic mutation analyses, GSEA and GSVA, based on ten mRNA datasets to identify and investigate novel hub genes involved in breast cancer. In parallel, RNA in situ detection technology was applied to validate those novel hub gene. EZH2 was recognized as a key gene by PPI network analysis. CENPL, ISG20L2, LSM4 and MRPL3 were identied as four novel hub genes through the WGCNA analysis and literate search. Among those ve hub genes, many studies on EZH2 gene in breast cancer have been but no studies are related to the roles of ISG20L2, and in These novel four hub genes were up-regulated in breast cancer tissues and associated with tumor progression. ROC and Kaplan-Meier indicated these four hub genes all showed good diagnostic performance and prognostic values for breast cancer. The preliminary analysis revealed those novel hub genes are four potentially candidate genes for further exploring the molecular mechanism of breast cancer. ISG20L2, four

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 bioinformatics method that can remove substantial inter-study variations and statistical analysis di culties existed in individual studies via integrating the gene expression pro les of different platforms [13]. It has been used in various malignant tumor studies, such as in hepatocellular carcinoma, colorectal cancer, lung cancer and thyroid cancer [14][15][16][17]. The heterogeneity is one of the characteristics of tumor cells, which is re ected genes have different RNA expression patterns at the transcriptome level [18]. 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 uncover the potential function of genes in disease state. RNA in situ detection technology is a new tool 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 quantitative analysis of intracellular RNA at the single cell level [19].
In this study, we rst integrated DEGs from multiple datasets based on RRA algorithms, and then constructed a weighted gene co-expression network using WGCNA algorithms to identify the key gene modules and novel hub genes related to breast cancer. The diagnostic performance and prognostic value of these hub genes were evaluated and their possible molecular mechanisms in breast cancer were explored by bioinformatics methods. At the same time, 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 to analyze the expression differences and correlations to validate those results from the above bioinformatics analysis.

Methods
Collection of breast cancer-related Gene expression pro le datasets Ten different datasets comprising eight datasets from GEO database, one dataset from TCGA database (TCGA_BRCA dataset) and one METABRIC dataset, a total of 3,414 breast cancer tissues and 280 normal breast tissues, were included in our study. The eight series matrix les and corresponding platform annotation information les in each GEO dataset were downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), and further processed using R package "GEOquery" [20]. 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.cbioportal.org/), of which the mRNA expression levels were determined by Illumina Human v3 microarray and normalized by logarithm [21]. 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 the data and remove the batch effect of merged dataset [22,23]. The merged dataset (GEO_BRCA dataset) was used to validate the diagnostic performance of single hub genes in breast cancer.

RRA analysis and identi cation of robust DEGs
In order to discern the DEGs between breast cancer and normal breast tissue in each dataset from GEO database. The"limma" package [24] in R language was adopted to normalize the gene expression data and conduct differential gene expression analysis. The differentially expressed genes (DEGs) in each dataset were sorted by their fold change value. Subsequently, R package "RobustRankAggreg" [13]was chosen to integrate the ranked DEGs of 8 datasets to nd the most important and robust DEGs that were ranked always higher than expected. Finally, according to the thresholds |log2fold change| ≥ 1 and false discovery rate (FDR) < 0.05 determine the robust DEGs.

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 "clusterPro ler" R package [25]. The GO terms or KEGG pathways with Adjusted P values less than 0.01 indicated statistical signi cance. Last, bubble plots were used for visualizing the top 20 enrichment results of GO terms and KEGG pathways.

PPI Network Construction and Analysis of Modules
Protein-protein interaction (PPI) network was constructed for identi cation of hub genes and pivotal gene modules using the search tool for the retrieval of interaction genes (STRING, https://string-db.org/, database version 11.0), which is a publicly accessible online database that provides critical assessment and integration of protein-protein interactions of physical and functional associations [26,27]. In this study, 512 robust DEGs were mapped to STRING database to construct a PPI network. Cytoscape software (version 3.7.1) [28]was used to visualize the PPI network of robust DEGs. Nine topological algorithms in Plug-in CytoHubba, 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 hub genes in PPI analysis. Moreover, the plug-in Molecular Complex Detection (MCODE) [29] in Cytoscape was 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 identi cation
The WGCNA algorithm [30] 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 3,776 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 R2 = 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 3,776 genes were classi ed 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 coe cient. The candidate gene modules were identi ed based on the degree of correlation between the module's ME values and samples. Genes with gene signi cance (GS) value greater than 0.3 and module membership (MM) value greater than 0.6 in candidate modules were de ned as hub genes for breast cancer. These genes may have stronger association with the progression and development of breast cancer. Finally, these hub genes were further ltered based on literature searches.

Correlation analysis of each hub gene with clinical variables
The hub gene expression pro le data and corresponding clinical variables of breast cancer patients were extracted from the METABRIC dataset. Association between each hub gene expression and clinical variables were analyzed using the Spearman correlation coe cient in R software [31]. P-value less than 0.05 was considered statistically signi cant.

Potential hub genes validation
The Gene Expression Pro ling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) database and The Human Protein Atlas (HPA; http://www.proteinatlas.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. With the aid of R package "pROC" [32], 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. To assess the prognostic value of each hub gene in METABRIC dataset and TCGA_BRCA dataset, the samples in both datasets were divided into high-expression group and lowexpression group based on the expression levels of each hub gene's best separation cut-off values. Using built-in "survminer" package and "survival" package in R software [33], 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 signi cant. In addition, hazard ratio (HR) value at 95% con dence 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 bene cial factor for breast cancer.

Correlation analysis of Methylation level and gene expression of hub genes
The human disease methylation database (DiseaseMeth, version 2.0, http://bioinfo.hrbmu.edu.cn/disease meth/) is a database that integrates massive methylation data from microarray and sequencing results, providing the methylation status annotation information of human diseases [34]. 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.cbioportal.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 In order to nd the potential biological functions of single hub genes in breast cancer, R package "clusterpro ler" was chosen to conduct GSEA analysis of KEGG pathway for each hub gene in METABRAC dataset. Refer to the split-group approach of OS analysis, the samples in METABRIC dataset 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 "highexpression 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 "clusterPro ler" R package, p.adjust < 0.05 was regarded as statistically signi cant. Moreover, the GSVA was implemented to discern the differential KEGG pathways of high-expression group and low-expression group via R package "GSVA" [35]. The reference gene set "c2.cp.kegg. v7.0. symbols" were obtained from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp). Cutoff value of differential KEGG pathways was set |logFC| 0.2, and P < 0.01 was regarded as statistically signi cant.

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. Speci c operations was performed referring to literature reported by Ruijie Deng et al [36] 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 speci c Splint R DNA ligase, the rolling-circle ampli cation (RCA) is initiated under the action of primer and Phi29 DNA polymerase, and Finally, uorescently labeled probes were added to achieve signal detection Supplementary Materials: RNA in situ detection protocol). For high detection e ciency, three padlock probes were designed for each hub gene in this study (Supplementary Materials: The padlock probe designed for ve hub genes). The cell lines used in this study include: MCF10A, MCF7, MD-MB-231 and SKBR3. Image analysis and quanti cation of signal intensity from each probe was performed in CellPro ler software. A minimum of 1000 cells was counted for each cell line probe set, and wilcox -test was conducted to compare the rolling circle products (RCPs) of each hub gene between human breast epithelial cell line (MCF10A) and different breast cancer cell lines (MCF7, MD-MB-231 or SKBR3). The expression correlation of these hub genes in breast cancer cell lines were analyzed with Spearman method. P-value less than 0.05 was considered statistically signi cant.

Results
Identi cation 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 was shown in Table 1. We used |log2FC| > 1 and FDR < 0.05 as screening criteria to obtain the robust (DEGs) between breast cancer tissues and normal tissues. A total of 512 robust DEGs were identi ed, containing 202 up-regulated genes and 310 down-regulated genes. 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 signi cant 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 signi cant down-regulated gene in breast cancer tissues. Fig 1 shows the top 20 most signi cant up-regulated and downregulated robust DEGs obtained by RRA methods from these eight different datasets.

GO functional enrichment analysis and KEGG pathway enrichment analysis of robust DEGs
To further understand the biological functions of those 512 robust DEGs, GO functional enrichment analysis and KEGG pathways analysis were performed. These DEGs were signi cantly enriched in GO terms related to biological process including extracellular structure organization, extracellular matrix organization, ossi cation, mitotic nuclear division and regulation of lipid metabolic process (Fig. 2a). Cellular component GO terms were mainly distributed in collagen−containing extracellular matrix, extracellular matrix component, lipid drople and brillar collagen trimer (Fig. 2b). The molecular function GO terms consisting of extracellular matrix structural constituent, glycosaminoglycan binding, heparin binding, sulfur compound binding growth factor binding those DGEs were signi cantly enriched (Fig. 2c). 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 signi cantly associated with these robust DEGs (Fig. 2d) PPI network analysis of robust DEGs Next, the PPI network of the 512 DEGs was constructed via STRING database, including 493 nodes and 2993 edges. To identify the hub genes in this PPI network, we ranked the network nodes using 9 topological analysis methods including both local-and global-based algorithms from cytoHubba plugin of Cytoscape, and nally found the EZH2 score ranked in the top 10 by 9 above-cited consisting of "MCC", "MNC", "Degree", "BottleNeck","EcCentricity", "Closeness", "Radiality", "Betweenness" and "Stress" algorithms (Table 2). Subsequently, gene module analysis was conducted using MCODE plugin in Cytoscape, and the hub gene EZH2 was also found in module 1 (Additional le 1: Figure S1), which is the most important module (MCODE score = 31.942) in all modules. Lastly, GEPIA database analysis showed that the mRNA expression levels of EZH2 were signi cantly higher in breast cancer tissues than normal breast tissues (Additional le 1: Figure S2) and the expression levels of EZH2 was signi cantly associated with the OS of patients in breast cancer (Additional le 1: Figure S3). GSEA demonstrated that high expression level group of EZH2 was signi cantly enriched in "Cell cycle" and "DNA replication", which have been con rmed as tumor cell proliferation related pathways [37](Additional le 1: Figure S4).

WGCNA identi es key gene module and four novel key genes in breast cancer
We performed WGCNA on the TCGA_BRCA dataset that incorporate 3,769 up-regulation DEGs (p value <0.05) derived from the RRA analysis to nd the key modules with breast cancer. After series of quality assessment for gene expression matrix, we set soft threshold as 5 (scale free R 2 = 0.97, slope= −1.92) to construct and validate a scale-free network Fig. 3a-b . By setting minimal module size as 50 genes and cut height as 0.25 to merge similar modules, seven modules were obtained eventually ( Fig. 3c; non-clustering DEGs shown in gray). From the heatmap of module-trait correlations (Fig. 3d), we identi ed that the blue module (cor = 0.44, p = 4e-55) and brown module (cor = 0.46, p = 3e-63) were most correlated with breast cancer (Fig. 3e-f). The blue module contained 920 genes and the brown module contained 730 genes. Next, we set the lter standard of hub gene associated with breast cancer with module membership (MM) value >0.6 and gene signi cance (GS) value >0.3, and 64 hub genes from the blue module and 38 hub genes from the brown module met the eligibility criteria (Table 3). Among these 102 hub genes de ned by GS and MM, four genes CENPL, ISG20L2, MRPL3, and LSM4 were selected for analysis due to they were never reported before in breast cancer.

Correlation analysis of the four novel hub genes with clinicopathological variables in breast cancer
In view of the METABRIC dataset 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 characteristics in this dataset. Clinicopathological variables in METABRIC dataset mainly include age at diagnosis, grade, tumor stage, tumor size and the number of positive lymph nodes. The results of spearman correlation analysis were shown in Fig. 4, showing the expression levels of hub genes (CENPL, ISG20L2, MRPL3, and LSM4) were notably correlated with clinicopathological variables in breast cancer samples (P < 0.05). Higher expression levels of hub genes of all four selected hub genes (CENPL, ISG20L2, MRPL3, and LSM4) were associated with higher grade or latter tumor stage. And higher expression levels of CENPL, ISG20L2, and LSM4 correlated with bigger tumor size. While higher expression of CENPL, ISG20L2 were correlated with more the number of positive lymph nodes.
Validation of the expression differences of the four novel key genes in breast cancercompared to normal tissues 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 signi cantly higher in breast cancer tissues than normal tissues (Additional le 1: Figure S5). Immunohistochemistry analysis from The Human Protein Atlas database also demonstrated the up-regulated of those hub genes expression in protein level in breast lobular carcinoma (Fig. 5). In order to elucidate the underlying mechanisms of abnormal up-regulation of these four hub genes in breast cancer, we 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 signi cantly reduced in breast cancer compared with normal breast tissues (p<0.05) (Fig. 6a, 6c and 6d). While the mean methylation levels of ISG20L2 signi cantly 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 cBiportal database, as shown in Fig. 6e, the four hub genes were altered in 570 (26%) of 2173 breast cancer patients. And CENPL and ISG20L2 altered the most (20%), showing gene ampli cation was the main alteration type.

Identifying the diagnostic performance and prognostic value of eachhub gene in breast cancer
We rst performed ROC analysis to assess the diagnostic performances of the four hub genes for detecting breast cancer in TCGA_BRCA dataset, and their AUC values (CENPL AUC: 0.934, LSM4 AUC: 0.948, MRPL3 AUC: 0.891, ISG20L2 AUC: 0.918) were showed in Fig. 7a, indicating good diagnostic performance. Subsequently, ROC analysis in GEO datasets further validate the diagnostic value of these four hub genes: AUC values showed in Fig. 7b, CENPL AUC: 0.83, LSM4 AUC: 0.913, MRPL3 AUC: 0.841, ISG20L2 AUC: 0.951, meaning that the four hub genes all yielded good prediction performance. Plus, we used TCGA_BRCA dataset and METABRC dataset to perform OS rates analysis to evaluate the prognostic values of these four novel hub genes in breast cancer patients. Although the two datasets composition were inconsistent, the Kaplan-Meier curves still showed that the difference between high expression groups and low expression groups were signi cant (all P < 0.05), and the higher expression levels of these four hub genes were signi cantly associated with the poor OS of breast cancer patients (HR 1, Fig. 8).

GSEA and GSVA exhibit a tight relationship between the four hub genes and tumor cell proliferation
To further elucidate the lurking biological functions of CENPL, ISG20L2, MRPL3 and LSM4 in breast cancer occurrence and development, we conducted GSEA and GSVA analyses based on METABRIC dataset. The results of GSEA were shown in Fig. 9, the genes in high expression groups of CENPL, ISG20L2, MRPL3, and LSM4 were all signi cantly enriched in tumor cell proliferation related pathways such as "cell cycle" and "DNA replication" pathways. Meanwhile, GSVA results substantiated that these cell proliferation-associated genesets were signi cantly up-regulated in the high-expression groups of CENPL, ISG20L2, MRPL3 and LSM4 (Fig. 10).

RNA in situ detection
We measured the expression abundance and spatial localization of the ve hub genes by RNA in situ detection technology. The results show that the signal of ve 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 ve genes, ISG20L2 showed the highest expression level, while CENPL and LSM4 have fewer signal (Fig. 11a-d). Compared to normal mammary epithelial cell (MCF10A), the RCPs of each hub gene in both breast cancer cell lines (MCF7 and MD-MB-231) show a signi cant increase, but the RCPs of ISG20L2 reduced and the RCPs of LSM4 no change was observed in SKBR3 cell (Fig. 11e-i). Considering the ve hub genes, which were identi ed based on PPI networks and WGCNA, share the same signaling pathways during breast cancer progression, we conducted correlation analysis between novel four hub genes (CENPL, ISG20L2, LSM4 and MRPL3) and EZH2. The results of correlations analysis are 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 (Additional le 1: Figure S6a-d, p < 0.0001).

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 potentials for the searching of new drug treatment targets. 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-related datasets to identify and validate potential hub genes related to breast cancer. First of all, 512 robust DEGs between breast tumor tissues and normal breast tissues were identi ed using RRA method. Among them, COL11A1, S100P, LEP and FGF2 are the most signi cant DEGs, which have been reported to be associated with the pathogenesis and development of breast cancer [38][39][40][41]. Functional annotation analysis revealed that the robust DEGs signi cantly enriched in many GO terms 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 were consistent with published data [42][43][44]. Moreover, the KEGG pathway enrichment analysis showed that 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 apoptosis and thus in breast cancer development [45,46]. In addition, we also found that those robust DGEs were signi cantly enriched in PPAR signaling pathway, EC−receptor interaction, AMPK signaling pathway. Multiple studies have shown these pathways are involved in the development and progression of breast cancer and affect the nal outcomes [47][48][49]. Based on the PPI networks analysis of robust DEGs, EZH2 was found to be a hub gene in the development of breast cancer, and this has been shown in various breast cancer-related studies [50][51][52].
Subsequently, we used "WGCNA" R package to construct a weighted gene co-expression network to identify key modules and novel hub genes associated with breast cancer. With the help of Pearson correlation coe cient analysis of gene modules with breast cancer sample traits, we rst identi ed brown module and blue module were two top signi cant modules that were associated with breast cancer. The brown module contained 730 genes and the blue module contained 920 genes, of which 104 hub genes (38 genes in brown module and 64 genes in blue module) in above modules were identi ed according to the lter standard (MM 0.6, GS 0.3). Literature search in PubMed suggested that the majority of the hub genes such as UHRF1, HMMR, AXIN1 and BAX, have been involved in the occurrence, progression and distant metastasis of breast cancer [53][54][55][56]. Among those 104 hub genes, we chose four never reported genes in breast cancer, namely CENPL, ISG20L2, MRPL3, and LSM4, to further explore their correlations with clinicopathological variables as well as their diagnostic and prognostic values.
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 [57]. 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 [58,59]. 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 [60]. 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 [61,62]. Interestingly, Joseph S. Baxter et al. have also identi ed 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 [63]. In this study, we determined that CENPL, ISG20L2, MRPL3 and LSM4 were not only signi cantly up-regulated in breast cancer tissues, but were also positively correlated with the grade, stage, size of tumor and the number of positive lymph nodes, suggesting their important contributions to the pathogenesis and progression of breast cancer. ROC analysis revealed that the mRNA levels of these four hub genes showed excellent diagnostic performance for distinguishing breast cancer from normal breast tissues. Prognosis analysis showed that these hub genes were high risk genes that the higher expression level 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 BC patients, these four hub genes thus also have potentially served as promising candidate diagnostic biomarkers and prognosis predictors for breast cancer.
Alternatively, we have also undertaken a preliminary analysis of up-regulation mechanism of four hub genes refer to DiseaseMeth 2.0 and cBiportal Database. In DiseaseMeth 2.0 database, the corelation analysis of DNA methylation patterns of hub genes with mRNA expression showed that CENPL, MRPL3 and LSM4 were hypomethylated in breast cancer samples compared with adjacent normal ones, whileISG20L2 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 hypermethylate can lead to mRNA upregulation [64]. Moreover, in cBiportal database, we also found that the abnormal expression of the above hub genes in breast cancer were associated with genetic alterations. In Summary, we think that the complexity of gene expression regulation and the mRNA up-regulation mechanism of above 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 signi cantly enriched in pathways related to cell proliferation such as "cell cycle" pathway. The results of GSVA was 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 [65,66]. Given that the speci c function that novel four hub genes have in cell cycle regulation is unknown, additional research is required to investigate those mechanisms in breast cancer in vitro and in vivo.
Finally, RNA in situ detection technology was applied to detect the ve hub genes screened based on bioinformatics methods, and the results showed that the expression of those ve hub genes have different expression abundance and similar spatial localization at the transcriptional level. 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 correlation with EZH2 expression, and RNA in situ detection veri ed this inference. As an oncogene, growing evidence has identi ed EZH2 was closely related to breast cancer development and progression through multiple molecular mechanisms [67,68]. but never studies are related to the roles of CENPL, ISG20L2, MRPL3 and LSM4 in breast cancer. Thus, these four potential hub genes serve as aberrant molecules in the maintenance of breast cancer progression, their exact functional mechanisms deserve further in-depth study.

Availability of data and materials
Gene expression microarray datasets were downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). RNA-seq data and corresponding clinical of TCGA_BRCA was acquired from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/) (up to May 01, 2020). The mRNA expression data and clinicopathological characteristics of METABRIC dataset was obtained from the cBioPortal website (https://www.cbioportal.org/).

Ethics approval and consent to participate
None.

Consent for publication
The work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part.