The gene expression datasets GSE20711 and GSE85871 were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo). GSE20711 was comprised of 88 breast cancer samples and 2 normal breast tissue samples, and used the platform GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) (13). GSE85871 was comprised of the gene expression profiles of MCF-7 cells, which were treated with 102 different molecules used in traditional Chinese medicine, and used the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array) (14). The Cancer Genome Atlas (TCGA) Breast Invasive Carcinoma dataset (including high-throughput sequencing (HTSeq) and clinical data of 1,104 breast cancer tissue samples and 113 normal breast samples) was downloaded using the R package ‘TCGAbiolinks (version 2.10.0)’ (15).

The ‘limma (version 3.36.2)’ package was used to load normalized data into R (version 3.3.3; http://www.r-project.org) software and screen the DEGs between breast cancer and non-tumor tissues (16). The genes with fold-change ≥2 and an adjusted P-value (false discovery rate) <0.05 were identified as DEGs (17).

Weighted correlation network analysis (WGCNA) is a commonly used systemic biological data mining method for describing the correlation patterns among genes and identifying the modules of highly correlated genes; it uses average linkage hierarchical clustering coupled with topological overlap dissimilarity based on high-throughput chip or RNA-Seq data (18). The ‘WGCNA (version 1.63)’ package in R was used to construct the co-expression network for the DEGs in the 88 breast cancer samples in GSE20711 (18). β is a soft-thresholding parameter that emphasizes strong correlations between genes and depreciates weak correlations (19). In the present study, β=18 (scale-free R²=0.8) was used to ensure a scale-free network. A cut height of 0.85 and effect size of ≥10 were used to identify the modules. Pearson's correlation matrices were calculated for the modules (20).

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the genes in the modules were performed using the ‘clusterProfiler (version 3.9.1)’ R package based on hypergeometric distribution algorithm (P<0.05), and ‘GOplot (version 1.0.2)’ was used for further analysis (21,22).

To validate the genes in the turquoise modules, the largest modules, the clinical information and RNA sequencing data (HTSeq-FPKM) of breast cancer were obtained from TCGA Project database (https://cancergenome.nih.gov). Kaplan-Meier survival analysis with the log-rank test was conducted to evaluate the association between the genes in the turquoise module and patient survival. P<0.05 was considered to indicate a statistically significant association. Univariate Cox analysis was used to test whether the genes may be used as independent prognostic factors. The data were randomly divided into two groups: A discovery cohort (n=458) and an internal testing cohort (n=457). The genes and clinicopathological characteristics that were significant in the univariate Cox analysis and the Kaplan-Meier survival analysis were used for the multivariate Cox regression analysis. The Akaike information criterion (AIC) value, which was calculated based on different influencing factors by the multivariate Cox regression analysis, was used to remove the confounding factors to obtain the best variable for data fitting, where the minimum AIC value has the best fit (23). Subsequently, a prognostic mRNA and clinical trait signature with min AUC value was constructed, which may be used to calculate a risk score for each individual patient with breast cancer. According to the median of the risk score (median value, 1.915), the patients were stratified into low-risk and high-risk groups. The risk groups from the two cohorts were evaluated using Kaplan-Meier analysis and the log-rank test (24). A nomogram combining the risk score with two other clinical factors (age and sex) was constructed to provide a graphic representation of the prediction model using the R package rms (version 5.1–3.1) (25).

The search tool for the retrieval of interacting genes/proteins (STRING) database (https://string-db.org/) was used to construct the PPI networks, and the results exported from STRING were imported into Cytoscape (version 3.4.0) for visualization (24).

The Human Protein Atlas (http://www.proteinatlas.org/) was used to validate the expression of the three genes in breast cancer tissue (26). The direct links to these images are as follows: RAD54B, https://www.proteinatlas.org/ENSG00000197275-RAD54B/pathology/tissue/breast+cancer#imid_2186145; KIF21A, https://www.proteinatlas.org/ENSG00000139116-KIF21A/pathology/tissue/breast+cancer#imid_17124866; and C8orf76, https://www.proteinatlas.org/ENSG00000189376-C8orf76/pathology/tissue/breast+cancer#imid_6235831.

GSEA was performed using KEGG pathway annotation data from the KEGG database (27). According to the median value of RAD54B expression, patients from the GSE20711 and TCGA datasets were divided into two groups. The clusterProfiler (version 3.9.1) package was used to analyze the data and construct the ridge plot (28).

The MCF-7 cell line was obtained from the American Type Culture Collection. The cells were cultured in complete Dulbecco's modified Eagle's medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) in 5% CO₂ stored at 37°C. The cell cultures were checked for Mycoplasma contamination using a Mycoplasma PCR Detection kit (Sigma-Aldrich; Merck KGaA) every 3 months. The cell line was not listed in the database of commonly misidentified cell lines (https://iclac.org/databases/cross-contaminations).

Japonicone A (>97% purity) was kindly gifted by Professor Weidong Zhang (School of Pharmacy, Second Military Medical University, Shanghai, China) (29–32). Japonicone A was dissolved in dimethyl sulfoxide (DMSO) and the solution was diluted in the cell culture media to ensure that the concentration of DMSO was <0.1%, as described in a previous study (31). The working concentration (10 µM) of Japonicone A was the same as the one applied in the GSE85871 dataset, which was initially used to identify that Japonicone A inhibited the expression of RAD5B in MCF-7 cells (15). In addition, Hu et al (29) used 10 µM Japonicone A in a treatment assay. Preliminary experiments were performed with 10 µM Japonicone A on the proliferation of MCF-7 cells in the present study, and the results suggested that the concentration was effective (data not shown). Therefore, MCF-7 cells were treated with 10 µM Japonicone A for 24 h prior to cell collection and analysis. DMSO without Japonicone A was used as the solvent control. The volume of the solvent control was the same with DMSO-dissolved Japonicone A used in the experimental group.

Total RNA was extracted from MCF-7 cells using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). RNA was reverse-transcribed into cDNA with High-Capacity cDNA Reverse Transcription kit (Applied Biosystems; Thermo Fisher Scientific, Inc.). qPCR was performed using SYBR^® Premix Ex Taq II (Takara Bio, Inc.). Relative mRNA expression levels were calculated by normalizing the relative quantitation cycle value to the control group following standardization to the internal control β-actin (33). The primers used were as follows: Human RAD54B forward, 5′-AAGAACCTGACTGCCTCACG-3′ and reverse, 5′-TCCACCACAGGTAAACCAGC-3′; and human β-actin forward, 5′-CAGGGCGTGATGGTGGGCA-3′ and reverse, 5′-CAAACATCATCTGGGTCATCTTCTC-3′. The thermocycling conditions for the RT-qPCR was as follows: 2 min at 95°C, 40 cycles at 95°C for 10 sec, 59.5°C for 10 sec, 68°C for 15 sec, and 72°C for 10 sec.

The R package ‘survivalROC (version 1.0.3)’ was used for the ROC curve analysis (34). The ROC curve was used to test the sensitivity and specificity of the variables in predicting overall survival, and to assess the predictive ability of the calculated prognostic signature for 5-year patient survival (25).

R package ‘rmda (version 1.6)’ was used for DCA, which estimates the net benefit of a signature by subtracting the false-positives from the true-positives (35,36).

Anti-mouse marker of proliferation Ki-67 (ki-67)-FITC flow cytometry antibody was purchased from Miltenyi Biotec (diluted 1:100 with 1X permeabilization buffer; catalog no. 130-117-691). The intracellular staining of ki-67 was performed using the Foxp3/Transcription Factor Staining Buffer Set (eBioscience; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Briefly, MCF-7 cells were collected and centrifuged at 150 × g for 5 min at room temperature. The cell pellets were mixed with 1 ml fixation-permeabilization buffer (a 1:3 mixture of fixation-permeabilization concentrate and diluent, which were included in the kit) at 4°C for 1 h. The cells were washed twice with 2 ml 1X permeabilization buffer (10X permeabilization buffer diluted to 1X with ddH₂O). The fixed cells were stained with the anti-ki-67-FITC antibody for 30 min at room temperature in the dark. Finally, the washed cells were resuspended in PBS prior to detection using BD FACSVerse (BD Biosciences). The flow cytometric data were analyzed with FlowJo software (version 10.3; FlowJo LLC).

A total of 5×10³ MCF-7 cells were seeded in 96-well plates and treated with the aforementioned amount of control (DMSO) or 10 µM Japonicone A for 24 h. Subsequently, 20 µl MTT (5 mg/ml) was added to the wells and incubated for 4 h. The wells were supplemented with 150 µl DMSO to dissolve the formazan crystals prior to optical density measurement by a microplate reader at 490 nm. There were three duplicate wells for each group, and the cell viability was normalized to the control group in each experiment.

The statistical analyses were performed using R (version 3.3.3). Student's t-test was used for comparisons between two independent groups. P<0.05 was considered to indicate a statistically significant difference. The experiments presented in Fig. 7 were replicated four times and each dot in the histogram represents one independent experiment; lines and error bars represent the mean ± SD.

The GSE20711 samples with complete clinical data were included in the WGCNA analysis (Fig. 1A). DEGs between breast cancer tissues and non-tumor tissues were identified (Table SI). Based on the DEGs, a co-expression network was constructed, and the modules were identified by WGCNA. The results revealed that the most significant DEGs could be grouped into three major modules (turquoise, blue and brown modules) (Fig. 1B). Further analysis demonstrated that the turquoise module exhibited the highest positive correlation with TNM stage and negative correlation with ER status compared with the other modules (Fig. 1C). Thus, this module was identified as the clinically significant module for subsequent analysis.

GO enrichment analysis was performed for the genes in the turquoise module, and the results were categorized into three functional groups: Biological process (BP), molecular function (MF) and cellular component (CC). The genes in the BP group were mainly enriched in ‘mitotic cell cycle process’, ‘cell cycle’, ‘cell cycle process’, ‘mitotic cell cycle’, ‘cell division’ and ‘nuclear division’ (Fig. 2A); the genes in the CC group were enriched in ‘condensed chromosome’, ‘chromosome’, ‘spindle’ and ‘chromosomal region’ (Fig. 2B); the genes in the MF group were mainly enriched in ‘protein binding’, ‘ATP binding’, ‘catalytic activity’ and ‘microtubule binding’ (Fig. 2C). In addition, the KEGG pathway analysis revealed that the DEGs were mainly involved in ‘cell cycle’, ‘oocyte meiosis’, ‘pyrimidine metabolism’, ‘p53 signaling pathway’ and ‘DNA replication’ (Fig. 2D). The PPI network based on the genes from the turquoise module demonstrated that the majority of the genes closely interacted with each other (Fig. 2E). Taken together, the results demonstrated that the DEGs in the clinically significant module mainly participated in the regulation of cancer cell proliferation, and that the majority of them interacted with each other.

The association of the turquoise module with tumor prognosis was analyzed based on TCGA data. Using univariate Cox analysis, several genes (SHCBP1, RAD54B, KIF21A and C8orf76) associated with prognosis were identified (Table I). The results of the Kaplan-Meier survival analysis identified RAD54B, KIF21A and C8orf76 as potential independent prognosis factors, which were negatively associated with overall patient survival time (Fig. 3A-C) and upregulated in late clinical tumor stages (III and IV) compared with early stages (I and II; Fig. S1).

The Human Protein Atlas database, the data of which are open access to all researchers, contains a systems-based analysis of protein expression in 17 types of cancer using data from 8,000 patients, as well as immunohistochemistry images directly demonstrating the indicated protein expression in tumor tissues (37,38). The database was searched to analyze the expression of RAD54B, KIF21A and C8orf76 in normal and breast cancer tissues. With the exception of C8orf76, the protein expression of RAD54B and KIF21A was increased in breast cancer tissues compared with that in normal tissues (Fig. 3D-F).

Since RAD54B, KIF21A and C8orf76 were significantly associated with the overall survival of patients with breast cancer, the patients from TCGA Project database were randomly divided into two groups: The discovery cohort and the internal cohort (Table SII). Combined with the clinicopathological characteristics in the datasets, a prognostic module was developed using forward conditional stepwise regression with multivariate Cox analysis in the discovery cohort. From the AIC values, a prognostic signature containing one gene (RAD54B) and one clinical trait (TNM stage) was identified. Based on this signature, the risk score was calculated for each sample in the discovery cohort using the following formula: Risk score=RAD54B expression × 0.236 + TNM stage (I/II=0 or III/IV=1) × 1.025. The samples were ranked by score and divided into high-risk and low-risk groups based on the median of the risk scores (median, 1.915), which was set as the cut-off point (Fig. 4A). The sample with the median value was assigned to the high-risk group. The patients in the high-risk group exhibited shorter overall survival times compared with the patients in the low-risk group (Fig. 4B and C). In addition, the ROC curve analysis revealed that the AUC was 0.687 based on the risk scores (Fig. 4D).

The risk score model was further evaluated using the internal cohort. Using the risk-score formula and cut-off point derived from the discovery cohort, the patients were divided into high-risk (n=256) and low-risk (n=201) groups. The overall survival of the patients in the internal cohort revealed a similar trend to that in the discovery cohort (Fig. 4E-H). In addition, the nomogram based on the signature and the decision curve analysis demonstrated that this prognostic signature was effective for evaluating the prognosis of patients with breast cancer (Fig. S2).

As RAD54B may serve an important role in breast cancer and may be a potential therapeutic target, the data from TCGA and GSE20711 were divided according to the expression of RAD54B into high and low expression groups (Fig. 5A and B). GSEA analysis demonstrated that the cell cycle function exhibited the strongest association with the expression of RAD54B in the two datasets (Fig. 5C and D). These results suggested that RAD54B may regulate the proliferation of breast cancer cells.

Based on the GSE85871 dataset, the expression of RAD54B in breast cancer MCF-7 cells treated with 102 different molecules used in traditional Chinese medicine was analyzed. The results demonstrated that treatment with a compound known as Japonicone A resulted in downregulated RAD54B expression in these cells (Fig. 6). To confirm this result, an inhibitory assay with Japonicone A on MCF-7 cells was performed. RT-qPCR analysis revealed that the expression of RAD54B was inhibited in MCF-7 cells treated with 10 µM Japonicone A for 24 h (Fig. 7A). In addition, the proportion of Ki-67⁺ MCF-7 cells, which was the indicator for cell proliferation used in the present study, was reduced (Fig. 7B). The viability of MCF-7 cells was determined by MTT assay (39); the results demonstrated that cell viability was decreased following treatment with10 µM Japonicone A (Fig. 7C).