Gene module analysis of juvenile myelomonocytic leukemia and screening of anticancer drugs.

Juvenile myelomonocytic leukemia (JMML) is a rare but severe primary hemopoietic system tumor of childhood, most frequent in children 4 years and younger. There are currently no specific anticancer therapies targeting JMML, and the underlying gene expression changes have not been revealed. To define molecular targets and possible biomarkers for early diagnosis, optimal treatment, and prognosis, we conducted microarray data analysis using the Gene Expression Omnibus, and constructed protein‑protein interaction networks of all differentially expressed genes. Modular bioinformatics analysis revealed four core functional modules for JMML. We analyzed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functions associated with these modules. Using the CMap database, nine potential anticancer drugs were identified that modulate expression levels of many JMML‑associated genes. In addition, we identified possible miRNAs and transcription factors regulating these differentially expressed genes. This study defines a new research strategy for developing JMML‑targeted chemotherapies.

Introduction

Juvenile myelomonocytic leukemia (JMML) is an infrequent but aggressive hematological system tumor of infancy and childhood with poor prognosis. Signs and symptoms include fever, thrombocytopenia, hemorrhage, lymphadenopathy, high fetal hemoglobin levels and progressive hematologic anemia, hepatosplenomegaly, and clonal proliferation of myelomonocytic cells (1–3). The most effective treatment currently available is allogeneic hematopoietic stem cell transplantation (HSCT). However, the overall survival at 5 years after transplant is only 52–63% due to treatment-related toxicity and frequent disease relapse (3–7). In addition, some patients lack a suitable human leukocyte antigen-matched donor, thus novel chemotherapeutic regimens are needed to further improve outcome.

Malignant tumors are a major cause of mortality in children, and traditional cancer treatments such as chemotherapy, radiotherapy and surgery carry severe side effects or do not markedly improve prognosis. As an aggressive myeloproliferative neoplasm, response to conventional chemotherapy is weak, thus most JMML patients require early allotransplantation (6,7). Alternatively, gene module analysis of JMML may facilitate the screening of anticancer drugs targeting specific anomalies in gene expression. In the present study, we performed gene module analysis to identify gene expression signatures of JMML and to illuminate the functions of these altered genes and protein networks in disease pathogenesis.

While there are voluminous studies on anticancer drug pharmacokinetics and efficacy against solid tumors, there have been few studies focusing on JMML. Drugs that target JMML-associated genes can be an effective anticancer approach, either alone or in combination with traditional modalities. Modular analysis of bioinformatic data is emerging as a valuable strategy for systematic and comprehensive analysis of regulatory signaling pathways relevant to disease etiology and drug treatment effects (8,9). In the present study, we used modular analysis bioinformatics methods to search for differentially expressed genes and their key functional group(s) and networks in JMML.

Genes differentially expressed in leukemic cells compared to normal cells were classified into upregulated and downregulated groups. These differentially expressed genes were then compared to differentially expressed genes in cells treated with small molecules from the CMap database by converting them into probe setters based on the HG-U133A platform. Genes were then assigned enrichment values between-1 and 1 according to whether specific small molecules stimulated the expression state in leukemia cells (closer to 1) or normal cells (closer to −1). Based on this analysis, we identified drugs that regulate key nodes of JMML-associated regulatory networks. The large number of genes and drugs identified provide many possible targets and effective chemotherapies for JMML.

Materials and methods

Screening for differentially expressed genes in JMML

