Identification of Novel Kinase-Transcription Factor-mRNA-miRNA Regulatory Network in Nasopharyngeal Carcinoma by Bioinformatics Analysis.

PURPOSE: Nasopharyngeal carcinoma (NPC) is one of the most common malignant tumors of the head and neck. This study aimed to investigate the crucial genes and regulatory networks involved in the carcinogenesis of NPC using a bioinformatics approach.
METHODS: Five mRNA and two miRNA expression datasets were downloaded from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) and miRNAs (DEMs) between NPC and normal samples were analyzed using R software. The WebGestalt tool was used for functional enrichment analysis, and protein-protein interaction (PPI) network analysis of DEGs was performed using STRING database. Transcription factors (TFs) were predicted using TRRUST and Transcriptional Regulatory Element Database (TRED). Kinases were identified using X2Kgui. The miRNAs of DEGs were predicted using miRWalk database. A kinase-TF-mRNA-miRNA integrated network was constructed, and hub nodes were selected. The hub genes were validated using NPC datasets from the GEO and Oncomine databases. Finally, candidate small-molecule agents were predicted using CMap.
RESULTS: A total of 122 DEGs and 44 DEMs were identified. DEGs were associated with the immune response, leukocyte activation, endoplasmic reticulum stress in GO analysis, and the NF-κB signaling pathway in KEGG analysis. Four significant modules were identified using PPI network analysis. Subsequently, 26 TFs, 73 kinases, and 2499 miRNAs were predicted. The predicted miRNAs were cross-referenced with DEMs, and seven overlapping miRNAs were selected. In the kinase-TF-mRNA-miRNA integrated network, eight genes (PTGS2, FN1, MMP1, PLAU, MMP3, CD19, BMP2, and PIGR) were identified as hub genes. Hub genes were validated with consistent results, indicating the reliability of our findings. Finally, six candidate small-molecule agents (phenoxybenzamine, luteolin, thioguanosine, reserpine, blebbistatin, and camptothecin) were predicted.
CONCLUSION: We identified DEGs and an NPC regulatory network involving kinases, TFs, mRNAs, and miRNAs, which might provide promising insight into the pathogenesis, treatment, and prognosis of NPC.

Introduction

Nasopharyngeal carcinoma (NPC) is a type of squamous cell carcinoma arising from the nasopharyngeal mucosal lining and is one of the most common malignancies in the head and neck.1 NPC is a relatively uncommon cancer, accounting for only 0.7% of all cancers diagnosed in 2018.2 However, the incidence of NPC is extremely inconsistent all over the world, with more than 70% of new cases occurring in East and Southeast Asia.2,3 Pathologically, most NPCs are characterized as poorly differentiated squamous cell carcinoma and undifferentiated carcinoma.4 NPC is usually located around the ostium of the eustachian tube in the lateral wall of the nasopharynx. The potential risk factors for the progression of NPC include Epstein-Barr virus (EBV) infection, genetic susceptibility, dietary factors, smoking, alcohol consumption, and workplace exposure.5,6 Because of the occult occurrence site and lack of characteristic early symptoms, NPC patients are often diagnosed at an advanced stage. Although great achievements have been made in recent NPC research, the molecular mechanisms underlying the oncogenesis and tumor progression of NPC have not yet been fully understood. Therefore, new targets for the early diagnosis and treatment of NPC are needed, in the hope of improving the prognosis and lowering mortality.

Accumulating evidence has verified that several differentially expressed genes (DEGs) and dysregulated regulatory signaling pathways are associated with oncogenesis.7,8 Transcription factors (TFs), kinases, and miRNAs are key modulators at the transcriptional and post-transcriptional levels. Changes in the expression levels of TFs, kinases, and miRNAs might cause the oncogenesis and progression of NPC.9–11 These regulators do not work alone, and the cross-talk between TFs, kinases, and miRNA could link an interconnected molecular network to execute a particular role in the pathological process of NPC. Kinases and TFs are upstream regulators, and deregulated kinases and TFs might lead to the up- or down-regulation of downstream mRNAs.12 In addition, TFs have been reported to regulate the expression of mRNA by interacting with miRNAs.13 To this end, an integrated study involving mRNA, TFs, kinases, and miRNA is more helpful, as this could provide more comprehensive and credible insights into the molecular regulatory mechanisms in NPC. Recently, a variety of studies have reported the function of single kinases, TFs, miRNAs, and mRNAs in NPC. However, there have been few studies on integrated regulatory network analysis in NPC.

