This article is mentioned in:

Introduction

Acute lymphoblastic leukemia (ALL) is the most frequent hematologic malignancy in pediatric patients, and ALL develops due to a combination of genetic susceptibility and environmental exposure in fetal life and infancy (1). In the United States, the incidence rate (age-adjusted) was reported to be 4.9 per 100,000 children from 2007 to 2011 (2). Despite advancements in the management of pediatric ALL (B-cell and T-cell ALL), that have occurred, the incidence of the disease has exhibited an increasing trend (3). Post-transplant relapse remains the primary causative factor for treatment failure in 30% of patients with ALL (4). The disease relapse may be due to the chemotherapy resistance of leukemic cells following the initial treatment in the majority of these patients (5). Therefore, there is a requirement for novel effective therapeutic methods for the treatment of ALL.

Gene expression microarrays have emerged as an effective technique for identifying differentially expressed genes (DEGs) in a comprehensive approach. Previously, by using higher density oligonucleotide arrays, ~60% of the newly identified subtypes of discriminating genes were recognized as novel markers (6). Additionally, 54 genes were identified with distinguished resistance from sensitive ALL samples by comparing gene expression signatures of pediatric ALL tissue samples with or without high minimal residual disease (MRD) load (7). Additionally, a number of genes (protein products) were identified with increased expression levels (1.5–5.8 fold) in leukemic cells compared with control samples, including cluster of differentiation (CD)58 and human DNAJ protein (HDJ)-2 (8). These previous studies demonstrated that a large cohort of DEGs may be involved in the progression of ALL. Small molecule drugs have been identified as potential inhibitors of ALL, including imatinib, an Abelson tyrosine-protein kinase 1 inhibitor (9). Additionally, borrelidin was established as a potent agent inhibiting the cell proliferation in ALL cell lines via the induction of apoptosis and the regulation of G1 arrest (10). A previous study which established the Gene Expression Omnibus (GEO) DataSet gene expression array GSE42221 identified that the epigenetic and transcriptional regulation of zinc finger protein (ZNF)423 was associated with B-cell differentiation in ALL (11,12). However, the study focused on the causative role of the gene ZNF423 in ALL. Associated genes and their potential interactions, as well as small molecule drugs, were not considered in Harder et al (11).

Therefore, the current study re-analyzed the published expression profiles (11) and consequently identified DEGs in pediatric B-cell ALL. Additionally, functional enrichment analysis of DEGs was performed to investigate the dysregulated biological processes using the altered DEGs in the progression of the disease. In addition, the protein-protein interactions (PPIs) of the DEGs were screened and small molecular drugs associated with DEGs were analyzed using a number of bioinformatic methods. All these analyses were aimed to provide effective and novel therapeutic methods for the management of pediatric ALL.

Materials and methods

Data source

Gene expression data with the accession number of GSE42221, which was deposited by Harder et al (11) in the public National Center for Biotechnology Information GEO database, were downloaded. The dataset consisted of 7 primary human B-precursor samples (ALL samples) and 4 normal B-cell progenitor lymphoblast samples (control samples). For paired expression analysis (pairs, n=4), single primary leukemic cells and normal lymphoblasts were isolated from the same patients with pediatric ALL. The platform was GPL96 (HG-U133A) using the Affymetrix Human Genome U133A Array (Affymetrix Inc., Santa Clara, California, USA). The annotation of the probes from Affymetrix, including all information of Affymetrix ATH1 chip (25K array) was downloaded, as well as the primary files.

Data preprocessing and differential analysis