We downloaded JMML microarray data from the GEO platform (https://www.ncbi.nlm.nih.gov/geo/) [data code GSE71449 (human)] (10).

The criteria for data selection were gene chip expression profiling data accrued over the last three years (since 2014) with clear disease and normal control groups in the data set. The disease queried was ‘Juvenile myelomonocytic leukemia’. The largest number of GSE71449 samples were obtained by such a standard. While our analysis is indeed limited to all disease-related data sets in the database, further analysis of other data sets, including miRNA chips, is already under way as no single data set represents all genes associated with a disease. The mechanisms of any disease are only revealed by direct experimental study, but these gene chip results may identify genes with greater probability of involvement. The single-chip DEG analysis can only reveal that in the current data set, these differentially expressed genes (DEGs) are most likely to be related to disease.

Array Express was retrieved, and the results retrieved with the same search criteria coincided with those in GEO. There was no larger sample size under the same screening criteria. GEO data sets are included in the EBI database, thus there is no mention of it.

The database provides tissue sample results from 44 JMML patients and seven healthy donors obtained using the Agilent-041648 CMGG Human V1.1 60k [Probe Name Version] microarray platform. We used the R software bagaffy (11) (Version1.50.0, http://www.bioconductor.org/packages/release/bioc/html/affy.html) to read the downloaded microarray matrix data and the robust multi-array average method (12,13) to perform standardized data preprocessing (including Background correction, Normalization, and Expression calculation). We performed annotation for probes in the platform annotation file and deleted probes that did not match the gene symbols used. When the same gene was reflected by different probes, we adopted the average as the final expression value.

For the two data groups (JMML and controls), we screened for differentially expressed genes using the R bag Limma microarray analysis package (14) and adopted the BH-corrected T test to identify genes significantly over- or under expressed. For every differentially expressed gene, we assumed a threshold P<0.05 and log|FC| >0.58. Differentially expressed genes were selected based on a threshold |logFC| >0.58, corresponding to a fold change >1.5 or ≤1.5 (15). We then used the common enrichment analysis tool DAVID (16) (version 6.8) and subjected genes with enrichment ≥2 and significance threshold of hypergeometric test P<0.05 to Gene Ontology (GO) (17) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) (18) pathway analyses.

For observing the function of differentially expressed genes intuitively, we used the ClueGO (19) plug-in (version 2.2.6, http://apps.cytoscape.org/apps/ClueGO) of Cytoscape software (20) and constructed separate GO BP and KEGG pathway cross-linked enrichment graphs for upregulated and downregulated genes using P<0.05 as the significantly enriched threshold.

Building and analyzing modules of differentially expressed genes

We used the STRING (21) database bank to predict possible protein-protein interaction (PPI) networks for all differentially expressed proteins and genes with parameter PPI score 0.4 (medium confidence). Networks with average interaction strengths at protein nodes were constructed using Cytoscape.

Interaction networks may contain nodes in which the interactions are closed among the constituent genes. These nodes thus represent distinct biological processes. There are many methods to define nodes by PPI cluster analysis but MCODE (22) was chosen because it is the most common. We applied MCODE to calculate the scores of every node, with higher scores indicating a greater degree of separation from other nodes and stronger association with specific processes. We choose nodes with score ≥5 and node number ≥5 for subsequent GO and KEGG pathway analyses.

Screening for drugs that regulate functional modules

The CMap database (23,24) stores the genome-wide expression profiles of human cells treated with various active small molecules. In total, CMap contains data from 6,100 small molecule interference experiments (with normal control groups) using 1,309 small molecules, for a total of 7,056 different expression profiles.

We analyzed gene expression differences between normal cells and leukemia cells and then compared the responses of differentially expressed genes to identify those with similar or opposite effects (upregulation vs. downregulation) on normal cells. We divided the genes differentially expressed between normal and leukemia cells into upregulated and downregulated subgroups. Then we transformed the HG-U133A platform probe set results and compared them to the set of differentially expressed genes under small molecule treatment from the CMap database to obtain enrichment values. In this case, enrichment values range from-1 and 1, with those closer to 1 influenced to a greater extent by small molecules in normal cells and those with values closer to −1 influenced to a greater extent by small molecules in leukemia cells.

For key modular gene screening, we constructed connectivity maps and used the gene expression differences in human cells treated with small molecules to identify drugs affecting the expression levels of genes associated with the disease.

Identification of possible miRNAs regulating target genes

We identified potential miRNA targeting genes within PPIs using Webgestalt (http://www.webgestalt.org/) (25) tools and overrepresentation enrichment analysis methods with a threshold P<0.05 by hypergeometric tests and BH correction for enriched gene number and count ≥5.

Identification of possible transcription factors regulating target genes

Based on the transcription factor-regulated network data in the ITFP database and TRANSFAC bank, we searched for transcription factors regulating differentially expressed genes and further screened the differentially expressed target genes regulated by these transcription factors for network integration. We performed network integration for the obtained TF-Target, miRNA-Target, and PPI network and constructed integration networks using Cytoscape.

Statistical analysis

t-tests were used to screen for differentially expressed genes and to evaluate the enrichment obtained by GO_BP- and KEGG-based methods. A P<0.01 was considered statistically significant. With respect to each herbal component, Fisher's exact test was used to evaluate the genes in modules. Based on the above data, we set P<0.01 for each component.

Results

We downloaded the JMML expression data from the GEO database, and identified 700 upregulated and 737 downregulated genes at threshold levels of P<0.05 and log|FC| >0.58. From these genes, we constructed heat maps (Fig. 1). We then performed GO and KEGG pathway enrichment analyses for upregulated and downregulated mRNA muster, and according to the screening threshold value obtained GO terms and KEGG pathways including these differentially expressed genes. Fig. 2A shows the top five GO Biological Process (BP) terms and KEGG pathways including upregulated mRNAs and Fig. 2B shows these same results for downregulated mRNAs.

Figure 1.

We downloaded data from the GEO database (experimental number GSE71449) and obtained 700 upregulated and 737 downregulated genes using threshold P<0.05 and log|FC| >0.58 was threshold. We constructed the heat map for the obtained differentially expressed genes through screening.

Figure 2.

GO and KEGG pathway enrichment analysis of the upregulated and downregulated mRNA sets. According to the screening threshold, the GO terms and KEGG pathways involving these differentially expressed genes were obtained. (A) Top five GO BP terms and KEGG pathways for upregulated mRNAs. (B) Top five GO BP terms and KEGG pathways for downregulated mRNAs. GO, Gene Ontology; BP, Biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene.

We used clueGO to perform GO and KEGG pathway analysis for all upregulated and downregulated differentially expressed genes (Fig. 3A and B). ClueGO provides kappa coefficients that can be used to divide the pathways into functional groups. In the figure, different colors indicate pathway enrichment results and correlations between two pathways are reflected by a connecting line. The sizes of nodes reflect the P-value. Fusion of enrichment results was performed when different terms were enriched for the same gene. In the figure, GO terms are indicated by different colors. Thus, one color indicates a functional group. The size of the term nodes depends on the threshold P-value, increasing as P gets smaller.

Figure 3.

The enrichment and crosslink figure of GO BP and KEGG pathways. We used clueGO to perform GO and KEGG pathway analysis for all upregulated (A) differentially expressed genes. We used clueGO to perform GO and KEGG pathway analysis for all downregulated (B) differentially expressed genes. ClueGO revealed correlations among channels by calculating the kappa coefficients. We divided the path into several functional groups according to the kappa coefficients to reveal several networks. Different colors represent different paths. When the two paths have relevance, they are connected by a straight line. The size of the node represents the P-value of the GO term, with smaller P-values indicated by larger nodes. When the same gene is enriched to different terms, the term is plotted in a variety of colors. The same color term represents a functional group. The minimum P-value of the same functional group is termed as the maximum node. GO, Gene Ontology; BP, Biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Key genes and PPI networks constructed from differentially expressed genes

Combined with the PPI database, we constructed the network figures for differentially expressed genes (Fig. 4). The PPI network of differentially expressed genes contained 908 nodes and 5,053 related pairs. The relationships among the nodes are closed. Table I lists the top 10 highest degree nodes obtained from the PPI network and Fig. 4 shows the PPI network constructed from the differentially expressed genes, with pink nodes indicating upregulated and blue nodes the downregulated genes.

Figure 4.

Combining the PPI database, the network diagram was constructed for the differentially expressed genes. The PPI network of differentially expressed genes contained 908 nodes and 5,053 related pairs. There is a close relationship between the nodes. Table I provides the degree node for the PPI network scoring top ten. The PPI network based on differential gene construction is shown. The pink nodes show the upregulated genes and the blue nodes the downregulated genes. PPI, protein-protein interaction.

Table I.

The top 10 highest degree nodes obtained from the PPI network.

Gene	Degree	Betweenness	Closeness
MAPK1	283.0	142874.9	0.54246414
CD4	218.0	93239.1	0.49266702
OAS3	212.0	110066.46	0.48373333
CD79A	173.0	59164.633	0.47862798
NFKBIA	157.0	37323.453	0.48424986
BCL6	154.0	42434.684	0.4748691
EGR2	138.0	33878.51	0.4512438
PLEK	134.0	53153.11	0.45509282
TRIB3	115.0	24556.303	0.44591936
RASL11A	112.0	16269.245	0.44923228

PPI, protein-protein interaction.

Modular analysis

We obtained a total of 10 modules from the PPI network, of which four had scores >3 (Fig. 5). We also performed functional enrichment analysis for the genes of modules. Fig. 6 shows the top five GO BP terms and KEGG pathways in module 1, the four KEGG pathways and top five GO BP terms in module 2, the four GO BP terms in module 3 (no KEGG pathways were enriched), and the top five GO BP terms and KEGG pathways in module 4.

Figure 5.

We obtained a total of 10 modules from the PPI network, including 4 modules (A and B) with scores >3. Upregulated genes are in red and downregulated genes in blue. We obtained a total of 10 modules from the PPI network, including 4 modules (C and D) with scores >3. Upregulated genes are in red and downregulated genes in blue.

Figure 6.

Functions of the genes in each module revealed by GO BP and KEGG analyses. (S) The top five GO BP terms and KEGG pathway results for module 1. (B) The four KEGG pathways and top five GO BP terms in module 2. (C) The four GO BP terms in module 3 (no enrichment to KEGG pathways). (D) Top five GO BP terms and KEGG pathway results in module 4. Pathways of greatest significance in each module were Regulation of cell proliferation in module 1, Pathways in cancer, Transcriptional dysregulation in cancer, and Signal transduction in module 2, Transcription, DNA-templated in module 3, and Immune response in module 4. GO, Gene Ontology; BP, Biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes.

For module 1, pathway analysis indicated that differentially expressed genes were most strongly (lowest P-value) related to Regulation of cell proliferation. For module 2, differentially expressed genes were most strongly related to Pathways in cancer. In module 3, the strongest relationship was with Transcription, DNA-templated, and in module 4 the strongest relationship was with Immune response.

Small drug molecules related to diseases

We performed small molecule drug analysis for the 4 modules using the CMap database and threshold P<0.05. The regulatory effects on each module are shown in Fig. 7.

Figure 7.

Small drugs affecting the 4 modules (A-D) according to the CMap database. Shown are the top 10 results with threshold P<0.05. Drug were found that affected each module.

Statistical relationships between modular genes and small drugs

We used local perl script and performed Pubmed searches for the analyzed modular genes and the corresponding small molecules, and constructed the network figures shown in Fig. 8. The P-values were considered in descending order, and the drugs with the 10 smallest P-values were used for subsequent analyses. A total of 40 drugs were identified as possible regulators of the 4 modules. To obtain novel anticancer drugs, we conducted a literature review of the cancer-related research on the 24 drugs with greatest effects on gene expression. In cancer treatment, the number of publications reflects the extent of research on a given drug's potential efficacy. Multiple drugs have been researched extensively and already certified to be effective anticancer drugs, such as etoposide and phosphonothreonine (the number of relevant publications was >100). Others, however, have not been fully studied for their anticancer efficacies, such as disulfiram, ursolic acid, miconazole, thioridazine, loperamide, monastrol (the number of relevant publications was ≤10), and others, such as nadolol, tetracaine and levomepromazine have never been studied with respect to cancer. Publication numbers for drug are shown in Fig. 9.

Figure 8.

We used local perl script and performed Pubmed searches for each small molecule, and constructed network figures (Pink genes in module 1, red genes in module 2, purple genes in module 3, and yellow genes in module 4). Arrows represent small molecule drug regulatory relationships. The line thickness represents the volume of existing research on these small molecule drugs and genes, with thicker lines indicating more extensive research (more publications). Blue V indicates small molecule drugs.

Figure 9.

To evaluate potential efficacy, we conducted a literature appraisal of 24 small molecules. The number of publications is indicative of the extent of study on cancer. Those with fewer publications may be novel anticancer drugs. Of those identified, many are confirmed to be anticancer drugs, for example etoposide and phosphonothreonine (>100 relevant publications), while others have not been studied extensively for anticancer efficacy, for example disulfiram, ursolic acid, miconazole, thioridazine, loperamide and monastrol (≤10 relevant publications). A few, such as nadolol, tetracaine, and levomepromazine, have never been examined for anticancer properties.

Identification of candidate miRNAs regulating differentially expressed genes

We identified candidate miRNAs that may regulate the differentially expressed genes using a screening threshold of P<0.05. For upregulated genes, we identified 10 candidate miRNAs (Table IIA) and for the downregulated genes four candidate miRNAs (Table IIB).

Table II.

Identification of candidate miRNAs regulating differentially expressed genes (DEGs).

A, miRNAs for upregulated DEGs
miRNA	P-value
AGTCAGC, MIR-345	0.002846737
CAGGTCC, MIR-492	0.003876902
ACCAATC, MIR-509	0.011602863
CCAGGTT, MIR-490	0.017555181
TGCACTG, MIR-148A, MIR-152, MIR-148B	0.023681008
TGAGATT, MIR-216	0.030489785
TCCAGAT, MIR-516-5P	0.033806763
AGTCTTA, MIR-499	0.034998955
GTGTCAA, MIR-514	0.037442107
AAAGGGA, MIR-204, MIR-211	0.040666691
B, miRNAs for downregulated DEGs
miRNA	P-value
GGCACTT, MIR-519E	0.002921036
TCCGTCC, MIR-184	0.037734308
GTAAGAT, MIR-200A	0.040359704
AGCATTA, MIR-155	0.04249545

Identification of candidate transcription factors regulating differentially expressed genes

For the upregulated differentially expressed genes, we identified 22 differentially expressed transcription factors (TFs). For the downregulated differentially expressed genes, we identified 19 differentially expressed TFs. The results are summarized in Table III.

Table III.

Identification of candidate transcription factors (TFs) regulating differentially expressed genes (DEGs).

TFs for upregulated DEGs	TFs for downregulated DEGs
ABTB2, AHCTF1, BATF, CORO7, ERG, FGD5, FOXN2, IFIH1, IFIT1, IFIT2, IFIT3, ING1, INPP5A, LASP1, PARP14, PLEK, PTBP2, SAMD9, STAT1, TANK, TTC7B, WDR1	CCNE1, CDC20, CKAP5, FAM129B, KRT7, LGR4, MKL2, NFIA, REXO1, SMARCA4, SPEN, TLE2, TRIM10, TRIM27, TTC16, ZBTB4, ZNF497, ZNF750, ZNF771

Discussion

Currently the only curative treatment option for JMML is hematopoietic stem cell transplantation (HSCT). However, disease recurrence remains a major cause of treatment failure (7). Clinical symptoms are caused by hematopoietic insufficiency and excessive proliferation of leukemic monocytes and granulocytes, leading to hepatosplenomegaly, lymphadenopathy, skin rash and respiratory failure (3,7,26). A serious obstacle to the research of JMML is the lack of suitable experimental models, impeding the development and pre-clinical evaluation of novel therapeutic approaches. Primary JMML leukemia cells cannot be maintained in culture as they differentiate and become apoptotic (27), while an immortalized cell line derived from JMML cells has not yet been successfully established (28). The generation of induced pluripotent stem cell lines originating from JMML cells was reported, but conceptually such systems are limited by their artificial nature and the risk of further transformation during reprogramming (29). Therefore, it is important to identify tumor-promoting processes, such as dysregulation of specific genes and networks, and develop targeted therapies.

In the present study, we combined the genome-wide gene chip information already available from JMML tissues and used modular analysis bioinformatics procedures to identify molecular targets and small molecule compounds that affect their expression. For this purpose, modular analysis is a valuable and well established bioinformatic method that is useful for analysis of both large-scale protein networks and single proteins and genes. Moreover, this method has been used to study the development and treatment of multiple diseases. After the enrichment analysis, we found that genes related to the hematopoietic cell lineage and to regulation of immune responses were significantly downregulated, while chemokine signaling pathway and inflammatory response genes were significantly upregulated compared to healthy tissue.

Hematopoietic progenitors produce a myriad of diverse lineages, including progenitors with lymphoid, myeloid, and erythroid potential, prior to hematopoietic stem cells. Previous pathway analyses revealed that genes related to hematopoietic cell lineages were significantly downregulated in JMML (30–34), and we further found that genes related to regulation of immune responses were also significantly downregulated. Anticancer treatment leads to dense immune infiltrate, including many CD8+ and CD4+ cells, concomitant with tumor regression, at both treated and untreated lesions, consistent with generation of a tumor-specific systemic immune response (35). Several studies have indicated the synergistic potential of various immunomodulatory agents (36–39). Therapies that enhance or impede immune responses are essential but optimal timing and the administration route are important detemising results, such as increased survival and immune responses, leading to complete tumor regression in some cases (40–42).

Rminants of efficacy. Preclinical studies combining these agents have shown proThe PI3K/Akt and MAPK signaling pathways have been shown to mediate chemokine-induced migration of multiple cell types. Modular analysis suggested that the Chemokine signaling pathway (module 1) was strongly related to JMML morbidity. In addition to regulation of cell proliferation, differentiation, invasion, and inflammation, chemokines are widely involved in the regulation of cancer. Analysis of module 1 revealed that regulation of cell proliferation was significantly increased and strongly related to JMML occurrence. Further, analysis of all 4 functional modules revealed complex interactions among a number of genes involved in Chemokine signaling pathway and Regulation of cell proliferation. Chemokines are small proteins expressed in response to injury or infection and during normal immune surveillance. Chemokines are involved in leukocyte trafficking and regulate tumor metastasis, proliferation, differentiation and angiogenesis (43–46).

Module 2 contained genes primarily involved in Pathways in cancer, Transcriptional dysregulation in cancer, and Signal transduction, suggesting that module 2 contains clusters of key tumor-promoting genes. Accumulation of driver somatic alterations in genes modifies critical cellular processes leading to cancer (47,48). In recent years, the catalog of driver genes known to take part in the development of malignancies has expanded due to whole-exome and whole-genome analyses of large tumor sets by large international consortia (49,50). In addition to upregulation of oncogenes, downregulation of anticancer proteins or dysregulation of subcellular localization may also contribute to JMML.

Cancer cells scatter from tumors in hoards and build new tumors in distant tissues and organs (51,52). In over 90% of fatal cancers, metastasis is the cause of death. Tumors contain a heterogeneous population of cells in an organized hierarchy akin to normal tissue. Tumor progression, invasiveness, and self-renewal are attributes of smaller subfractions of cancer cells. Tumor initiation by disseminated cancer cells relies on their ability to self-renew and initiate metastatic tumors.

Genes in module 3 were primarily involved in Response to hypoxia, Negative regulation of transcription from RNA polymerase II promoter, Regulation of transcription, DNA-templated, and Transcription, DNA-templated. We found that among these pathways, transcription, DNA-templated was most strongly related to JMML occurrence. The identification and detection of specific nucleic acids (either DNA or RNA) is an enabling technology for forensic analysis (53), recognition of genetic mutations (54) and pathogen identification (55). For many methods, quantification of nucleic acids is important for identifying a given DNA sequence. For example, in many forms of cancer (56), oncogene copy number increases and deletion of tumor-suppressor genes can be found. In addition to measuring changes in gene copy number, quantitative detection of circulating DNAs can be a useful diagnostic tool for identifying cancer (57).

The genes in module 4 are closely associated with extracellular matrix organization and immune responses. Much effort has been devoted to determining how cellular components of the tumor promote cancer development and initiate formation of a niche conducive to growth (58). Alternatively, the importance of non-cellular components of the niche during cancer progression is well documented, particularly extracellular matrix (ECM) organization (59–62). In addition to a stable structure with supportive functions in maintaining tissue morphology, the ECM is a surprisingly dynamic and versatile part of the cell milieu influencing fundamental aspects of cell biology (63). For major developmental processes, the ECM directly or indirectly regulates almost all cellular behaviors (64–67). Indeed, in diseases such as cancer, abnormal ECM dynamics may be a strong determinant of clinical outcome (68). How disruption of ECM dynamics contributes to tumor development is a challenging issue in cancer biology.

The local immune response has also emerged as an important element in the multistep process of cancer development (69). Observations that some tumors arise from chronic inflammation sites have led to speculation of a strong connection between tumor onset and inflammatory pathways. In addition, some tumors are infiltrated by both the innate and adaptive arms of the immune system, thus many different immune cells are present within the tumor microenvironment (70). Elements of both the innate and adaptive immune systems have been reported to act both as pro- or anti-tumorigenic factors depending on the relative balance. The intercellular communication between infiltrating immune cells and cancer cells modulates this immune response so as to positively influence tumor development (71).

We found the module 2 was most strongly related to JMML occurrence. There are likely synergistic interactions among module 2 genes.

Many investigations have attempted to screen for well-tolerated and affordable anti-neoplastic medications, and a myriad of drugs with cytotoxicity against human cancer cell lines have been described (72,73). In our study, we also screened for drugs with curative potential against JMML through modulation of associated gene pathways. In 1999, Mayer et al (74) first identified monastrol, a cell-permeable small molecule with antimitotic activity but without neuronal cytotoxicity. Through inhibition of kinesin Eg5, monastrol can induce the mono-astral conformation of microtubules (75,76). Several subsequent studies have clarified the anti-mitotic mechanisms of monastrol (77–79), but few studies have investigated its anticancer activity (80–82).

In both canine and human cancer patients, the peripherally acting µ-opiate receptor agonist loperamide hydrochloride is recommended as a treatment for chemotherapy-related diarrhea (83,84). Loperamide was shown to dose-dependently induce apoptosis and suppress the proliferation of human liver, lung, bone, and breast cancer cell lines (85). In human cancer cell lines, the mechanism underlying apoptosis induction has not been fully elucidated, although the caspase 3 pathway has been implicated (85) and loperamide will exert anticancer properties at clinically relevant doses. In the clinic, loperamide would be an attractive drug for JMML because of its minimal side effects and low price. Our bioinformatic analysis revealed a potential role for loperamide in JMML treatment. Further studies are needed to establish mechanisms of antitumor activity.

Originally, thioridazine (TDZ) was used as a therapy for psychotic disease (86,87). In addition, it has been used to treat drug-resistant microorganisms (88,89). Recently TDZ was reported to have potent effects on various types of cancer cells, including anti-angiogenesis and apoptosis promotion of breast and ovarian cancers cells (90,91). In addition, TDZ induced cytotoxicity of cervical (92), prostate (93), gastric (94), and pancreatic cancer cells (95). TDZ also has selectivity in leukemia as a dopamine receptor inhibitor (96,97).

A few publications have suggested that miconazole has anticancer effects (98). Miconazole is a common treatment for superficial fungal infection and a prominent systemic antifungal agent. In different human neoplastic cell lines, Wu et al demonstrated that miconazole could induce cell cycle arrest. This growth arrest was dose-dependent and related to the p53 signaling pathway (99). Anticancer effects of miconazole were reported on 4T1 (breast cancer) and 5637 (bladder cancer) cell lines.

Ursolic acid (UA) was found to have a biphasic response against three breast cancer cell lines (100). Another study (101) reported that disulfiram suppressed tumor growth by killing cancer cells and was even effective in combination with DHA. Indeed, UA could act as a substitute for clioquinol.

Compared to treatments targeting individual genes, treatment with agents affecting larger gene groups may have better efficacy as anticancer therapy. The aim of the present research was to identify possible molecular targets for cancer treatment and potential anticancer drugs. We searched for drugs that regulate the essential functional modules of JMML using the DrugBank Small Molecule database, and found several drugs already shown to have anticancer effects, including nadolol, disulfram, ursolic acid, micronazole, thioridazine, loperamide, monastrol, tetracaine, and levomepromazine. In contrast, other drugs identified have never been examined for anticancer efficacy, such as tetracaine, levomepromazine and nadolol. Nonetheless, effects on differentially expressed genes in JMML suggest therapeutic potential.

The nine drugs identified are potential therapies for JMML. Many widely used and studied drugs have never been examined for effects against cancer, such as tetracaine, levomepromazine and nadolol. According to our enrichment analyses, these agents hold potential for JMML treatment. We searched the Pubmed database and found several RCTs on these drugs, but no published clinical RCTs supporting anticancer effects on JMML. In the future, we plan to perform RCTs on these nine drugs examining possible therapeutic efficacy against JMML.

When designing therapeutic regimens for JMML, possible adverse reactions must be considered. It is thus noteworthy that the drugs identified here are in clinical use with well described safety profiles. When tumor cell proliferation reaches the highest activity, drugs interfering with the cell cycle or targeting proliferation pathways are significantly more potent.

MicroRNAs are endogenous 20–25 nucleotide non-coding RNAs found in eukaryotes that regulate gene expression at the level of translation. Mature miRNAs are produced by sequential cleavage of longer primary transcripts. Then the RNA-induced silencing complex RISC induced by RNA is assembled to recognize target mRNAs by base complementary pairing, and the target mRNA is degraded or suppressed according to the degree of complementarity. Recent studies have shown that the silencing complex degrades the target mRNA or inhibits the translation of the target mRNA. Through gene silencing, miRNAs regulate diverse processes including growth, virus defense, hematopoiesis, organ formation, cell propagation, apoptosis, and fat metabolism, among others. Transcription factors (TFs) are proteins with special structures that regulate gene expression at the level of the genome. The transcription initiation of eukaryotes is complex and requires the assistance of multiple protein cofactors. Transcription factors form a transcription initiation complex with RNA polymerase II. Both miRNAs and TFs regulate genes, thus differentially expressed genes (DEGs) are predicted by miRNAs and TFs, which further extends the number of potential targets. Like differentially expressed genes, the prediction of miRNAs and TFs involved in differential expression also provides clues to disease mechanisms.

In conclusion, our study used multiple bioinformatic methods to identify 4 gene modules associated with JMML. In addition, we identified 40 drugs based on the CMap database that can alter expression and function of the 4 modules. We performed Fisher's exact test to precisely screen for drugs with the strongest effects on module regulation. Through our study, we identified nine drugs that have considerable potential as new anticancer drugs. Moreover, the study provides a new research template for future research on JMML-targeted anticancer treatments through detailed analyses of core functions. In the future, we will conduct follow-up studies to verify the anticancer effects of the selected drugs.