In the recent years, microarray technology has been widely applied in the examination of general genetic abnormalities.14 Based on microarray expression data, a large number of DEGs between NPC and normal nasopharyngeal tissue could be identified.15,16 However, as differences exist in the specimen sources, sample processing, and experimental batches, the reported DEGs vary from study to study. In addition, few studies have analyzed the relationship between DEGs and their transcriptional regulatory networks in NPC. In this study, to provide multiscale insight into NPC, we identified DEGs and differentially expressed miRNAs (DEMs) between NPC and normal nasopharyngeal tissue. Subsequently, we analyzed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and protein–protein interaction (PPI) network of the DEGs. We also predicted the key TFs, kinases, and miRNAs of the DEGs and then constructed a kinase–TF–mRNA–miRNA integrated network, and hub nodes were selected. Finally, candidate small-molecule agents were predicted based on DEGs using CMap, which might provide novel insights into the treatment of NPC.

Materials and Methods

Data Resources

The NPC mRNA and miRNA expression profiling dataset search was conducted using the NCBI Gene Expression Omnibus (GEO) database (). The datasets were searched using the following four criteria: (1) NPC, (2) tumor and normal control samples included, (3) sample resource, Homo sapiens, and (4) at least six samples included. Five mRNA microarray datasets (GSE53819, GSE12452, GSE64634, GSE40290, and GSE13597) and two miRNA microarray datasets (GSE22587 and GSE43039) were downloaded. GSE53819 included 18 NPC and 18 normal nasopharyngeal samples; GSE12452 included 31 NPC tissue samples and ten normal nasopharyngeal tissues; GSE64634 included 12 NPC tissue and four normal nasopharyngeal samples; GSE40290 contained 25 non-keratinizing NPC samples and eight normal nasopharyngeal tissue samples; GSE13597 included 25 NPC samples and three normal nasopharyngeal tissues; GSE22587 contained 8 NPC tissues compared to four normal nasopharyngeal tissues; GSE43039 included 20 poor or undifferentiated NPC samples and 20 normal nasopharyngeal specimens. As the difference in the sample number between NPC and control is large in GSE13597, this dataset was used as validation cohort in this study. Detailed information on these datasets is provided in Table 1. As all relevant datasets were publicly available through open access, ethics approval was not required in this study.

Table 1

Characteristics of mRNA and miRNA Expression Profiles of Nasopharyngeal Carcinoma (NPC)

GEO Series	Expression Type	Platform	Sample Number		Reference	Serve as
			Tumor	Normal
GSE53819	mRNA	GPL6480 Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version)	18	18	Bao et al.	Derivation
GSE12452	mRNA	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	31	10	Dodd et al.	Derivation
GSE64634	mRNA	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	12	4	Bo et al.	Derivation
GSE40290	mRNA	GPL8380 Capitalbio 22K Human oligo array version 1.0	25	8	Yao et al.	Derivation
GSE13597	mRNA	GPL96 [HG-U133A] Affymetrix Human Genome U133A Array	25	3	Bose et al.	Validation
GSE22587	miRNA	GPL8933 Illumina Human Beta-version microRNA expression BeadChip	8	4	LI et al	Derivation
GSE43039	miRNA	GPL16414 CCDTM-miRNA850-version 4p1.4	20	20	Lyu et al	Derivation

DEG and DEM Identification and Data Processing

Differential expression between NPC and normal control samples was analyzed using R software (version 3.5.3, ). Datasets GSE53819, GSE12452, GSE64634, and GSE40290 were used for screening the DEGs, and datasets GSE22587 and GSE43039 were used for DEM screening. DEG and DEM identification were conducted using the limma and GEOquery package in R. |log2FC(fold change)| > 1.0, and adjusted P-value (Adj. P) < 0.05, which served as the cut-off criteria for DEG and DEM screening in these microarray datasets. Volcano plots of DEGs and DEMs in each dataset were constructed using the ggplot2 package in R. Venn diagrams were produced using the Jvenn online tools ().17

GO and KEGG Pathway Enrichment Analyses of DEGs

GO is a database that provides gene annotation and gene function analysis in three subsets: biological process (BP), cellular component (CC), and molecular function (MF).18 The KEGG database contains gene functional pathway information.19 The GO and KEGG pathway enrichment analyses were performed using the WebGestalt tool (). WebGestalt is an online database that uses multiple algorithms for enrichment analysis involving a variety of functions.20 The Benjamini–Hochberg false-discovery rate (FDR) was adopted for screening statistically significant functional modules and pathways. Statistical significance was set at FDR < 0.05. The goplot and ggplot2 R packages were used to visualize GO and KEGG enrichment results, respectively.

Establishment of PPI Network and Topological Analysis

