Super-Enhancer-Associated Hub Genes In Chronic Myeloid Leukemia Identified Using Weighted Gene Co-Expression Network Analysis

Introduction

Chronic myeloid leukemia (CML), a homogeneous genetic disease and clonal myeloproliferative disorder of pluripotent hematopoietic stem cells, is mainly triggered by a Philadelphia (Ph) chromosome encoding the BCR-ABL oncogenic fusion protein with constitutive and aberrant tyrosine kinase activity.^1–3 While the causes and mechanisms underlying progression of CML in most cases remain unknown, it is suggested that genetic variations and aberrant expression of key genes play a role in the pathogenesis of the disease.⁴

Recent studies have identified hundreds of hub genes in critical signaling pathways involved in oncogenesis and that may be effectively utilized as therapeutic targets of CML.⁵ BCR-ABL, a fusion tyrosine kinase, was initially determined as a therapeutic target for patients with CML.⁶ Several other genes and pathways have been reported as potential prognostic markers and drug-sensitive indicators for CML, including ephrin type-B receptor 4 (EPHB4),⁷ Janus kinase 2 (JAK2),⁸ epidermal growth factor receptor (EGFR),⁹ β-catenin (CTNNB1),¹⁰ vascular endothelial growth factor A (VEGFA),¹¹ KIT proto-oncogene receptor tyrosine kinase (c-Kit), and tumor protein p53 (TP53).¹² Due to the small sample sizes of the studies performed to date, the molecular mechanisms associated with clinical traits, such as development and progress of CML, are yet to be fully elucidated.

Weighted gene co-expression network analysis (WGCNA) is a widely used systems genetic data analysis strategy based on pairwise correlations between variables.^13–15 WGCNA is used to define modules, intramodular hubs, and network nodes with regard to module membership to determine the relationships between co-expression modules and compare the topology of different networks, thereby defining the significant eigengenes related to clinical traits.¹⁴ WGCNA has been extensively applied for analyzing genomics and metabolomics data, including microarray data,^16,17 single-cell RNA-Seq data,¹⁸ DNA methylation, data¹⁹ and non-coding RNA data,^20,21 peptide counts,^22,23 and microbiota 16sRNA data,²⁴ and shown to be suitable for investigating integrated co-expression networks obtained with large-scale samples.

Super-enhancers (SEs) are a special group of enhancers, which have since been identified in many cell types.^25–27 Compared with typical enhancers (TEs), SEs recruit an exceptionally large number of transcription factors/co-factors and consequently induce transcription of multiple target oncogenes. H3K27ac and H3K4me1 are both indicators for identifying active enhancers. Compared to H3K4me1, mediators (CDK7, BRD4) and other biomarkers, H3K27ac is a common and efficient biomarker for identifying super-enhancers.^25–27 Aberrant expression of oncogenes triggered by SEs is a key factor underlying tumorigenesis, so the screening and identification of hub oncogenes-driven by SEs attract many researchers.

In the present study, we developed a WGCNA-based method to effectively identify the SE-associated hub genes involved in clinical characteristics of CML. The CML microarray datasets and corresponding clinical traits were applied to perform WGCNA and determine key eigengenes. SE-associated genes were obtained using the dbSuper database based on H3K27ac ChIP-seq.²⁸ Node genes were selected on the basis of gene overlap in the yellow module and SE-associated genes. SE-associated hub genes were further predicted according to the signal peaks of H3K27ac, H3K4me1 and BRD4 ChIP-seq. The BRD4-specific inhibitor, JQ1, was additionally applied to identify SE-associated genes in the CML cell line K562. Elucidation of the mechanisms associated with upregulation of these hub genes by SEs should facilitate our understanding of the tumorigenesis process and aid in improving therapeutic options.

Materials And Methods

Data Collection

