1Department of Otolaryngology, Head and Neck Surgery, Graduate School of Medicine, Kyoto University, Kyoto, Japan; 2Department of Otolaryngology-Head and Neck Surgery, Shanghai Zhongshan Hospital Affiliated to Fudan University, Shanghai, People’s Republic of China
Correspondence: Xinsheng Huang
Department of Otolaryngology-Head and Neck Surgery, Shanghai Zhongshan Hospital Affiliated to Fudan University, Shanghai, 200032, People’s Republic of China
Email [email protected]
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.
Keywords: nasopharyngeal carcinoma, differentially expressed genes, microRNA, transcription factor, kinase, bioinformatics
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
The NPC mRNA and miRNA expression profiling dataset search was conducted using the NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). 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)
DEG and DEM Identification and Data Processing
Differential expression between NPC and normal control samples was analyzed using R software (version 3.5.3, https://www.r-project.org/). 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 (http://jvenn.toulouse.inra.fr/app/index.html).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 (https://www.webgestalt.org). 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; http://www.string-db.org/).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, http://www.grnpedia.org/trrust/) database.23 The Transcriptional Regulatory Element Database (TRED, http://rulai.cshl.edu/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) (http://mirwalk.umm.uni-heidelberg.de).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 (http://www.oncomine.org).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, http://www.broadinstitute.org/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 (https://pubchem.ncbi.nlm.nih.gov/), a public repository of chemical information, including substance information, compound structures, and biological activities, was used to analyze the molecules obtained.
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 (Supplementary Table 1).
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 and KEGG Pathway Enrichment Analyses of DEGs
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 Supplementary Table 2).
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 Supplemental Table 3). The top 20 identified kinases are listed in Table 3.
Table 2 The Predicted Transcription Factors (TFs) and the TF-Targeted Differentially Expressed Genes (DEGs)
Table 3 Top 20 List of the Identified Kinases and Kinase-Targeted Transcription Factors (TFs)
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.
In further analysis, we used the miRWalk online database to predict potential miRNAs for DEGs, and 2499 candidate miRNAs were retrieved (Supplementary Table 4). 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
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.
Validation of Hub Gene Expression
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 (Supplementary Table 5). 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 (Supplementary Table 6). 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)
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.
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.
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.
We thank International Science Editing for editing this manuscript. We would also like to thank Wenhao Zhang for the technical assistance in this research.
This research was funded by the National Natural Science Foundation of China (82000980) and Shanghai Shen Kang Hospital Development Center Clinical Research Plan (SHDC12018118).
The author reports no conflicts of interest in this work.
1. Wei WI, Sham JST. Nasopharyngeal carcinoma. Lancet (London, England). 2005;365(9476):2041–2054. doi:10.1016/S0140-6736(05)66698-6
2. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424. doi:10.3322/caac.21492
3. Chen Y-P, Chan ATC, Le Q-T, Blanchard P, Sun Y, Ma J. Nasopharyngeal carcinoma. Lancet (London, England). 2019;394(10192):64–80. doi:10.1016/S0140-6736(19)30956-0
4. Tian Y, Tang L, Yi P, et al. MiRNAs in radiotherapy resistance of nasopharyngeal carcinoma. J Cancer. 2020;11(13):3976–3985. doi:10.7150/jca.42734
5. Chou J, Lin Y-C, Kim J, et al. Nasopharyngeal carcinoma–review of the molecular mechanisms of tumorigenesis. Head Neck. 2008;30(7):946–963. doi:10.1002/hed.20833
6. Liu J-P, Cassar L, Pinto A, Li H. Mechanisms of cell immortalization mediated by EB viral activation of telomerase in nasopharyngeal carcinoma. Cell Res. 2006;16(10):809–817. doi:10.1038/sj.cr.7310098
7. Zeng Z, Huang H, Huang L, et al. Regulation network and expression profiles of Epstein-Barr virus-encoded microRNAs and their potential target host genes in nasopharyngeal carcinomas. Sci China Life Sci. 2014;57(3):315–326. doi:10.1007/s11427-013-4577-y
8. Carroll PA, Freie BW, Mathsyaraja H, Eisenman RN. The MYC transcription factor network: balancing metabolism, proliferation and oncogenesis. Front Med. 2018;12(4):412–425. doi:10.1007/s11684-018-0650-z
9. Lahiry P, Torkamani A, Schork NJ, Hegele RA. Kinase mutations in human disease: interpreting genotype-phenotype relationships. Nat Rev Genet. 2010;11(1):60–74. doi:10.1038/nrg2707
10. Chen H-C, Chen G-H, Chen Y-H, et al. MicroRNA deregulation and pathway alterations in nasopharyngeal carcinoma. Br J Cancer. 2009;100(6):1002–1011. doi:10.1038/sj.bjc.6604948
11. Razak ARA, Siu LL, Liu F-F, Ito E, O’Sullivan B, Chan K. Nasopharyngeal carcinoma: the next challenges. Eur J Cancer. 2010;46(11):1967–1978. doi:10.1016/j.ejca.2010.04.004
12. Mastrogamvraki N, Zaravinos A. Signatures of co-deregulated genes and their transcriptional regulators in colorectal cancer. NPJ Syst Biol Appl. 2020;6(1):23. doi:10.1038/s41540-020-00144-8
13. Martinez NJ, Walhout AJM. The interplay between transcription factors and microRNAs in genome-scale regulatory networks. Bioessays. 2009;31(4):435–445. doi:10.1002/bies.200800212
14. Lou Y, Jiang H, Cui Z, Wang X, Wang L, Han Y. Gene microarray analysis of lncRNA and mRNA expression profiles in patients with high-grade ovarian serous cancer. Int J Mol Med. 2018;42(1):91–104. doi:10.3892/ijmm.2018.3588
15. Jiang X, Feng L, Dai B, Li L, Lu W. Identification of key genes involved in nasopharyngeal carcinoma. Braz J Otorhinolaryngol. 2017;83(6):670–676. doi:10.1016/j.bjorl.2016.09.003
16. Zhang J-Z, Wu Z-H, Cheng Q. Screening and identification of key biomarkers in nasopharyngeal carcinoma: evidence from bioinformatic analysis. Medicine (Baltimore). 2019;98(48):e17997. doi:10.1097/MD.0000000000017997
17. Bardou P, Mariette J, Escudié F, Djemiel C, Klopp C. jvenn: an interactive Venn diagram viewer. BMC Bioinform. 2014;15(1):293. doi:10.1186/1471-2105-15-293
18. Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–29. doi:10.1038/75556
19. Altermann E, Klaenhammer TR. PathwayVoyager: pathway mapping using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. BMC Genomics. 2005;6:60. doi:10.1186/1471-2164-6-60
20. Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019;47(W1):W199–W205. doi:10.1093/nar/gkz401
21. Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–D613. doi:10.1093/nar/gky1131
22. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–2504. doi:10.1101/gr.1239303
23. Han H, Cho J-W, Lee S, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res. 2018;46(D1):D380–D386. doi:10.1093/nar/gkx1013
24. Jiang C, Xuan Z, Zhao F, Zhang MQ. TRED: a transcriptional regulatory element database, new entries and other development. Nucleic Acids Res. 2007;35(Database issue):D137–40. doi:10.1093/nar/gkl1041
25. Chen EY, Xu H, Gordonov S, Lim MP, Perkins MH, Ma’ayan A. Expression2Kinases: mRNA profiling linked to multiple upstream regulatory layers. Bioinformatics. 2012;28(1):105–111. doi:10.1093/bioinformatics/btr625
26. Dweep H, Sticht C, Pandey P, Gretz N. miRWalk–database: prediction of possible miRNA binding sites by “walking” the genes of three genomes. J Biomed Inform. 2011;44(5):839–847. doi:10.1016/j.jbi.2011.05.002
27. Dweep H, Gretz N, Sticht C. miRWalk database for miRNA-target interactions. Methods Mol Biol. 2014;1182:289–305. doi:10.1007/978-1-4939-1062-5_25
28. Rhodes DR, Yu J, Shanker K, et al. ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia. 2004;6(1):1–6. doi:10.1016/s1476-5586(04)80047-2
29. Lamb J, Crawford ED, Peck D, et al. The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science. 2006;313(5795):1929–1935. doi:10.1126/science.1132939
30. Chua MLK, Wee JTS, Hui EP, Chan ATC. Nasopharyngeal carcinoma. Lancet (London, England). 2016;387(10022):1012–1024. doi:10.1016/S0140-6736(15)00055-0
31. Cho WC-S. Nasopharyngeal carcinoma: molecular biomarker discovery and progress. Mol Cancer. 2007;6:1. doi:10.1186/1476-4598-6-1
32. Li X-J, Peng L-X, Shao J-Y, et al. As an independent unfavorable prognostic factor, IL-8 promotes metastasis of nasopharyngeal carcinoma through induction of epithelial-mesenchymal transition and activation of AKT signaling. Carcinogenesis. 2012;33(7):1302–1309. doi:10.1093/carcin/bgs181
33. Song L, Liu H, Ma L, Zhang X, Jiang Z, Jiang C. Inhibition of autophagy by 3-MA enhances endoplasmic reticulum stress-induced apoptosis in human nasopharyngeal carcinoma cells. Oncol Lett. 2013;6(4):1031–1038. doi:10.3892/ol.2013.1498
34. Hu H, Tang KF, Chua YN, et al. Expression of interleukin-18 by nasopharyngeal carcinoma cells: a factor that possibly initiates the massive leukocyte infiltration. Hum Pathol. 2004;35(6):722–728. doi:10.1016/j.humpath.2004.01.026
35. Peng T, Hu M, Wu -T-T, et al. Andrographolide suppresses proliferation of nasopharyngeal carcinoma cells via attenuating NF-κB pathway. Biomed Res Int. 2015;2015:735056. doi:10.1155/2015/735056
36. Guo C, Huang Y, Yu J, et al. The impacts of single nucleotide polymorphisms in genes of cell cycle and NF-kB pathways on the efficacy and acute toxicities of radiotherapy in patients with nasopharyngeal carcinoma. Oncotarget. 2017;8(15):25334–25344. doi:10.18632/oncotarget.15835
37. Zhao W, Ma N, Wang S, et al. RERG suppresses cell proliferation, migration and angiogenesis through ERK/NF-κB signaling pathway in nasopharyngeal carcinoma. J Exp Clin Cancer Res. 2017;36(1):88. doi:10.1186/s13046-017-0554-9
38. Caldwell RG, Wilson JB, Anderson SJ, Longnecker R. Epstein-Barr virus LMP2A drives B cell development and survival in the absence of normal B cell receptor signals. Immunity. 1998;9(3):405–411. doi:10.1016/s1074-7613(00)80623-8
39. Young LS, Dawson CW. Epstein-Barr virus and nasopharyngeal carcinoma. Chin J Cancer. 2014;33(12):581–590. doi:10.5732/cjc.014.10197
40. Trifan OC, Hla T. Cyclooxygenase-2 modulates cellular growth and promotes tumorigenesis. J Cell Mol Med. 2003;7(3):207–222. doi:10.1111/j.1582-4934.2003.tb00222.x
41. Li Z-L, Ye S-B, OuYang L-Y, et al. COX-2 promotes metastasis in nasopharyngeal carcinoma by mediating interactions between cancer cells and myeloid-derived suppressor cells. Oncoimmunology. 2015;4(11):e1044712. doi:10.1080/2162402X.2015.1044712
42. Tan K-B, Putti TC. Cyclooxygenase 2 expression in nasopharyngeal carcinoma: immunohistochemical findings and potential implications. J Clin Pathol. 2005;58(5):535–538. doi:10.1136/jcp.2004.021923
43. Zhong Q, Wang Z, Liao X, Wu R, Guo X. LncRNA GAS5/miR-4465 axis regulates the malignant potential of nasopharyngeal carcinoma by targeting COX2. Cell Cycle. 2020;19(22):3004–3017. doi:10.1080/15384101.2020.1816280
44. Wilson CB, Leopard J, Cheresh DA, Nakamura RM. Extracellular matrix and integrin composition of the normal bladder wall. World J Urol. 1996;14(Suppl 1):S30–7. doi:10.1007/BF00182062
45. Ma L-J, Lee S-W, Lin L-C, et al. Fibronectin overexpression is associated with latent membrane protein 1 expression and has independent prognostic value for nasopharyngeal carcinoma. Tumour Biol. 2014;35(2):1703–1712. doi:10.1007/s13277-013-1235-8
46. Wang J, Deng L, Huang J, et al. High expression of Fibronectin 1 suppresses apoptosis through the NF-κB pathway and is associated with migration in nasopharyngeal carcinoma. Am J Transl Res. 2017;9(10):4502–4511.
47. Song L, Liu H, Liu Q. Matrix metalloproteinase 1 promotes tumorigenesis and inhibits the sensitivity to 5-fluorouracil of nasopharyngeal carcinoma. Biomed Pharmacother. 2019;118:109120. doi:10.1016/j.biopha.2019.109120
48. Lee DCW, Chua DTT, Wei WI, Sham JST, Lau ASY. Induction of matrix metalloproteinases by Epstein-Barr virus latent membrane protein 1 isolated from nasopharyngeal carcinoma. Biomed Pharmacother. 2007;61(9):520–526. doi:10.1016/j.biopha.2007.08.007
49. Wang M-H, Zhou X-M, Zhang M-Y, et al. BMP2 promotes proliferation and invasion of nasopharyngeal carcinoma cells via mTORC1 pathway. Aging (Albany NY). 2017;9(4):1326–1340. doi:10.18632/aging.101230
50. Xu T, Huang Z, Su B, et al. Prognostic significance of circulating CD19+ B lymphocytes in EBV-associated nasopharyngeal carcinoma. Med Oncol. 2014;31(10):198. doi:10.1007/s12032-014-0198-y
51. Qi X, Li X, Sun X. Reduced expression of polymeric immunoglobulin receptor (pIgR) in nasopharyngeal carcinoma and its correlation with prognosis. Tumour Biol. 2016;37(8):11099–11104. doi:10.1007/s13277-016-4791-x
52. Ai C, Zhang J, Lian S, et al. FOXM1 functions collaboratively with PLAU to promote gastric cancer progression. J Cancer. 2020;11(4):788–794. doi:10.7150/jca.37323
53. Xiao M, Zhang L, Zhu X, et al. Genetic polymorphisms of MDM2 and TP53 genes are associated with risk of nasopharyngeal carcinoma in a Chinese population. BMC Cancer. 2010;10:147. doi:10.1186/1471-2407-10-147
54. Chung GT-Y, Lou WP-K, Chow C, et al. Constitutive activation of distinct NF-κB signals in EBV-associated nasopharyngeal carcinoma. J Pathol. 2013;231(3):311–322. doi:10.1002/path.4239
55. Rao M, Zhu Y, Cong X, Li Q. Knockdown of CREB1 inhibits tumor growth of human gastric cancer in vitro and in vivo. Oncol Rep. 2017;37(6):3361–3368. doi:10.3892/or.2017.5636
56. Milde-Langosch K. The Fos family of transcription factors and their role in tumourigenesis. Eur J Cancer. 2005;41(16):2449–2461. doi:10.1016/j.ejca.2005.08.008
57. He S-W, Xu C, Li Y-Q, et al. AR-induced long non-coding RNA LINC01503 facilitates proliferation and metastasis via the SFPQ-FOSL1 axis in nasopharyngeal carcinoma. Oncogene. 2020;39(34):5616–5632. doi:10.1038/s41388-020-01388-8
58. Tulalamba W, Janvilisri T. Nasopharyngeal carcinoma signaling pathway: an update on molecular biomarkers. Int J Cell Biol. 2012;2012:594681. doi:10.1155/2012/594681
59. Hirota Y, Yamashita S, Kurihara Y, et al. Mitophagy is primarily due to alternative autophagy and requires the MAPK1 and MAPK14 signaling pathways. Autophagy. 2015;11(2):332–343. doi:10.1080/15548627.2015.1023047
60. Bonney EA. Mapping out p38MAPK. Am J Reprod Immunol. 2017;77:5. doi:10.1111/aji.12652
61. Li L, Guo L, Tao Y, et al. Latent membrane protein 1 of Epstein-Barr virus regulates p53 phosphorylation through MAP kinases. Cancer Lett. 2007;255(2):219–231. doi:10.1016/j.canlet.2007.04.014
62. Kockeritz L, Doble B, Patel S, Woodgett JR. Glycogen synthase kinase-3–an overview of an over-achieving protein kinase. Curr Drug Targets. 2006;7(11):1377–1388. doi:10.2174/1389450110607011377
63. Hu W, Xiao L, Cao C, Hua S, Wu D. UBE2T promotes nasopharyngeal carcinoma cell proliferation, invasion, and metastasis by activating the AKT/GSK3β/β-catenin pathway. Oncotarget. 2016;7(12):15161–15172. doi:10.18632/oncotarget.7805
64. Zhu H-M, Jiang X-S, Li H-Z, et al. miR-184 inhibits tumor invasion, migration and metastasis in nasopharyngeal carcinoma by targeting Notch2. Cell Physiol Biochem. 2018;49(4):1564–1576. doi:10.1159/000493459
65. Yan H-L, Li L, Li S-J, Zhang H-S, Xu W. miR-346 promotes migration and invasion of nasopharyngeal carcinoma cells via targeting BRMS1. J Biochem Mol Toxicol. 2016;30(12):602–607. doi:10.1002/jbt.21827
66. Blanchard P, Lee A, Marguet S, et al. Chemotherapy and radiotherapy in nasopharyngeal carcinoma: an update of the MAC-NPC meta-analysis. Lancet Oncol. 2015;16(6):645–655. doi:10.1016/S1470-2045(15)70126-9
67. Seelinger G, Merfort I, Schempp CM. Anti-oxidant, anti-inflammatory and anti-allergic activities of luteolin. Planta Med. 2008;74(14):1667–1677. doi:10.1055/s-0028-1088314
68. Ong C-S, Zhou J, Ong C-N, Shen H-M. Luteolin induces G1 arrest in human nasopharyngeal carcinoma cells via the Akt-GSK-3β-Cyclin D1 pathway. Cancer Lett. 2010;298(2):167–175. doi:10.1016/j.canlet.2010.07.001
69. Wu -C-C, Fang C-Y, Hsu H-Y, et al. EBV reactivation as a target of luteolin to repress NPC tumorigenesis. Oncotarget. 2016;7(14):18999–19017. doi:10.18632/oncotarget.7967
70. Khaiwa N, Maarouf NR, Darwish MH, et al. Camptothecin’s journey from discovery to WHO essential medicine: fifty years of promise. Eur J Med Chem. 2021;223:113639. doi:10.1016/j.ejmech.2021.113639
71. Li B-S, Huang J-Y, Guan J, Chen L-H. Camptothecin inhibits the progression of NPC by regulating TGF-β-induced activation of the PI3K/AKT signaling pathway. Oncol Lett. 2018;16(1):552–558. doi:10.3892/ol.2018.8688