Interactions between DEGs were evaluated using the Search Tool for the Retrieval of Interacting Genes (STRING) database (Version 11.0; ).21 A combined score of ≥0.4 was set as the threshold. The constructed PPT network was visualized using Cytoscape software (Version 3.7.2).22 To screen the most significant modules in the PPI network, we used the molecular complex detection (MCODE) plugin in Cytoscape, which is based on the graph-theoretic clustering algorithm. The screening criteria were set as follows: degree cut-off = 2, maximum depth = 100, node score cut-off = 0.2, and k-Core = 2.

Prediction and Analysis of the Upstream TFs and Kinases of DEGs

In this study, potential TFs targeted to DEGs were inferred using the Transcriptional Regulatory Relationships Unraveled (TRRUST version 2, ) database.23 The Transcriptional Regulatory Element Database (TRED, ) was used for TF annotation.24 Additionally, Expression2Kinases (X2Kgui, Version 1.0) was applied to determine the association between the DEGs and upstream protein kinases.25 X2Kgui is software that uses a network-based reverse engineering approach to make predictions based on prior knowledge network data. To improve the prediction reliability, we analyzed the acquired results comprehensively in the above three databases and cross-referenced the TFs, which were targeted by kinases acquired in X2Kgui, with the TFs identified by TRRUST and TRED. A P-value < 0.05, was used to screen the predicted TFs and kinases. The interaction network of mRNA–TF and TF–kinase was visualized using Cytoscape.

Prediction and Analysis of the Potential miRNA of DEGs

The potential miRNAs of the DEGs were predicted using the miRWalk online database (Version 3.0) ().26,27 The cut-off value was a binding score > 0.95. The miRNAs predicted in the miRWalk were cross-referenced with the DEMs identified in the two miRNA microarray datasets (GSE22587 and GSE43039). Based on these overlapped DEMs, an mRNA–miRNA regulatory network was constructed using Cytoscape.

Construction of a Kinase–TF–mRNA–miRNA Integrated Network

According to the three networks acquired above, that is, the mRNA–TF network, TF–kinase network, and mRNA–miRNA network, an integrated kinase–TF–-mRNA–miRNA network was established to reveal the potential molecular mechanisms involved in the  carcinogenesis of NPC. Cytoscape was used to visualize and analyze the constructed network. The cytoHubba plugin in Cytoscape was used to identify the hub nodes (hub genes, hub TFs, hub kinases, and hub miRNAs) in the integrated network for further analysis.

Validation of Hub Gene Expression

To validate the expression levels of the hub genes selected in the kinase–TF–mRNA–miRNA network, another independent NPC mRNA expression profiling dataset, GSE13597, was used for analysis. The cut-off criteria were set as adjusted P-value <0.05, and |logFC| >1. We also validated the expression levels of the hub genes in the Oncomine database ().28 The threshold was P-value <1E−4, fold change > 2, and gene rank = top 10%. The results were screened by selecting Sengupta Head-Neck: nasopharyngeal carcinoma vs normal tissue.

Identification of Candidate Small Molecules Agents

The Connectivity Map database (CMap, ), a gene-expression-based drug development platform, was used to investigate potential small-molecule agents that might reverse the biological status of NPC.29 Based on the uploaded upregulated and downregulated DEGs, gene set enrichment analysis was performed to obtain small-molecule agents. Enrichment scores ranged from −1 to +1. The positive connectivity score (closer to +1) indicates that small-molecule agents could facilitate gene expression in NPC, and the negative connectivity value (close to −1) indicates a higher similarity between the gene and the small molecule, which represents the potential therapeutic value of this small molecule in NPC. We ranked the results by P-value, and the top six molecules were selected. PubChem (), a public repository of chemical information, including substance information, compound structures, and biological activities, was used to analyze the molecules obtained.

Results

Identification of DEGs and DEMs in NPC

A flowchart of the study design is presented in Figure 1. Four mRNA (GSE53819, GSE12452, GSE64634, and GSE40290) and two miRNA (GSE22587 and GSE43039) expression profile datasets, which contained 139 NPC samples and 67 normal control samples, were downloaded from the GEO database. For each dataset, differential expression analysis was performed between NPC and normal control samples based on the cut-off criteria |log2FC| > 1.0, and adjusted P-value < 0.05. DEGs and DEMs in each dataset were identified, and the volcano plots are shown in Figure 2A and C. Furthermore, co-DEGs and co-DEMs are presented by Venn diagrams (Figure 2B and D). We identified 122 DEGs, including 37 upregulated genes and 85 downregulated genes, and 44 DEMs between NPC and normal nasopharyngeal samples ().

Figure 1

A flowchart of the study design.

Figure 2