CML datasets of GSE77191 were available from the NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). GSE77191 is composed of 4 normal and 36 CML samples, including 26 chronic phase (CP), 7 accelerated phase (AP), and 3 blast phase (BP) samples. Annotation of the probe and clinical information was downloaded from the [HTA-2_0] Affymetrix Human Transcriptome Array 2.0 platform [transcript (gene) version]. Raw data were processed to generate an expression matrix and all probes matched to their gene symbol using the R package “limma” and annotation information on the GPL17586 array platform.

SE-associated genes in K562 CML cells were downloaded from the dbSuper database (http://asntech.org/dbsuper/).²⁸ Master transcriptional factors and mediators (SOX2, OCT4, Nanog and BRD4) as biomarkers were reported to be used to identify SEs.^25–27 H3K27ac and H3K4me1 as active enhancer markers are frequently used as biomarkers for the identification of SEs.^25–27 In the present study, according to the algorithm of rank ordering of super-enhancers (ROSE), H3K27ac ChIP signal peaks (GSM733656) within 12,500 bp were stitched into one peak and considered as a super-enhancer, while peaks entirely contained within a window of ± 2,000 bp around transcription start site (TSS) were excluded from stitched regions.²⁸ SE-associated genes were in the closest proximity to their corresponding SEs. We also use BRD4 (GSM2700494) and H3K4me1 (GSM788085) ChIP-seq to verify super-enhancer.

Construction Of WGCNA

The R package of “WGCNA” provides functions for carrying out all aspects of weighted network analysis (module construction, hub gene selection, module preservation statistics, differential network analysis, and network statistics). Here, WGCNA was used to construct a co-expression network for whole genes in samples from GSE77191 with the corresponding clinical information. Firstly, we employed Pearson’s correlation analysis to identify the outlying microarray samples. A matrix of similarity using Pearson’s correlation analysis of all gene pairs was constructed, which excluded outlying samples. The power value of the appropriate soft-thresholding was estimated for network construction by calculating the scale-free topology fit index (0.9) for several powers through the pickSoftThreshold function of WGCNA. The best-fit power value was selected to raise the matrix of similarity to achieve a scale-free co-expression network. Cluster analysis was subsequently performed with flashClust. Moreover, adjacency matrix was converted into a topological overlap matrix (TOM). According to the TOM-based dissimilarity measure with a minimum size (gene group) of 30 for the gene dendrogram, average linkage hierarchical clustering was conducted and genes with similar expression profiles classified into the same modules using the dynamic tree cut algorithm (deep split = 2, cut height = 0.25; default values were used for the other parameters).²⁹ Module eigengenes (MEs) with dissimilarities were merged into the same modules under the cut line of 0.25.

Identification Of Clinically Significant Modules

Two approaches were applied to establish the relationships between modules and clinical phenotypes of CML via WGCNA. Firstly, the Pearson’s correlation coefficient was calculated to analyze the correlation between MEs and clinical traits. In addition, gene significance (GS) was calculated as the log₁₀ transformation of the p-value of each gene (GS = LgP) in the linear regression slope between gene expression and clinical information. Module significance (MS) representing the average GS value for all the genes in a selected module was measured to incorporate clinical information into the co-expression network. The module with the highest MS value was considered the key module related to clinical traits.

Functional Annotations Of Clinically Significant Modules And SE-Associated Genes

Metascape (http://metascape.org/gp/index.html) is a web-based tool for gene annotation and functional enrichment analysis of pathways.³⁰ Eigengenes within the clinically significant modules and SE-associated genes were annotated and analyzed using Metascape with default parameters, and p-values < 0.05 were regarded as significant.

Identification Of SE-Associated Hub Genes In Clinically Significant Modules

Venn diagrams were obtained with the online tool Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html) to identify the genes intersecting between eigengenes from clinically significant module and SE-associated genes. A protein–protein interaction (PPI) network was constructed using overlapping genes via Networkanalyst (https://www.networkanalyst.ca/). The network diagram of all PPI node genes was visualized with Cytoscape.³¹ SE-associated hub genes were further selected based on the rank of node degree. Differential expression of node genes in CML patients was examined via Heatmap and Dotplot. Due to broad usage of the predicted markers of SE, signal peaks of H3K27ac, H3K4me1 and BRD4 ChIP-seq were downloaded from the dbSuper database and visualized using Integrative Genomics Viewer (IGV).³²

Identification Of SE-Associated Hub Genes In CML Samples And K562 Cells

Thirty healthy and thirty CML blood samples were collected with the approval of Medical Ethics Committee of the Xiangya Hospital of Centre South University (201502032). The blood donors have written informed consent in accordance with the declaration of Helsinki. K562 cells were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) where they were characterized by mycoplasma negative, DNA-Fingerprinting, isozyme detection and cell vitality detection. K562 cells cultured in RPMI-1640 were treated with 1 or 10 μM JQ1 (CSNpharm, Shanghai, China, Cat # CSN13058) for 24 hrs. Total RNA was extracted using TRIzol (Thermol Fisher Scientific, Shanghai, China, Cat # 15596026) and expression levels of SE-associated hub genes identified using SYBR Green quantitative real-time PCR (TOYOBO, shanghai, China, Cat # QPK-201). Primer sequences were downloaded from PrimerBank.³³ The antibodies, including BUB1 (Cat # 14H5), KIF2C (Cat # 2488C3a), ORC1 (Cat # F10), RRM2 (Cat # ab57653) and CENPO (Cat # ab173489) were purchased from Santa Cruz Biotechnology (Shanghai, China) and Abcam (Shanghai, China). Protein levels of SE-associated genes were analyzed via Western blotting according to the previous methods in our publications.³⁴

Statistical Analysis

Statistical analysis of differential gene expression was conducted by a one-way ANOVA and post-hoc analysis by Student–Newman–Keuls test using SPSS version 13.0 (SPSS, Inc., Chicago, IL, USA). Biological differences were considered significant at p-values < 0.05.

Results

Construction Of The Weighted Gene Co-Expression Network

To build the weighted gene co-expression network, raw data of GSE77191 were downloaded from the GEO database. In total, 40 samples were enrolled for study, including 26 CML CP, 7 CML AP, 3 CML BP, and 4 normal samples. A total of 3538 probes were annotated and gene expression preprocessed identically using the “limma” R package for background correction and normalization. Expression for each gene was calculated and ranked from large to small, and the top 3000 genes selected for WGCNA. Cluster analysis was further performed using the flashClust function of the WGCNA package (Figure 1A). GSM2045756 as an outlier sample was removed in the following analysis.

Figure 1 Clustering of samples and determination of soft-thresholding power. (A) The clustering heat map was based on expression data of GSE77191, including 26 CML CP, 7 CML AP, 3 CML BP, and 4 normal samples. The top 10,000 genes with the highest SD values were used for WGCNA analysis. Color intensity was proportional to sample outliers, disease status (CML and normal control), phase (normal control, CP, AP, BP), sex, and age. (B) Analysis of the scale-free fit index (y-axis) for various soft-thresholding powers (β value of x-axis). (C) Analysis of the mean connectivity (y-axis) for various soft-thresholding powers (β value of x-axis).

To construct the gene co-expression network, the power value of appropriate soft-thresholding was obtained through prediction of the scale-free topology fit index (scale-free R²=0.90) for several powers. We selected β = 10 as the soft-thresholding power to ensure a scale-free network as the same analysis of the mean connectivity for various soft-thresholding powers (Figure 1B and C). A power of β = 10 was applied to produce a hierarchical clustering tree of module eigengenes (Figure 2A). We set the cut line at 0.25 to merge similar modules (Figure 2A), consequently generating 6 merged modules (Figure 2B).

Figure 2 Construction of co-expression modules using WGCNA. (A) Cluster dendrogram of module eigengenes. (B) Dendrograms of differentially expressed genes clustered based on dissimilarity measures (1-TOM). Each short vertical line represents one gene, and the area of each color signifies the number size of each cluster of co-expression modules.

Identification Of Clinically Significant Modules

We further analyzed the interactions of eigengenes in six modules and mapped the network heat map based on randomly 500 selected genes (Figure 3A). A module–trait relationship was calculated according to the Pearson correlation coefficient and displayed in the heat map. The yellow module displayed a higher correlation with CML status (r = 0.44, p = 0.005) and phase (r = 0.46, p = 0.003), compared with other modules (Figure 3B), with the highest module gene significance (GS) in relation to CML status (p = 3.5e-156) and phase (p = 8e-86) (Figure 3C, D). A scatter plot in Figure 3E illustrates the significant relationship between yellow module membership and GS. The heat map and column chart of gene expression in the yellow module further confirmed the significance of the relationship (Figure 3F). Accordingly, the yellow module was identified as the most significant clinically related unit and selected for subsequent analyses.

Figure 3 Identification of key modules associated with the clinical traits of CML. (A) Topological overlap matrix plot. The different colors on the horizontal and vertical axes represent different modules. The brightness of the yellow in the middle represents the degree of connectivity of different modules. We observed no significant differences in interactions among the different modules, indicating a high-scale independence degree. (B) Heat map of the correlation between modules and clinical traits of CML. Numbers in the table signify the correlations of the corresponding module eigengenes and clinical traits, with the p-values specified below the correlations in parentheses. (C and D) Distribution of the average gene significance and errors in the modules associated with status and phase of CML. (E) Scatter plot of module eigengenes related to status and phase in the yellow module. (F) Heat map and column chart of gene expression in the yellow module. Intensity and direction of the correlations are indicated in the heat map (red, positively correlated; green, negatively correlated).

Functional Annotation Of Clinically Significant Modules And SE-Associated Genes

SEs in K562 were downloaded from dbSuper and 806 SE-associated genes estimated according to the “ROSE” algorithm based on the H3K27ac ChIP-seq signal (Figure 4A). To identify the functions of eigengenes within the yellow module and SE-associated genes in K562, 254 eigengenes (Supplementary Table 1) and 806 SEs-associated genes (Supplementary Table 2) were annotated for pathway analysis via Metascape. Heat map was applied to determine the top enrichment clusters for SE-associated genes and eigengenes (Figure 4B) (Supplementary Table 3). The rows represent SE-associated genes and eigengenes, with a discrete color scale used to indicate statistical significance (-Log P). Both SE-associated genes and eigengenes in the yellow module were significantly enriched in eight pathways (red box). Metascape enrichment network visualization disclosed intra-cluster and inter-cluster similarities of enriched terms between SE-associated genes and eigengenes (Figure 4C and D). One term from each cluster was selected for description, shown as a color label. Common pathways are additionally depicted in the figure, whereby nodes are represented by pie charts, indicating their associations with each input gene.

Figure 4 Pathways and functional annotations for SE-associated genes and eigengenes in the yellow module. (A) Prediction of SEs and SE-associated genes. SEs were predicted using the ROSE algorithm, and SE-associated genes were located closest to their respective SEs. A dot plot of the rank of enhancers based on the H3K27ac ChIP-seq signal was generated. In total, 806 SEs and SE-associated genes were predicted (red dots). Blue dots represent typical enhancers (TE). (B) Heat map showing top enrichment clusters for SE-associated genes and eigengenes, one row per cluster, with a discrete color scale representing statistical significance (-Log P). Gray color indicates lack of significance. (C) Metascape enrichment network visualization showing intra-cluster and inter-cluster similarities of enriched terms for SE-associated genes and eigengenes in the yellow module. Each term is represented by a circle node, where its size is proportional to the number of input genes falling into that term, and its color represents cluster identity. Terms with a similarity score > 0.3 are linked by an edge (the thickness of the edge represents the similarity score). One term from each cluster was selected to have its description shown as a label. (D) Enrichment network visualization of results from SE-associated genes and eigengenes in the yellow module, whereby nodes are represented by pie charts indicating association with each input study.

Identification Of SE-Associated Hub Genes Related To CML Status And Phase

To identify the common genes between SE-associated genes and genes from yellow module, 25 node genes of the intersections were identified in the Venn diagram (Figure 5A). 25 genes were further uploaded into the STRING database for PPI analysis (Figure 5B). We visualized the 25 node genes and other connecting genes in the PPI network using cytoscape and identified 12 hub genes by sorting the node degree of candidate genes (Figure 5B and C). UBC was excluded in the following study. UBC is a selected hub gene in the PPI network, but is not SE-associated gene. The hierarchical clustering heat map further disclosed expression levels of 11 node genes in different disease status groups, including CP, AP, and BP, compared to normal control (Figure 5D). Five SE-associated hub genes showing significantly different expression levels were selected to analyze ChIP-seq tracks of SE. As shown in Figure 6, signal peaks of H3K27ac in blue were used to predict SE regions. H3K4me1 and BRD4 ChIP-seq were also used to display in pink and green, respectively. In Figure 6, predicted SEs and several constituent active enhancers are shown in red. The five SE-associated hub genes identified (BUB1, CENPO, KIF2C, ORC1, and RRM2) were located close to the predicted SE. Strong peak signal of H3K4me1 and BRD4 is also observed at the predicted active enhancer. Thus, these biomarkers are favorable biomarkers for the identification of super-enhancer. In the boxplot diagram (Figure 7A), expression levels of the five genes from GSE77191 in the normal control group were markedly lower than those of CML (CP, AP, and BP) groups (p < 0.05). These hub genes also showed higher mRNA expression in CML patients, compared to healthy controls (Figure 7B). Expression of the genes gradually decreased after treatment with JQ1 (1 or 10 μM) at both the mRNA and protein levels (Figure 7C and D). Our results clearly demonstrate that expressions of SE-associated hub genes in both CML patients and K562 cells are probably upregulated by SEs.

Figure 5 Screening of SE-associated hub genes. (A) Venn plot of SE-associated genes and yellow module eigengenes. Genes intersecting within the Venn plot were considered candidate hub genes. (B) Plot of ranking node genes based on node degree. Genes with a node degree of >3 were considered hub genes. (C) PPI network of SE-associated genes constructed using NetworkAnalyst. The square and oval boxes represent node and hub genes, respectively. The continuous color map depicts the node degree. (D) Heat map of SE-associated hub gene expression in CML patients and control groups.

Figure 6 Signal tracks for H3K27ac (blue), H3K4me1 (pink) and BRD4 (green) ChIP-seq profiles of SE-associated hub genes in K562 cells visualized using IGV. The stitched regions of active enhancers and SEs are shown in red.

Figure 7 Identification of SE-associated hub genes in CML patients and K562 cells. (A) Boxplots showing mRNA expression of SE-associated hub genes in CML patient groups (CP, AP, and BP) (GSE), compared to the control group. Differences were considered significant at p-values < 0.05 (*). Compared to different phase groups, p-values < 0.05 were considered significant (#). (B) Boxplots and dotplots showing mRNA expression of SE-associated hub genes in CML patients from Xiangya Hospital. Differences were considered significant at p-values < 0.05 (*), compared to the control group. (C and D) mRNA and protein levels of SE-associated genes in K562 cells after treatment with JQ1 (1 and 10 μM). Differences were considered significant at p-values < 0.05 (*), compared to the control group.

Discussion

CML is characterized by accumulation of myeloid cells in blood and caused by uncontrolled growth of myeloid cells in the bone marrow.⁴ Despite improvements in CML treatments during the last decade, the ability to treat advanced cases remains limited due to the lack of precise molecular targets. Thus, identification of novel biomarkers potentially involved in specific prognosis and progression of CML is essential to improve therapeutic options.

In this study, we identified the key modules and hub genes in tumorigenesis and phases of CML using WGCNA. While genome-wide gene expression results from CML datasets are available, analyses that can effectively link expression patterns to clinical traits remain a major challenge.³⁵ WGCNA methods avoid this limitation by focusing on a group of genes rather than individual genes through screening the modules and eigengenes associated with clinical traits. The top 3000 genes based on SD ranking were used to construct a co-expression network and six significant modules were subsequently identified. Based on the analysis of module–trait relationships, the yellow module was most significantly associated not only with CML disease status but also with CML phase. Within the yellow module, 254 genes with high connectivity were screened and enrichment analysis of signaling pathways and GO function conducted. Genes in the yellow module and SE-associated genes were mostly enriched in the pathways involved in cell cycle phase transition, megakaryocyte differentiation, Ub-specific processing proteases, hemostasis, regulation of chromosome organization, signaling by interleukins, and centromere complex assembly.

Super-enhancers are an emerging sub-class of regulatory regions that recruit master regulators and co-activators, and control cell identity and disease-related genes.^36,37 Recent findings indicate that SE-associated specific oncogenes may be required for tumor cell proliferation but are not observed in normal cells.³⁷ The utility of SEs in the therapeutic targeting strategies has been minimally explored to date.³⁸ Given the pivotal roles of SE-associated genes in determining cancer cell identity and supporting cancer cell growth, we speculated that the SE-associated genes screened from eigengenes in the clinically related module are also critical for CML tumorigenesis and development. In our study, 25 SE-associated genes from K562 cells were selected and further subjected to PPI analysis.

Based on the ranking of node genes in the PPI network and differences in expression levels between the normal control and CML groups, five SE-associated genes were identified as hub genes. In particular, KIF2C and RRM2 displayed significant differences between the CP and AP phases of CML while the other genes were associated with status but not phase. Following the treatment of K562 with a specific BRD4 inhibitor, JQ1, the levels of five SE-associated genes were reduced. These findings imply that five hub genes significantly related to CML status and phase are upregulated by SEs and suppressed by JQ1. Although KIF2C and BUB1 connected each other in Figure 5, there is no report about this connection and binding partners. In future, the interaction of KIF2C and BUB1 is probably interfered by drugs as a therapeutic target. Other interactions between the five SE-associated genes were not found. Based on the signal peaks of H3K27ac ChIP-seq shown in blue, five SE-associated hub genes (BUB1, CENPO, KIF2C, ORC1, and RRM2) were identified close to their respective predicted SEs. The high signal peaks of H3K4me1 and BRD4 ChIP-seq also showed near SE-associated hub genes. Our results support the enhancement of expression of the hub genes identified in CML patients. While several critical SEs have been previously reported, it remains unclear whether the regulation of these genes is relevant to the development and status of CML.^26,39 The SE-associated gene CCND3 was reported in an earlier study, but was not significantly changed in CML patients in GSE77191 dataset (p > 0.05).^26,39 Thus, these SE-associated genes probably play pivotal roles in K562 cell line but not in CML patients. Here, we focused solely on the most important module and genes displaying the most significant fold changes. Thus, several SE-associated genes in K562 cells previously reported as important were excluded from our study due to the indistinct expression changes.

BUB1 encodes a serine/threonine–protein kinase that functions, in part, by phosphorylating members of the mitotic checkpoint complex and activating the spindle checkpoint.⁴⁰ BUB1 additionally plays a role in abrogating activation of the anaphase-promoting complex/cyclosome.⁴⁰ Germline mutations or aberrant expression in BUB1 have been associated with aneuploidy and several forms of cancer.^41–43 The BUB1 inhibitor, BAY 1816032, is reported to sensitize tumor cells to taxanes, ATR Serine/Threonine Kinase (ATR), and Poly(ADP-Ribose) Polymerase (PARP) inhibitors in human triple-negative breast xenograft models.⁴⁴ These findings highlight a role of BUB1 in inducing a state of tumorigenesis and suggest promising new therapeutic avenues for treating highly invasive cancer.

Kinesin Family Member 2C (KIF2C), also known as mitotic centromere-associated kinesin (MCAK), functions as a microtubule-dependent molecular motor involved in related cell cycle, mitotic, and reelin pathways.⁴⁵ The function of KIF2C is correlated to its conformational changes and depolymerization activity.⁴⁵ A model of its regulation by multiple mitotic kinases has been proposed and its potential involvement in oncogenesis and drug resistance highlighted.⁴⁵ KIF2C serves as a novel prognostic marker in human gliomas, and suprarenal epithelioma.^46,47 Expression of this protein is also correlated to poor prognosis of operable esophageal squamouscell carcinoma in male patients.⁴⁸ Thus, repression of KIF2C may present a novel therapeutic strategy for combating malignancy in some tumor entities.

Ribonucleotide reductase regulatory subunit M2 (RRM2) encodes one of two non-identical subunits for ribonucleotide reductase involved in interrelated cell cycle, mitotic, and pyrimidine metabolism pathways. Integration of transcriptomic data revealed significantly elevated RRM2 in glioblastoma,⁴⁹ breast cancer,^50–53 prostate cancer,⁵⁴ neuroblastoma,⁵⁵ colorectal cancer,^17,56 pancreatic cancer,⁵⁷ hepatocellular carcinoma,^58,59 and adrenocortical cancer.⁶⁰ Knockdown of RRM2 with siRNA effectively impeded pancreatic tumor growth either alone or synergistically with doxorubicin.⁶¹ Delivery of RRM2 siRNA to vascular smooth muscle cells through liposome–polycation–DNA complex conjugated with cell-penetrating peptides could also markedly inhibit RRM2 mRNA and protein expression, resulting in significant suppression of cellular proliferation and migration.⁶² Trans−4, 4ʹ-dihydroxystilbene (DHS) was shown to suppress DNA replication by inhibiting RRM2, thereby decreasing the growth of pancreatic, ovarian, and colorectal tumor.⁶³ In summary, RRM2 as a hub gene and related pathways may serve as biomarkers or therapeutic targets for cancer.

ORC1, a highly conserved six-subunit protein complex, binds specifically to origins of replication and serves as a platform for the assembly of additional initiation factors, such as Cdc6 and Mcm proteins.⁶⁴ ORC1 is reported to be involved in cell cycle, mitosis, and E2F-mediated regulation of DNA replication.^65,66 Recent findings revealed that knocking down of ORC1 memorably suppressed cell proliferation, blocked cell cycle, decreased the expression of Bcl-2 while increased the apoptosis rate and the expression of c-caspase3 and cleaved PARP in cervical cancer cell lines HeLa and C33A.⁶⁷ CENPO encodes a component of the interphase centromere complex, which localizes to the centromere throughout the cell cycle and is required for bipolar spindle assembly, chromosome segregation, and checkpoint signaling during mitosis.⁶⁸ Although few reports have linked ORC1 and CENPO to cancer treatment to date, the findings that these genes are enriched in the signaling pathways of cell cycle and mitosis support their utility as novel therapeutic biotargets for CML.

Conclusion

This is the first study to apply WGCNA-based methods to identify SE-associated hub genes expressed in association with CML clinical traits. Further investigation of these genes should not only enhance our understanding of CML tumorigenesis and development but also provide potential therapeutic targets. Our study presents proof of principle applicable to other diseases and clinical traits, which can be used to redefine molecular mechanisms and disease phenotypes for optimization of clinical practice.