According to the annotation platform, a cohort of 22,283 probes was mapped to the corresponding genes. For each gene, the average value of multiple probes corresponding to a single gene was calculated to obtain gene expression value, which was then transformed into log2 ratio (13). Consequently, a total of 20,967 gene expression levels were obtained. DEGs were screened for each sample using Linear Models for Microarray Analysis (limma) package of R language (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (14), in which the Benjamini and Hochberg method was applied for multiple testing correction (15). False discovery rate <0.05 and |log fold change >1 were set as the threshold values for DEG selection.

Clustering analysis of DEGs

Biclustering analysis of gene expression data directly identified whether the DEGs that were screened had specificity for the samples (16,17). In the present study, pheatmap package of R language (http://cran.r-project.org/web/packages/pheatmap/index.html) was used for biclustering analysis of DEGs on the basis of Euclidean distance method (18).

Functional enrichment analysis of DEGs

Database for annotation, visualization, and integrated discovery (DAVID; http://david.abcc.ncifcrf.gov/) provided a set of data-mining tools that systematically combined functionally descriptive data with intuitive graphical displays (19). Therefore, DAVID tool was used for the Gene Ontology (GO) functional enrichment analysis of upregulated and downregulated DEGs, based on hypergeometric distribution algorithm. P<0.05 was considered to indicate a statistically significant difference.

Interaction analysis of DEGs

To investigate the association between DEGs at the protein level, the Search Tool for the Retrieval of Interacting Genes (String, version 9.1; http://string-db.org/) (20) was used to screen the PPIs with the threshold of combined score >0.4.

Screening of small-molecule drugs

The upregulated and downregulated DEGs were mapped onto the connectivity map (CMAP) database to obtain the small-molecule drugs targeting the selected DEGs with the criterion of correlation score >0.8 (21). Thereafter, combined with the PPI information, an integrated network between DEGs and small molecular drugs was constructed, which was visualized using the Cytoscape software (version 3.2.1, http://cytoscape.org/) (22).

Results

Screened DEGs between ALL samples and control samples

Following the preprocessing and normalization for gene expression profile data, the medians of the ALL samples were similar (Fig. 1), indicating the validity of these results for the subsequent analysis. Based on the aforementioned criteria, a set of 116 genes was identified as DEGs consisting of 56 downregulated and 60 upregulated genes.

Figure 1. Box diagram of normalization of samples included in the current study. Blue boxes represent normal samples, and pink boxes represent ALL samples (primary). The middle black line represents the median of expression value of each sample. ALL, acute lymphoblastic leukemia; GSE, Gene Expression Omnibus DataSet; GPL, Gene Expression Omnibus platform; GSM, Gene Expression Omnibus.

Clustering analysis of the DEGs

Using the pheatmap package of R language, the cluster heat map of the DEG expression was presented in Fig. 2. By performing biclustering analysis, it was identified that gene expression in pediatric ALL samples and control samples were distinct from each other suggesting that the DEGs that were selected were significant and reliable for distinguishing diseased and control samples.

Figure 2. Heat map of clustering analysis of differentially expressed genes. The color change from blue to orange represents the expression value change from downregulation (blue) to upregulation (red). The sample numbers are presented below, including 7 disease samples (primary) and 4 normal samples.

Altered functions of the DEGs

The upregulated and downregulated DEGs were imported to the DAVID database to perform functional enrichment analysis. Accordingly, upregulated DEGs were established to be significantly enriched in GO functionality, including ‘mesenchymal cell differentiation’, ‘mesenchymal cell’ and ‘mesenchymal development’ The upregulated DEGs included bone marrow morphogenetic protein 2 (BMP2), semaphoring 3F (SEMA3F) and ephrin B1 (EFNB1) (Table I). By contrast, the downregulated DEGs were associated with the terms ‘cell cycle phase’, ‘cell cycle process’ and ‘cell cycle’. The downregulated genes included marker of proliferation Ki-67 (MIK67), cyclin F (CCNF) and nucleolar and spindle associated protein (NUSAP1) (Table II).

Table I. Significantly enriched Gene Ontology terms of upregulated genes.

Table I.

Significantly enriched Gene Ontology terms of upregulated genes.

Gene ontology term	Count	P-value	Genes
0048762~mesenchymal cell differentiation	3	1.35×10−2	BMP2, SEMA3F, EFNB1
0014031~mesenchymal cell development	3	1.35×10−2	BMP2, SEMA3F, EFNB1
0060485~mesenchyme development	3	1.40×10−2	BMP2, SEMA3F, EFNB1
0045765~regulation of angiogenesis	3	2.02×10−2	HHEX, NPR1, EPHB4
0000902~cell morphogenesis	5	3.46×10−2	BMP2, LST1, EFNB1, L1CAM, POU4F1
0006469~negative regulation of protein kinase activity	3	3.67×10−2	SH3BP5, SPRY1, DUSP6
0033673~negative regulation of kinase activity	3	3.90×10−2	SH3BP5, SPRY1, DUSP6
0051348~negative regulation of transferase activity	3	4.38×10−2	SH3BP5, SPRY1, DUSP6
0032989~cellular component morphogenesis	5	4.84×10−2	BMP2, LST1, EFNB1, L1CAM, POU4F1

[i] Count, number of genes enriched in the term.

Table II. Significantly enriched Gene Ontology terms of downregulated genes.

Table II.

Significantly enriched Gene Ontology terms of downregulated genes.

Gene ontology term	Count	P-value	Genes
0022403~cell cycle phase	17	1.31×10−10	CCNF, MKI67, CENPE, NUSAP1, CDC20
0022402~cell cycle process	18	1.06×10−9	MKI67, CCNF, NUSAP1, ESPL1, CENPE
0007049~cell cycle	20	1.18×10−9	MKI67, CCNF, NUSAP1, CDC20, CENPE
0000279~M phase	15	1.78×10−9	KIF11, MKI67, CCNF, NUSAP1, CDC20
0000280~nuclear division	13	5.26×10−9	KIF11, CCNF, NUSAP1, CENPE, CDC20
0007067~mitosis	13	5.26×10−9	KIF11, CCNF, NUSAP1, CENPE, CDC20
0000087~M phase of mitotic cell cycle	13	6.51×10−9	KIF11, CCNF, NUSAP1, CENPE, CDC20
0048285~organelle fission	13	8.44×10−9	KIF11, CCNF, NUSAP1, CENPE, CDC20
0000278~mitotic cell cycle	15	8.72×10−9	KIF11, CCNF, NUSAP1, CDC20, CENPE
0051301~cell division	13	1.63×10−7	KIF11, CCNF, NUSAP1, CENPE, CDC20
0007059~chromosome segregation	6	9.05×10−3	SPC25, NUSAP1, BIRC5, CENPE, ESPL1, PTTG1
0000226~microtubule cytoskeleton organization	7	1.16×10−2	SPC25, CAV1, KIF11, NUSAP1, ESPL1, UBE2C, TACC3
0007051~spindle organization	5	2.05×10−2	SPC25, KIF11, ESPL1, UBE2C, TACC3
0007017~microtubule-based process	8	2.43×10−2	SPC25, CAV1, KIF11, NUSAP1, CENPE, ESPL1, UBE2C, TACC3

[i] Count, number of genes enriched in the term.

Integrated network between small-molecule drugs and targeted DEGs

Following the mapping of DEGs into the String database, a subset of 169 PPIs were identified. In addition, the upregulated and downregulated genes were mapped onto the CMAP database to identify a list of small-molecule drugs that meet the predefined criterion (Table III). Quinostatin and carbimazole achieved high enrichment scores of 0.904 and −0.859, respectively. These findings suggest that theses drugs were associated with ALL and may be potentially effective for the treatment of pediatric ALL.

Table III. Small molecular drugs significantly associated with differentially expressed genes.

Table III.

Small molecular drugs significantly associated with differentially expressed genes.

Cmap name	Enrichment score	P-value
Carbimazole	−0.859	0.00561
Nadolol	−0.854	0.00084
Prestwick-691	−0.845	0.00741
Trihexyphenidyl	−0.829	0.01002
Pargyline	0.833	0.00115
Blebbistatin	0.858	0.04094
Quinostatin	0.904	0.01889

[i] Cmap, connectivity map.

Combining the PPIs with the DEGs and small-molecule drugs, an integrated network was constructed between associated DEGs and the small-molecule drugs (carbimazole and quinostatin) and associated DEGs (Fig. 3). In this integrated network, the target genes of carbimazole and quinostatin were predicted, including signal transducer and activator of transcription 3 (STAT3), nucleolar and spindle associated protein 1 (NUSAP1) and cell division cycle 20 (CDC20).

Figure 3. The integrated network between the two small-molecule drugs (quinostatin and carbimazole) and differentially expressed genes. Red spots indicate significantly upregulated genes and green spots indicate downregulated genes. Orange triangles represent the small-molecule drugs.

Discussion

In the present study, 116 DEGs were identified in pediatric ALL based on gene expression profile data, including 56 downregulated and 60 upregulated genes. Functional enrichment analysis demonstrated that upregulated DEGs, including BMP2, SEMA3F and EFNB1, were significantly associated with mesenchymal cell differentiation and development processes, whilst downregulated DEGs, including MKI67, CCNF and NUSAP1, were enriched for cell cycle processes.

In a previous study, mesenchymal cells were demonstrated to regulate the response of ALL cells to asparaginase, which is a major component of ALL therapy (23). Additionally, the differentially methylated genes between the two leukemia subgroups (CEBPAsil and normal CD34+ bone marrow cells) were also associated with the GO term mesenchymal cell development (24), suggesting that mesenchymal cell differentiation and development may have important roles in regulating ALL. BMP2 was reported to be associated with the cardiac cushion epithelial-mesenchymal transition and the initiation of chondrogenic lineage development of mesenchymal stem cells (25,26). In T-cell ALL, BMP2 was established to be an important regulator by acting upstream of the NOTCH1 signaling pathway (27). The SEMA3F encoded protein is primarily involved in signaling transduction. It has also been demonstrated that SEMA3F is expressed in T-cell ALL, and it also functions in T-cell ALL (28). EFNB1 is another gene expressed in pediatric T-cell ALL (29). Although T-cell and B-cell leukomonocytes mature in the thymus and bone marrow, respectively, immune responses are implicated in both T-cell and B-cell ALL (30,31). Therefore, the regulation of these genes in B-cell ALL may be similar with the gene regulation in T-cell ALL. Considering that the 3 upregulated genes in ALL samples were significantly enriched in the mesenchymal development process, it may be inferred that these genes are involved in the mesenchymal development process during pediatric B-cell ALL progression.

However, the downregulated genes, including MKI67, NUSAP1 and CCNF, were significantly associated with the cell cycle. Cell cycle (including G0/G1 arrest) is a process, which is typically dysregulated in ALL, (32). ALL is a result of the clonal expansion of hematopoietic progenitors that have undergone malignant transformation at distinct stages of differentiation (33). The B-cell lineage ALL is derived from B cell precursors. The Ikaros-mediated cell cycle was determined as the endpoint in the tumor suppression pathway that was modulated by the pre-B cell receptor (34). Aberrations of genes that are involved in cell cycle control, including the cytokine receptor-like factor 2 gene, have been implicated in the progression and relapse of B-cell ALL (35). The majority of cell cycle genes were deregulated in the Bruton tyrosine kinase (BTK)C481S mantle cell lymphoma (MCL) cells, especially the expression of MKI67 was remarkably decreased in the bone marrow at the point of relapse (36). A previous study demonstrated that microRNA-519d is a tumor suppressor that targets MKI67 in the hepatocellular carcinoma cell line (37), which may account for the downregulation of MKI67. NUSAP1 is a protein that is key for spindle movement (38). MRD is considered to be an important predictor of relapse in ALL (39), and reduced expression of NUSAP1 was identified in patients with MRD (40), which is consistent with the finding of downregulated NUSAP1 in the current study. CCNF encodes a member of the cyclin family that serves important roles in cell cycle transitions (41). The downregulation of CCNF was also established in B-cell ALL (42). NUSAP1 and CCNF were associated with altered cell cycle in acute myeloid leukemia (43), therefore providing further evidence of the role of the two genes in cell cycle modulation. These results collectively suggest that the 3 genes (MKI67, NUSAP1 and CCNF) mediate cell cycle processes and may have regulatory roles in the B-cell ALL progression.

The present study also identified a number of important small molecule drugs, including quinostatin and carbimazole. NOTCH1 is a member of the Notch family that has vital roles in a number of cellular processes. The activation of NOTCH1 signaling promotes important functions in the pathogenesis of T-cell ALL (44). NOTCH1 inhibits cell apoptosis via the transcriptional suppressor, hairy and enhancer of split 1 (HES1)-mediated repression downstream of NOTCH signaling in T-cell ALL (45). The small molecule perhexiline was identified as a HES1 modulator drug with potent anti-leukemic effects (45). Notably based on the information in CMAP, quinostatin was identified as another important small molecule drug that may target HES1 with a high enrichment score (0.983) (45), suggesting a potential inhibitory role of quinostatin in ALL. Additionally, quinostatin was reported to potently inhibit the mammalian target of rapamycin signaling network by directly targeting lipid-kinase activity (46). Combining these findings with the results of the present study that quinostatin was associated with DEGs (enrichment score, 0.904), quinostatin may be a potential therapeutic drug for treating ALL. However, further experimental validations are required to confirm this prediction.

Carbimazole has been widely used in the treatment of ALL (47,48). In the present study, it was demonstrated that quinostatin and carbimazole share a number of common target genes, including STAT3, NUSAP1 and CDC20. Excluding the aforementioned functions, STAT3 was revealed to mediate an oncogenic association with two characterized transcription factors, Ets-related protein Tel1 (TEL) and Runt related transcription factor 1 (AML1), in t(12;21) translocation, which has been reported to be the most frequent chromosomal abnormality in pediatric leukemia (49). NUSAP is an important mitotic regulator, and the activity of NUSAP is vital for a number of cellular processes during mitosis (50). CDC20, as a mitotic complex/cyclosome activator, directly interacts with anaphase-promoting complex, which leads to the promotion of onset of anaphase and mitotic exit via the ubiquitination of securin and cyclin B1 (51,52). Taken together, it may be hypothesized that the two small molecule drugs, particularly quinostatin, may be an effective agent for the management of ALL, as it targets STAT3, NUSAP1 and CDC20 genes. However, further experiments are required to validate these targets. Although the current study performed a comprehensive bioinformatics analysis, it should be noted that the data were derived from a previously published study with a small sample size, and the present study did not consider the subgroups of ALL. All the conclusions of the current study require further experimental validation.

In conclusion, by using a bioinformatic methodology, the current study identified that the dysregulated biological processes (including mesenchymal cell differentiation and development) that were mediated by the upregulated genes and the cell cycle processes that were mediated by downregulated genes may have significant roles in the pathogenesis and prognosis of pediatric B-cell ALL. In addition, quinostatin may be a potent novel agent for the management of pediatric B-cell ALL as it targets a large number of DEGs, including STAT3, NUSAP1 and CDC20. However, the results obtained in the present study require further validation.

References