Identification of differentially expressed genes (DEGs) and miRNAs (DEMs). (A) The volcano plots for four mRNA expression datasets (GSE53819, GSE12452, GSE64634, and GSE40290); red dots indicate up-regulated DEGs; blue dots present down-regulated DEGs; and gray dots indicate nonsignificant DEGs. (B) Venn diagrams exhibit a total of 122 DEGs, consisting of 37 up-regulated and 85 down-regulated genes were identified in the four mRNA expression datasets. (C) The volcano plots for two miRNA expression datasets (GSE22587 and GSE43039). (D) Venn diagram shows 44 DEMs were identified in the two miRNA expression datasets.

GO functional and KEGG pathway enrichment analyses were performed using the WebGestalt database to further understand the biological processes and pathways involved in the DEGs. As shown in Figure 3, GO analysis indicated that DEGs were significantly enriched in the immune response, leukocyte activation, movement of cell or subcellular components, endoplasmic reticulum lumen, complex of collagen trimers, and extracellular matrix structural constituent (Figure 3A–C). KEGG analysis revealed that the DEGs were mainly enriched in the NF-κB signaling pathway, B cell receptor signaling pathway, hematopoietic cell lineage, small cell lung cancer, and microRNAs in cancer and other pathways (Figure 3D). The results indicated that DEGs were mostly involved in extra- and intracellular components, immune function, and cancer- and immune-related pathways.

Figure 3

GO and KEGG pathway enrichment analyses of differentially expressed genes (DEGs). (A) Circular dendrogram depicts the relationship between genes and GO terms of biological process. (B) Circular dendrogram depicts the relationship between genes and GO terms of cellular component. (C) Circular dendrogram depicts the relationship between genes and GO terms of molecular function. (D) Bubble plot exhibits the relationship between genes and KEGG pathways.

PPI Network and Module Analysis

PPIs among the DEGs were analyzed using the STRING online database. A total of 122 nodes and 198 edges were included in the PPI network (Figure 4A). The average node degree was 3.25, the PPI enrichment P-value was <1.0e−16, and the threshold of the minimum required interaction score was >0.4. Subsequently, we analyzed and found four significant modules from the PPI network using the MNC algorithm in the MCODE plugin (Figure 4B–E and ).

Figure 4

Protein–protein interaction (PPI) network of the differentially expressed genes (DEGs). (A) PPI network of the DEGs. The nodes represent proteins encoded by genes and the edges represent connections between the nodes. Red-colored nodes represent up-regulated genes, blue-colored nodes represent down-regulated genes. (B–E) Significant modules in the PPI network obtained by the MNC algorithm.

Construction of the Kinase–TF–mRNA–miRNA Interaction Network

To evaluate the function of DEGs in NPC, regulatory associations among DEGs, kinases, TFs, and miRNAs were analyzed. We predicted potential TFs targeted to DEGs using the TRRUST database and annotated them using TRED. In total, 26 TFs and 38 TF-targeted DEGs were identified (Figure 5A and Table 2). Subsequently, the upstream kinases of the DEGs were predicted using X2Kgui. To improve the prediction reliability, we analyzed the acquired results in the above three databases comprehensively and retrieved the kinases that also targeted TFs identified by TRRUST and TRED. Then, 73 kinases and 14 kinase-targeted TFs were identified, and an upstream kinase-downstream TF network was constructed (Figure 5B and ). The top 20 identified kinases are listed in Table 3.

Figure 5

Construction of the kinase-transcription factor (TF)–mRNA–miRNA interaction network. (A) The TF–mRNA network. (B) The kinase–TF network. (C) Venn diagram shows seven overlapping differentially expressed miRNA (DEMs) were identified between the miRWalk database and the DEMs identified in the two miRNA microarray datasets (GSE22587 and GSE43039). (D) The mRNA–miRNA network. (E) The kinase–TF–mRNA–miRNA interaction network. Red-colored circles: up-regulated genes; blue-colored circles: down-regulated genes; green rectangles: TFs; purple hexagon represent kinase; kelly rhombus: miRNA.

Table 2

The Predicted Transcription Factors (TFs) and the TF-Targeted Differentially Expressed Genes (DEGs)

Key TF	P value	Q value	Gene Counts	DEGs
NFKB1	2.39E-05	0.00066	10	MMP3,PTGDS,PIGR,MMP1,PLAU,FN1,BMP2,PTGS2,LCN2,CR2
FOS	3.07E-05	0.00066	5	CLU,MMP3,MMP1,PTGS2,PLAU
RELA	0.000131	0.00172	9	PTGS2,BMP2,PLAU,PIGR,MMP1,CR2,MMP3,FN1,LCN2
CEBPD	0.00016	0.00172	3	C3,TNFAIP6,PTGS2
CEBPB	0.000576	0.00351	4	MMP1,CR2,PTGS2,TNFAIP6
SPI1	0.000652	0.00351	4	IFIT3,MS4A1,CD22,ALOX15
ATF2	0.00109	0.00435	3	FN1,PTGS2,PLAU
JUND	0.0013	0.00435	3	SOX4,PTGS2,CLU
ATF4	0.00141	0.00435	3	PTGS2,S100P,PLAU
STAT1	0.00203	0.00571	4	IFIT3,CD22,LCN2,PTGS2
STAT2	0.00212	0.00571	2	PTGS2,IFIT3
EGR1	0.0024	0.00608	4	FN1,LTB,PLAU,PTGS2
JUN	0.00264	0.0063	5	MMP3,GJA1,PLAU,PTGS2,MMP1
TP53	0.00398	0.00858	5	GPNMB,KLF2,MMP1,PRC1,SSTR2
JUNB	0.00399	0.00858	2	PTGS2,CLU
PAX5	0.00778	0.0152	2	MMP1,CD19
USF1	0.00822	0.0153	3	PTGS2,PTGDS,PIGR
ETV4	0.00853	0.0153	2	MMP1,PTGS2
SP1	0.0106	0.0168	8	GJA1,ADAM23,PLAU,CD1D,PTGS2,CXCL14,LTF,PROM1
ETS1	0.014	0.0207	3	MMP1,PLAU,MMP3
ETS2	0.0176	0.0244	2	MMP1,MMP3
CREB1	0.0197	0.0265	3	PTGS2,MSMB,PLAU
AR	0.0215	0.027	3	FN1,PIGR,PTGS2
STAT6	0.0219	0.027	2	ALOX15,PTGS2
POU2F1	0.0243	0.029	2	PLAU,POU2AF1
USF2	0.036	0.0418	2	PIGR,PTGS2

Table 3

Top 20 List of the Identified Kinases and Kinase-Targeted Transcription Factors (TFs)

Kinase	P-value	Z-Score	TF Count	TF
MAPK1	5.22E-18	−2.2592031	8	FOS,RELA,SP1,CEBPB,JUN,JUND,AR,TP53
MAPK3	2.72E-12	−2.1817021	7	TP53,FOS,SP1,CEBPB,JUND,JUN,AR
GSK3B	4.18E-12	−2.1347149	9	JUN,JUND,CEBPB,NFKB1,JUNB,AR,TP53,RELA,CREB1
MAPK14	1.58E-11	−2.2678044	9	FOS,TP53,JUN,RELA,JUNB,CREB1,AR,USF1,JUND
CDK1	7.06E-11	−2.0346413	5	CEBPB,TP53,POU2F1,AR,JUND
CSNK2A1	7.85E-10	−1.9961411	7	FOS,NFKB1,SP1,TP53,RELA,JUN,CREB1
CDK2	8.35E-10	−1.9852625	5	CEBPB,FOS,TP53,POU2F1,JUN
AKT1	8.74E-10	−2.0687625	2	AR,CREB1
MAPK8	3.25E-09	−1.8922677	7	AR,TP53,SP1,JUNB,JUN,CEBPB,JUND
CHUK	4.73E-09	−1.6282563	3	NFKB1,RELA,TP53
IKBKB	9.02E-09	−1.5826631	3	NFKB1,RELA,FOS
PRKDC	7.63E-07	−1.5704482	7	POU2F1,JUN,NFKB1,SP1,RELA,TP53,USF1
MAP3K7	1.45E-06	−1.4165878	2	RELA,FOS
ALK	1.58E-06	−0.9642052	2	TP53,RELA
CDK7	2.29E-06	−0.6642677	3	AR,TP53,POU2F1
ABL1	2.89E-06	−1.3842712	4	JUN,TP53,CREB1,CEBPB
IKBKE	2.89E-06	−1.1826175	2	RELA,NFKB1
MAP3K8	2.92E-06	−0.9652278	3	RELA,NFKB1,CREB1
SRC	5.82E-06	−1.4174873	1	AR
PRKCD	1.02E-05	−1.2348518	1	TP53

In further analysis, we used the miRWalk online database to predict potential miRNAs for DEGs, and 2499 candidate miRNAs were retrieved (). The miRNAs predicted in miRWalk were cross-referenced with the DEMs identified in the two miRNA microarray datasets (GSE22587 and GSE43039). Seven overlapping miRNAs (hsa-miR-184, hsa-miR-30a-3p, hsa-miR-346, has-miR-380-3p, has-miR-448, hsa-miR-518b, and hsa-miR-521) were identified (Figure 5C and D, Table 4).

Table 4

The List of Overlapping miRNAs

miRNA	mRNA Count	mRNA
hsa-miR-184	8	RAB37,XG,CHST9,DNAI1,GALNT11,CLU,ADAM23,SPA17
hsa-miR-30a-3p	2	CHST9,CD1D
hsa-miR-346	13	GJA1,RAB37,RIBC1,RRAD,PROM1,EFHC1,CHST9,BCAS1,COL8A1,TNFSF11,CH25H,CAPS
hsa-miR-380-3p	1	LRRC6
hsa-miR-448	5	MEF2C,C22orf23,MUC20,CH25H,HOXA10
hsa-miR-518b	4	PROM1,RALGPS2,CDH26,GCNT3
hsa-miR-521	3	DST,PRC1

Finally, we built a kinase–TF–mRNA–miRNA interaction network in Cytoscape (Figure 5E). This network revealed hub nodes and their interactions associated with the molecular mechanisms of NPC. The properties of each node within the integrated network were analyzed using the cytoHubba plugin. Based on the degree of node algorithm, we selected the hub nodes using the top eight hub genes (PTGS2, FN1, MMP1, PLAU, MMP3, CD19, BMP2, and PIGR), the top six TFs (TP53, RELA, CREB1, JUN, FOS, and AR), the top three kinases (MAPK14, GSK3B, and MAPK1), and the top two miRNAs (hsa-miR-184 and hsa-miR-346). These hub nodes may play a crucial role in the pathological process of NPC.

To confirm the reliability of our results, we analyzed another independent NPC mRNA expression profiling dataset GSE13597 and NPC mRNA expression data from Oncomine. We found that all hub genes were expressed in the GSE13597 dataset and Oncomine. In the GSE13597 dataset, the expression of PTGS2, FN1, MMP1, PLAU, MMP3, and BMP2 were significantly upregulated, whereas CD19 and PIGR expression were significantly downregulated in NPC samples. These results were consistent with the analyzed DEGs (). Moreover, we obtained the same outcomes by validating the results in the Oncomine database, as shown in Figure 6A–H. The validation results indicate the stability and reliability of our analysis.

Figure 6

Validation of the expression of hub genes in the Oncomine database. (A) PTGS2; (B) FN1; (C) MMP1; (D) PLAU; (E) MMP3; (F) CD19; (G) BMP2; (H) PIGR. Under these box plots, 0 indicates normal; 1 indicates nasopharyngeal carcinoma (NPC).

Prediction of Candidate Small-Molecule Agents in NPC

CMap was used to identify candidate small molecules that have potential for NPC treatment. Based on gene set enrichment analysis, 98 molecules were obtained from a total of 6100 results (). We ranked the results by P-values. The top six molecules, phenoxybenzamine (enrichment score =  −0.917), luteolin (enrichment score =  −0.902), thioguanosine (enrichment score =  −0.869), reserpine (enrichment score =  −0.863), blebbistatin (enrichment score =  −0.857), and camptothecin (enrichment score =  −0.847), exhibited highly significant negative scores and showed potential for reversing NPC tumor status (Table 5). Figure 7 shows the 3D structures of these six molecules.

Table 5

Prediction of Candidate Small Molecule Agents in Nasopharyngeal Carcinoma (NPC)

CMap Name	Mean	Enrichment	p-value	Specificity
Phenoxybenzamine	−0.743	−0.917	0.0001	0.0091
Luteolin	−0.752	−0.902	0.00016	0
Thioguanosine	−0.707	−0.869	0.00058	0.007
Reserpine	−0.71	−0.863	0.00519	0.0201
Blebbistatin	−0.643	−0.857	0.04072	0.0608
Camptothecin	−0.707	−0.847	0.00709	0.1818

Figure 7

Three-dimensional structure of six candidate small molecule agents in nasopharyngeal carcinoma (NPC). (A) phenoxybenzamine; (B) luteolin; (C) thioguanosine; (D) reserpine; (E) blebbistatin; and (F) camptothecin.

Discussion

Because of the lack of early characteristic symptoms and the occult occurrence site, most NPC patients (> 70%) are diagnosed at a locoregionally advanced stage.30 Therefore, the screening of biomarkers for early diagnosis and in-depth study of the regulatory pathways involved in the oncogenesis and tumor progression of NPC are crucial research aims. Meanwhile, the rapid development of recent bioinformatic technologies and high-throughput technologies, for example, microarrays and RNA sequencing, have facilitated the exploration of mechanisms in numerous diseases, including NPC.31 In the present study, we identified 122 DEGs and 44 DEMs between NPC and normal samples by comprehensive bioinformatics analysis. We then constructed an integrated regulatory network involving kinases, TFs, mRNAs, and miRNAs in NPC. Eight genes (PTGS2, FN1, MMP1, PLAU, MMP3, CD19, BMP2, and PIGR) were identified as hub genes. We also predicted six candidate small-molecule agents (phenoxybenzamine, luteolin, thioguanosine, reserpine, blebbistatin, and camptothecin) with potential therapeutic capacity against NPC.

In this study, enriched GO terms and KEGG pathways demonstrated numerous differences between the NPC and normal tissues. GO terms “immune response,” “leukocyte activation,” and “endoplasmic reticulum lumen” were enriched in NPC samples. Previous studies have clarified that changes in the immune response, leukocyte activation, and endoplasmic reticulum stress may affect the oncogenesis, metastasis, and apoptosis of NPC.30,32,33 Moreover, Hu et al indicated that IL-18 could stimulate T cells and NK cells to produce IFN-γ, which consequently activates macrophages and other immune cells to secrete chemokines that initiate a leukocyte recruitment cascade.34 Furthermore, the KEGG pathways NF-κB signaling pathway and B cell receptor signaling pathway were enriched and regulated by DEGs. The NF-κB signaling pathway has been shown to promote cell proliferation, migration, and angiogenesis in NPC, and this pathway also exhibits a predictive effect in the clinical response to radiotherapy for NPC patients.35–37 EBV infection is thought to be a potential risk factor for the progression of NPC. EBV infection has been reported to affect the function of the B cell receptor signaling pathway in NPC.38,39 These findings confirmed that the oncogenesis and progression of NPC is closely related to defects in the integrity of the cellular structure and function, and the dysregulated pathways.

Eight genes, including PTGS2, FN1, MMP1, PLAU, MMP3, CD19, BMP2, and PIGR, were selected as hub genes in this study. PTGS2 (cyclooxygenase-2, COX-2), is a crucial enzyme in the conversion course of arachidonic acid to prostaglandins.40 It has been reported that PTGS2 promotes metastasis in NPC by regulating the interaction between cancer cells and myeloid-derived suppressor cells.41 Tan et al showed that PTGS2 expression increases as nasopharyngeal epithelium progresses from normal to dysplastic to carcinoma by immunohistochemistry and semiquantitative assessment.42 Zhong et al demonstrated that the lncRNA GAS5/miR-4465 axis regulated the malignant potential of NPC by targeting PTGS2.43 These findings suggest that PTGS2 might be associated with the metastasis and progression of NPC. FN1 (Fibronectin 1) is a member of the glycoprotein family located on chromosome 2q35. The up-regulation of FN1 is related to cell adhesion and migration in many malignancies.44 The overexpression of FN1 in NPC is associated with an advanced disease stage and short survival.45 Wang et al indicated that FN1 suppressed apoptosis via the NF-κB pathway, which might be related to migration in NPC.46 MMP1 and MMP3 belong to the matrix metalloproteinase family.47,48 MMP1 is reported to promote NPC cell proliferation, suppress apoptosis, and decrease the sensitivity of NPC cells to 5-FU.47 MMP3 is a lytic transactivator of EBV, which promotes the progression and invasion of NPC.48 BMP2 (Bone Morphogenetic Protein 2) is a secreted protein that promotes proliferation and invasion of NPC.49 Unlike the aforementioned upregulated DEGs, CD19 and PIGR (Polymeric Immunoglobulin Receptor) were downregulated in NPC tissue compared with normal tissue. Xu et al reported that a low level of circulating CD19+ B cells meant poorer prognosis in EBV-associated NPC patients.50 Qi et al indicated that serum PIGR expression could serve as an independent prognostic predictor for NPC.51 PLAU is reported to be associated with tumor cell migration, invasion, and metastasis in several cancers.52 However, the role of PLAU in NPC has not yet been studied, and we will focus on this gene in our future study.

Six TFs (TP53, RELA, CREB1, JUN, FOS, and AR) were identified as hub TFs based on integrated network analysis. All these TFs have been verified as key modulators and potential therapeutic targets in a wide variety of human malignancies.53–57 TP53 (tumor protein p53) is a classical tumor suppressor. High frequencies of TP53 mutations and/or deletions have been found in NPC tissues.53 RELA (v-rel reticuloendotheliosis viral oncogene homolog A, p65) is an NF-κB family member and a subunit of the NF-κB TF complex. Activation of the NF-κB signaling pathway could significantly affect the development and progression of NPC.54 In the present study, our KEGG analysis also showed the enriched NF-κB signaling pathway in NPC samples. CREB1 (cAMP-responsive element-binding protein 1) is a nuclear TF that belongs to the basic leucine zipper family. Although the role of CREB1 in NPC has not been reported, CREB1 has been proven to promote tumor proliferation, invasion, and metastasis in several malignancies.55 JUN (jun proto-oncogene) and FOS (FBJ murine osteosarcoma viral oncogene homolog) belong to the AP1 family. The regulatory effects of JUN and FOS in tumorigenesis and progression have been reported in numerous cancers, including NPC.56,58 AR (androgen receptor) is a ligand-dependent TF that belongs to the nuclear receptor superfamily. A recent study indicated that AR-induced long non-coding RNA LINC01503 promotes NPC proliferation and metastasis through the SFPQ-FOSL1 axis.57

In addition, three kinases (MAPK14, GSK3B, and MAPK1) and two miRNAs (hsa-miR-184 and hsa-miR-346) were identified as hub modulators in the regulatory network. MAPK14 (P38α) and MAPK1 (ERK2) are important members of the MAPK/ERK pathway and are associated with tumor cell proliferation, differentiation, and apoptosis in numerous cancers, including NPC. Activated MAPK/ERK was found in the NPC samples, while the blockage of this pathway induced NPC cells apoptosis.59,60 MAPK14 and MAPK1 have been reported to promote the carcinogenesis of NPC by activating p53.61 GSK3B (glycogen synthase kinase 3 beta) is a widely expressed serine/threonine kinase and play a crucial regulatory role in the metabolism, oncogenesis, and neurological diseases.62 Hu et al reported that the AKT/GSK3B/beta-catenin pathway is associated with cell proliferation, invasion, and metastasis in NPC.63 A recent study indicated that as a tumor-suppressive miRNA, hsa-miR-184 could inhibit the tumor invasion, migration, and metastasis in NPC by targeting Notch2.64 Hsa-miR-346 also could promote the migration and invasion of NPC.65 These findings verified the reliability of our prediction and provide several new experimental directions for future NPC research.

As a radiosensitive carcinoma, radiotherapy is the primary treatment for NPC.3 While NPC is also chemosensitive, adjuvant, or concurrent chemotherapy along with radiotherapy could improve the prognosis of NPC.66 In this study, we explored promising small-molecule agents for NPC using the CMAP database and screened 98 small-molecule agents. Among the top 6 small-molecule agents, luteolin is a plant flavonoid that owns multiple biological properties such as anti-oxidant, anti-allergy, anti-inflammation, and anti-cancer.67 Luteolin has been verified to induce G1 arrest and inhibit the reactivation of EBV in NBC cells.68,69 Camptothecin is a cytotoxic alkaloid that has been proven to attenuate the replication of cancer cells by inhibiting DNA topoisomerase 1 in human malignancies.70 Li et al reported that camptothecin could inhibit the progression of NPC by activation of the PI3K/AKT signaling pathway.71 However, further in vivo and in vitro studies, as well as clinical trials, are needed to validate the safety and therapeutic efficacy of these small-molecule agents on NPC.

In the present study, we identified candidate mRNAs, TFs, kinases, and miRNAs involved in the oncogenesis and tumor progression of NPC through a series of comprehensive bioinformatics analyses. To the best of our knowledge, our study is the first to integrate a regulatory network involving kinases, TFs, mRNA, and miRNA in NPC. Moreover, to improve the prediction reliability, in the present study, we adopted vigor selection criteria to select the datasets for DEG and DEM screening. The candidate DEGs, TFs, kinases, and miRNA were identified by cross-referencing the data from at least two datasets or databases, and the hub genes were validated using two databases, and consistent results were obtained. However, our study had some limitations. First, as an integrated bioinformatics analysis, the heterogeneity among different datasets, data platforms, and statistical analysis is almost inevitable, which might impact the reliability of this study. Second, although the algorithms and R packages used in this study are the latest and most reasonable, random errors and selection bias are hardly avoided. Third, the predictive findings of this study have not been validated by in vitro or in vivo experiments. Further experimental investigations are needed to corroborate these novel targets and explore the potential underlying mechanisms involved.

Conclusion

In summary, we conducted a comprehensive bioinformatics analysis of DEGs between NPC and normal tissues. We also integrated a regulatory network involving kinases, TFs, mRNAs, and miRNAs in NPC. Eight hub genes (PTGS2, FN1, MMP1, PLAU, MMP3, CD19, BMP2, and PIGR) were screened in the integrated network. Six candidate small-molecule agents (phenoxybenzamine, luteolin, thioguanosine, reserpine, blebbistatin, and camptothecin) were predicted. These findings might provide information for a better understanding of the oncogenesis of NPC and reveal promising insights into the pathogenesis, treatment, and prognosis of NPC. Further experimental validation will aid in the confirmation of our findings and molecular understanding of the underlying mechanisms.