This article is mentioned in:

Introduction

Nasopharyngeal carcinoma (NPC) is a common malignancy occurring in the head and neck. It is prevalent in the eastern and southeastern parts of Asia, especially in southern China, which accounts for 71% of new cases worldwide (1). In 2014, there were an estimated 44,600 new cases and 24,200 cases of NPC-associated mortality in China, and the crude incidence and mortality rates of NPC were 3.26/100,000 and 1.77/100,000, respectively (2). One of the possible reasons for this may be related to the occult characteristics of the occurrence site of NPC. NPC is commonly located in the pharyngeal recess (fossa of Rosenmüller), which is superior and posterior to the torus tubarius (3); therefore, many patients with NPC lack early characteristic symptoms and are diagnosed at an advanced stage. Although previous research has demonstrated that inherited factors, Epstein-Barr virus (EBV) infection and environmental factors are important causative factors of NPC (4), the molecular mechanisms in the pathogenesis of NPC have not yet been fully explained. Thus, studies in search of new targets for early diagnosis and treatment of NPC may be beneficial to improve patient quality of life and survival time.

At present, based on gene expression profiling, a variety of studies have been performed to explore the molecular mechanism underlying NPC. Using mRNA expression profiling, Bose et al (5) have demonstrated that loss of ATM serine/threonine kinase function may be an important step in the pathogenesis of NPC. Cai et al (6) discovered EBV-microRNA (miRNA)-binder of Arl2 induces metastasis of NPC by regulating phosphatase and tensin homolog-dependent pathways by miRNA profiling. With the help of long non-coding (lnc)RNA and mRNA expression profiling, ENST00000470135 was demonstrated to be upregulated in highly metastatic NPC cell lines and is related to the inhibition of migration, invasion and proliferation of NPC cells in vitro (7).

In the present study, an integrated bioinformatics approach was used to analyze the gene expression profiles of NPC using existing data. Differentially expressed genes (DEGs) between NPC and normal nasopharyngeal tissues were screened using R software. Distributions of DEGs in chromosomes were explored based on the annotation file and the CYTOBAND database of DAVID. The functions of DEGs were explored by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Additionally, a protein-protein interaction (PPI) network for DEGs was constructed to identify hub genes, key modules and important transcription factors (TFs). The results of the present study may be conducive to understanding the occurrence and development of NPC and may provide candidate molecular targets for diagnosis and treatment of NPC.

Materials and methods

Microarray data collection

NPC gene expression profiling datasets GSE12452 (8) and GSE34573 (9) were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) database of The National Center for Biotechnology Information. GSE12452 contained 31 NPC specimens and 10 normal nasopharyngeal tissues, whereas GSE34573 included 13 NPC specimens and 3 normal nasopharyngeal tissues excised by laser microdissection. The datasets were based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform. The annotation file (HG-U133_Plus_2 Annotations, CSV format, Release 36) is available for download from the Affymetrix website (https://www.thermofisher.com/cn/zh/home/life-science/microarray-analysis/microarray-data-analysis/genechip-array-annotation-files.html).

Data processing and DEG analysis

Data processing and statistical analysis were performed using packages from Bioconductor (http://www.bioconductor.org) for R software (version 3.4.1; http://www.r-project.org). A quality assessment was conducted by ‘affyPLM’ and ‘RColorBrewer’ packages to remove any disqualified samples (10,11). The subsequent processing steps included robust multiarray average normalization, calculation of missing values, and removal of batch effects using the ‘affy’, ‘impute’ and ‘sva’ packages (12–14). The expression data at the probes level were converted to gene symbols level data. The average value was used for genes with multiple probes. Subsequently, the limma package was used to screen DEGs between NPC specimens and normal nasopharyngeal tissues (15). The criteria for DEGs were false discovery rate <0.05 and |log fold change|>1. The DEGs were identified and their distributions on chromosomes were explored based on the annotation file and the CYTOBAND database of DAVID (v6.8; http://david.ncifcrf.go) (16).

GO terms and KEGG pathway enrichment analysis

The GO resource (http://www.geneontology.org) is an up-to-date biological database that developed a set of structured, controlled vocabularies to describe gene function analysis in terms of biological processes (BP), molecular function (MF) and cellular components (CC) (17). GO enrichment determines whether some gene sets defined by GO terms have high or low expression. KEGG (http://www.genome.ad.jp/kegg) is a knowledge database for systematic analysis of gene function, linking genomic information to higher order functional information (18). It stores graphical representations of cellular processes, such as metabolism, membrane transport, signal transduction and the cell cycle (18), and is widely used to identify clusters of co-expressed genes sharing the same pathway. GO and KEGG pathway enrichment analyses were performed to explore the biological function of DEGs using the ‘clusterProfiler’ package (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) in R (19). For GO analysis, the cut-off criteria were P<0.01 and q<0.05. For KEGG analysis, the cut-off criterion was P<0.05.

PPI network construction and network module analysis

To explore the functional associations of DEGs, a PPI network of NPC was predicted by the STRING database (version 10.5; http://www.string-db.org) (20). The highest confidence (0.9) was selected for the minimum required interaction score. Subsequently, the PPI network was visualized by Cytoscape software (version 3.6.0; http://www.cytoscape.org) (21). The Molecular Complex Detection (MCODE) app (22) was used to identify gene modules with the cut-off criterion of score ≥5. The ‘clusterProfiler’ package was used to compare the modules by their enriched GO terms and KEGG pathways. The scores of gene nodes were calculated by the CytoHubba app (23) to identify hub genes using 12 algorithms: Degree, Clustering Coefficient, Edge Percolated Component; Maximum Neighborhood Component, Density of Maximum Neighborhood Component, Maximal Clique Centrality, Bottleneck, Eccentricity, Closeness, Radiality, Betweenness and Stress. The 50 top-ranked genes were screened by each algorithm, and the genes common in >8 algorithms were considered to be candidate hub genes. The GSE13597 dataset, based on the platform of GPL96, submitted by Bose et al (5), was used to verify the differential expression levels of candidate hub genes between NPC and normal nasopharyngeal tissues.

Upstream transcriptional regulator analysis

To investigate the potential regulatory mechanism of DEGs, the iRegulon app (24), which performs motif discovery and track discovery, was used to predict upstream TFs targeting the DEGs. Normalized enrichment score (NES) >3 was used as the cut-off value.

Results

DEG analysis

A total of 906 DEGs were identified between NPC and normal nasopharyngeal tissues, including 434 (47.90%) upregulated genes and 472 (52.10%) downregulated genes. The top five significant genes ranked by adjusted P-value (Padj) were mucin 4, cell surface associated (MUC4; Padj=3.65×10−20), uroplakin 1B (Padj=2.33×10−16), tetratricopeptide repeat domain 9 (Padj=2.33×10−16), mucin 1 (MUC1; Padj=2.33×10−16) and β3 glycosyltransferase 7 (Padj=2.33×10−16). A heatmap of DEGs (top 25 upregulated and top 25 downregulated genes ranked by log fold change) is presented in Fig. 1. The DEGs were distributed on all human chromosomes except for the Y chromosome (Fig. 2A). Notably, among the DEGs, the location of LOC105371220 was undefined by the annotated file of the Affymetrix official website, so it could not be included in Fig. 2A. The peak count of DEGs was on chromosome 1 (117 DEGs), followed by chromosome 2 (76 DEGs) and chromosome 12 (60 DEGs). Chromosome 21 and chromosome 22 contained the fewest DEGs, 10 and 7, respectively. Chromosome 1 and chromosome 2 each contained 50 upregulated DEGs, which represented the peak of upregulated DEGs. The highest count of downregulated DEGs was in chromosome 1, which contained 67 downregulated DEGs. The CYTOBAND database of DAVID revealed that the DEGs were located on chromosome 7p15-p14 [homeobox A3, gamma-glutamylcyclotransferase, homeobox A10 (HOXA10), anillin actin binding protein (ANLN) and nuclear factor erythroid 2 like 3], 2q31 [FAST kinase domains 1, cell division cycle-associated 7, collagen type III α1 chain (COL3A1) and glutamate decarboxylase 1 (GAD1)], 1q21-q22 [small proline rich protein (SPRR)1A, SPRR1B, SPRR2B, SPRR3 and selenium binding protein 1 (SELENBP1)], 1q21 [MUC1, cysteine rich C-terminal 1, S100 calcium binding protein (S100)A8, S100A7, cingulin, S100A9, cornulin, S100A2 and S100A12], 4q21 [C-X-C motif chemokine ligand (CXCL)1, CXCL3, CXCL9 and CXCL10] and 1p31-p22 [chloride channel accessory (CLCA)2, CLCA4 and crystallin ζ) (P<0.025; Fig. 2B).

Figure 1. Heatmap of differentially expressed gene expression. Top 25 upregulated (red) and top 25 downregulated (green) genes. Orange, normal nasopharyngeal tissues; light gray, NPC specimens. NPC, nasopharyngeal carcinoma.

Figure 2. Chromosomal distribution of DEGs in nasopharyngeal carcinoma. (A) The peak count number of DEGs was the on chromosome 1, followed by chromosome 2 and chromosome 12. (B) CYTOBAND database of DAVID revealed that the DEGs were located on 7p15-p14, 2q31, 1q21-q22, 1q21, 4q21 and 1p31-p22; P<0.025. DEGs, differentially expressed genes.

GO terms and KEGG pathway enrichment analysis

To explore the biological function of DEGs, GO and KEGG pathway enrichment analyses were conducted using the ‘clusterProfiler’ package in R. The results revealed that the DEGs were significantly enriched in 68 GO terms for BP, 12 GO terms for MF and 44 GO terms for CC with the cut-off criteria of p<0.05 and q<0.05. For BP, the three enriched terms with the lowest P-value were ‘cilium movement’, ‘microtubule bundle formation’ and ‘axoneme assembly’ (Fig. 3A). For MF, the enriched terms included ‘peptidase regulator activity’, ‘motor activity’ and ‘microtubule motor activity’ (Fig. 3B). For CC, the enriched terms were ‘motile cilium’, ‘axoneme’ and ‘ciliary plasm’ (Fig. 3C). Additionally, KEGG pathway enrichment analysis indicated that five pathways were significantly enriched, including ‘cell cycle’ (hsa04110), ‘DNA replication’ (hsa03030), ‘interleukin (IL)-17 signaling pathway’ (hsa04657), ‘amoebiasis’ (hsa05146) and ‘glutathione metabolism’ (hsa00480) (Fig. 3D).

Figure 3. Top 10 enriched GO terms and top 5 KEGG pathways for differentially expressed genes. (A-C) GO term enrichment analysis for (A) biological process, (B) molecular function, (C) cellular component. (D) KEGG pathway analysis. Node size represents gene ratio; node color represents Padj. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; Padj, adjusted P-value.

PPI network construction and network module analysis

To investigate the functional associations of DEGs, a PPI network containing 867 nodes and 1,241 edges was constructed with the screened DEGs using the STRING database (Fig. 4). After analyzing the network with the MCODE app in Cytoscape software, eight modules with a score ≥5 were identified (Table I). Then, the DEGs in the modules were used to perform GO and KEGG pathway enrichment analysis, and to make a comparison using the ‘compareCluster’ function of the ‘clusterProfiler’ package. This differed from the aforementioned method of enrichment analysis using all DEGs. The results indicated that the modules, with the exception of module 6, which did not contain enough enriched genes in the analysis, were involved in different functional categories (Fig. 5). The results also demonstrated that cyclin-dependent kinase 1 (CDK1) and exportin 1 (XPO1) in module 1 were enriched in the pathway of EBV infection (Table II). Additionally, 14 candidate hub genes were identified using the CytoHubba app; the genes were distributed in module 1 [cyclin B1 (CCNB1), cell division cycle-associated 5 (CDCA5), CDK1, mitotic arrest deficient 2-like 1 (MAD2L1), nucleoporin 107 (NUP107)], module 4 [aurora kinase A (AURKA), cyclin A2 (CCNA2), protein kinase cAMP-dependent type II regulatory subunit β (PRKAR2B), DNA topoisomerase 2α (TOP2A), TPX2 microtubule nucleation factor (TPX2)], module 7 [cell division cycle 6 (CDC6), Rac GTPase-activating protein 1 (RACGAP1), ribonucleotide reductase regulatory subunit M2 (RRM2)] and module 8 [H2A histone family member Z (H2AFZ)]. The majority of these genes were upregulated in NPC in the GSE12452 and GSE34573 datasets (data not shown); however, in the GSE13597 dataset, only five genes (AURKA, CDC6, MAD2L1, TOP2A and TPX2) exhibited significant differences between NPC and normal nasopharyngeal tissues (Fig. 6). CDCA5 was not included in the Affymetrix annotation file of the GSE13597 dataset, so it was not included in Fig. 6. Therefore, five hub genes that might serve an important role in the initiation and development of NPC were obtained.

Figure 4. Protein-protein interaction network of differentially expressed genes constructed using Cytoscape. Purple nodes represent upregulated genes. Green nodes represent downregulated genes.

Figure 5. Enriched GO terms and KEGG pathways for seven modules. (A) GO term enrichment for biological processes, molecular functions and cellular components. (B) KEGG pathway analysis. The number on the × axis is the number of enriched genes. The circle size indicates the ratio of genes in each term; colors correspond to Padj. P-values were adjusted by the Benjamini-Hochberg method. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; M, module; Padj, adjusted P-value.

Figure 6. Expression levels of candidate hub genes in the GSE13597 dataset. Five candidate hub genes were upregulated in NPC (P<0.05). AURKA, aurora kinase A; CCNA2, cyclin A2; CCNB1, cyclin B1; CDC6, cell division cycle 6; CDK1, cyclin-dependent kinase 1; H2AFZ, H2A histone family member Z; MAD2L1, mitotic arrest deficient 2-like 1; NUP107, nucleoporin 107; PRKAR2B, protein kinase cAMP-dependent type II regulatory subunit β; RACGAP1, Rac GTPase-activating protein 1 (RACGAP1); RRM2, ribonucleotide reductase regulatory subunit M2; TOP2A, DNA topoisomerase 2α; TPX2, TPX2 microtubule nucleation factor.

Table I. Modules identified from the protein-protein interaction network of nasopharyngeal carcinoma.

Table I.

Modules identified from the protein-protein interaction network of nasopharyngeal carcinoma.

Module	Score	Nodes	Edges	Node IDs
1	19	19	171	KIF18A, MAD2L1, CENPE, ZWINT, CCNB1, MLF1IP, CENPH, NUP107, AHCTF1, XPO1, CENPK, CDCA5, CENPF, BUB1B, CDK1, NDC80, ZWILCH, NUF2, KNTC1
2	11	11	55	ANXA1, CCR1, CCR5, CCL5, C3, CXCL11, NMU, CXCL3, CXCL9, CXCL10, CXCL1
3	10.909	12	60	MUC13, MUC19, MUC4, B3GNT7, MUC5AC, GALNT11, MUC16, GALNT12, MUC20, GCNT3, MUC1, MUC15
4	9.524	22	100	KIF15, DLGAP5, PRC1, CETN2, HAUS1, HAUS6, TUBA1A, TOP2A, AURKA, TPX2, HMMR, PRKAR2B, TUBB4B, KIF23, KIF11, KIF13B, NEDD1, KIF18B, KIF4A, NCAPG, CCNA2, TTK
5	8	8	28	COL4A2, COL4A1, COL1A2, COL3A1, COL5A1, COL15A1, COL5A2, COL6A3
6	6.571	8	23	IFI44, IFIT1, XAF1, IFIT3, RSAD2, IFIT2, IFI44L, BST2
7	5.429	22	57	CDC6, CEP55, NUSAP1, CDC7, TOPBP1, PSMD14, KIF20A, ASPM, ANLN, MELK, DBF4, PRIM1, RACGAP1, ORC6, MCM3, MCM7, TIPIN, PSMD12, CDC45, PBK, PSMA1, RRM2
8	5	9	20	BRIP1, PCNA, H2AFZ, DNA2, HIST2H2BE, SMC4, RBBP8, RMI2, RAD51AP1

Table II. Enriched KEGG pathways in module 1.

Table II.

Enriched KEGG pathways in module 1.

ID	Description	Count	P-value	Gene symbol
hsa04110	Cell cycle	4	2.48×10−6	CCNB1/BUB1B/CDK1/MAD2L1
hsa04914	Progesterone-mediated oocyte maturation	3	7.70×10−5	CCNB1/CDK1/MAD2L1
hsa04114	Oocyte meiosis	3	1.51×10−4	CCNB1/CDK1/MAD2L1
hsa05166	Human T-cell leukemia virus 1 infection	3	1.26×10−3	XPO1/BUB1B/MAD2L1
hsa04115	p53 signaling pathway	2	1.88×10−3	CCNB1/CDK1
hsa04218	Cellular senescence	2	8.99×10−3	CCNB1/CDK1
hsa03013	RNA transport	2	1.02×10−2	NUP107/XPO1
hsa05169	Epstein-Barr virus infection	2	1.39×10−2	XPO1/CDK1
hsa05170	Human immunodeficiency virus 1 infection	2	1.54×10−2	CCNB1/CDK1

Upstream transcriptional regulator analysis

The iRegulon app in Cytoscape was used to explore the potential regulatory mechanism of DEGs. A total of 14 TFs predicted to regulate 736 DEGs were identified (Table III). Among them, regulatory factor X4 (RFX4) regulated more DEGs compared with the others (n=362), but the TF with the highest NES score was E2F transcription factor (E2F) 4 (NES score=6.925). In addition, forkhead box M1 (FOXM1), E2F6 and E2F5 were DEGs. Additionally, the results revealed that 10 TFs [BCL6 transcription repressor, E2F1, E2F2, E2F4, E2F7, FOXM1, HNF1 homeobox B, SIN3 transcription regulator family member A (SIN3A), sex determining region Y box 8 (SOX8) and transcription factor Dp-1 (TFDP1)] were predicted to regulate 14 candidate hub genes. The possible regulatory relationships between the TFs and the candidate hub genes are presented in Fig. 7.

Figure 7. Possible regulatory relationships between the TFs and the hub genes of nasopharyngeal cancer. Green represents TFs, purple represents candidate hub genes. TF, transcription factor.

Table III. Transcription factors targeting differentially expressed genes of nasopharyngeal carcinoma predicted by iRegulon.

Table III.

Transcription factors targeting differentially expressed genes of nasopharyngeal carcinoma predicted by iRegulon.

ClusterCodea	Transcription factor	NES scoreb	Target genesc (n)	Motifs/tracks (n)
T1	E2F4	6.925	173	3
M1	E2F5	6.301	217	26
T2	TFDP1	6.258	107	1
T3	FOXM1	4.926	161	5
T4	E2F7	4.824	47	1
T5	SIN3A	4.600	187	2
M3	RFX4	4.485	362	8
T6	E2F1	3.945	87	1
M4	NFKB1	3.935	320	8
M5	SOX8	3.837	252	3
M8	E2F6	3.377	33	1
T7	E2F2	3.177	196	1
M10	BCL6	3.122	140	2
M13	HNF1B	3.003	161	1

a Each cluster code contains enriched motifs that are clustered by similarity or contains enriched tracks that are clustered by transcription factor name.

b NES represents the highest enrichment score for a motif/track that is part of the current cluster code.

c Targets represent the number of motifs/tracks that are be associated with the given transcription factor.

Discussion

Gene chips have been increasingly used in NPC research; however, the DEGs screened by different platforms or researchers vary significantly. Even studies on the same platform may yield different results, which may be associated with the differences in specimen sources, sample processing and experiment batches. Integration and analysis of gene chip data to obtain more consistent results has become problematic. Some researchers choose to do meta-analysis of multiple gene expression profiles based on different platforms or from a number of types of cells, whereas others chose to integrate data based on the same platform (25–27). Additional molecular biology experiments may be needed to conclude which method provides more accurate results. Prior to the present study, a number of NPC datasets in GEO databases were selected to perform DEG filtering separately; using the same method, the number of screened DEGs from two datasets (GSE12452 and GSE34573) based on the same platform was similar (data not shown). GSE12452 is a classical gene expression profile of NPC. GSE34573 was submitted recently, and to the best of our knowledge, no other researchers have integrated GSE34573 with other chips for bioinformatics research; in the present study, GSE12452 and GSE34573 were chosen as the research objects. In addition, the annotations to the gene chips and the databases of gene functions and pathways, including GO and KEGG databases, are constantly updated. More meaningful results may be obtained by using the updated database for additional gene function analysis in the future.

In this study, the raw data of GSE12452 and GSE34573 from GEO database were downloaded. After data processing and statistical analysis with R software, 906 DEGs were obtained. These DEGs were distributed on all human chromosomes except the Y chromosome. Chromosome 1 contained the largest number of DEGs, while chromosome 21 and chromosome 22 contained fewer DEGs. The DEGs were identified to be enriched on chromosome 7p15-p14, 2q31, 1q21-q22, 1q21, 4q21 and 1p31-p22 based on the CYTOBAND database of DAVID. Among the genes which were located on these chromosomal regions, HOXA10, ANLN, COL3A1, GAD1, SELENBP1, S100A8, S100A9, S100A2, CXCL1, CXCL9, CXCL10 and CLCA2 have been reported to be associated with NPC. In the original publications for GSE12452, only two chromosomal regions, 4p15.1-4q12 and 14q32-14q32.33, which were previously demonstrated to be linked to NPC, were evaluated (8,28). The genes in these two regions were found to be differentially expressed between tumor and normal tissue specimens. However, other chromosome regions were not identified in the previous publications. GSE34573 was designed to examine both genetic copy number changes and gene expression changes in the same samples (9). In the original publication for GSE34573, 1,049 a priori determined tumor-related genes were examined for differential expression (9), but in the present study, all genes on the chip were examined.

To explore the biological functions of the DEGs, GO and KEGG enrichment analyses were conducted. The enriched GO terms and pathways demonstrated a number of differences between NPC specimens and normal nasopharyngeal tissues. The human nasopharynx is lined with ciliated pseudostratified columnar epithelium, but the cilia are absent in NPC cells (29). In the present study, GO terms ‘motile cilium’, ‘ciliary plasm’ and ‘axoneme’ were enriched, and the associated genes exhibited a consistent pattern of low expression (data not shown). In addition, the ‘IL-17 signaling pathway’ was enriched in the KEGG analysis. IL-17 is a proinflammatory cytokine predominantly produced by T-helper 17 cells. In NPC, IL-17A has been reported to promote the migration and invasion of NPC cells via the p38 mitogen-activated protein kinase/NF-κB) signaling pathway, which subsequently upregulates the expression of matrix metallopeptidase (MMP) 2 and MMP9 and enhances epithelial-mesenchymal transition (EMT) (30).

To reveal the interactions between DEGs of NPC, a PPI network was constructed using the STRING database. By means of the MCODE app in Cytoscape software, eight modules that might serve an important role in the development of NPC were detected. Subsequently, GO and KEGG pathway enrichment analyses were performed to reveal the biological functions of these modules. The results demonstrated that the genes in the seven modules displayed different functions. For example, in module 2, the enriched GO terms included ‘chemokine activity’, ‘G-protein coupled receptor binding’ and ‘chemokine receptor binding’, and the enriched KEGG pathways included ‘Chemokine signaling pathway’, ‘Cytokine-cytokine receptor interaction’ and ‘Kaposi sarcoma-associated herpesvirus infection’; while in module 5, the enriched GO terms included ‘extracellular matrix structural constituent’, ‘extracellular matrix structural constituent conferring tensile strength’ and ‘platelet-derived growth factor binding’, and the enriched KEGG pathways included ‘Protein digestion and absorption’, ‘ECM-receptor interaction’ and ‘Amoebiasis’. In addition, the DEGs in module 1 were similar to previously identified modules (27), and ‘cell cycle’, ‘cell division’, ‘oocyte meiosis’ and ‘p53 signaling pathway’ were also identified to be enriched in enrichment analysis; however, compared with previous studies (25,27), more modules and information regarding GO terms and KEGG pathways of the modules were provided in the present study. Of note, CDK1 and XPO1 were demonstrated to be involved in the pathway of EBV infection. CDK1 encodes a catalytic subunit of M-phase promoting factor, which is necessary for G1/S and G2/M phase transitions during the eukaryotic cell cycle (31). EBV nuclear antigen (EBNA) leader protein (LP) and EBNA1 are important for EBV-induced B cell immortalization; it has been reported that CDK1 mediates the phosphorylation of EBNA1 and EBNA-LP at Ser35 and Ser393 sites, respectively (32,33). The protein encoded by XPO1 is involved in the nuclear export of cellular proteins; it binds to EBV SM protein, the spliced BSLF2/BMLF1 protein (34), which is an essential process for activity and cytoplasmic localization of SM protein (35). Therefore, these results suggested that CDK1 and XPO1 may be important genes in the pathogenic process of NPC induced by EBV.

There are a number of algorithms used for the identification of the hub genes. To improve the reliability of the results of the present study, the top 50 genes of 12 algorithms were investigated with the CytoHubba app. The genes common in >8 algorithms were regarded as the hub genes. The GSE13597 dataset was used to determine the expression levels of the candidate hub genes. A total of 14 candidate hub genes of NPC were identified. Among these, seven have been reported previously: AURKA, CCNA2, CCNB1, CDC6, CDK1, MAD2L1 and TOP2A (27). The other seven genes were identified as candidate hub genes for the first time: CDCA5, H2AFZ, NUP107, PRKAR2B, RACGAP1, RRM2 and TPX2. H2AFZ encodes protein H2A.Z.1, which is an isoform of H2A.Z; as a protein incorporated into chromatin, H2A.Z is associated with the orchestration of DNA double-strand break repair, maintaining the integrity of telomeres, genome stability and chromosome segregation (36). In hepatocellular carcinoma, H2A.Z.1 serves a pivotal role in liver tumorigenesis by selectively regulating key molecules in the cell cycle and EMT (37). NUP107, a component of the nuclear pore complex, regulates the exchange of molecules between the nucleoplasm and cytoplasm in a selective manner (38). A variant in the 3′-UTR region of NUP107 has been identified to be associated with sensitivity to platinum-based chemotherapy in ovarian cancer (39). RACGAP1 encodes a GTPase-activating protein, which is a component of the central spindlin complex; it is involved in tumor cell proliferation, migration, invasion and metastasis and serves an oncogenic role in a number of malignancies such as breast cancer, hepatocellular carcinoma and colorectal cancer (40). RRM2 encodes a subunit of ribonucleotide reductase, which catalyzes the conversion of ribonucleotides into deoxyribonucleotides; overexpression of RRM2 predicts a poor prognosis for patients with NPC (41). The protein encoded by TPX2 is associated with the assembly of microtubules and mitotic spindles during mammalian cell mitosis; it is necessary for targeting AURKA to the microtubules of the mitotic spindle, and elevated TPX2 expression is correlated with the aggressiveness of ovarian and salivary gland cancer (42). However, among the 14 candidate hub genes, only five (AURKA, CDC6, MAD2L1, TOP2A and TPX2) were demonstrated to have expression differences between NPC and normal nasopharyngeal tissues, which may be associated with the small number of samples in the GSE13597 dataset, which contained 25 NPC specimens and three normal nasopharyngeal tissues. Nevertheless, except for CDCA5 and PRKAR2B, the candidate hub genes exhibited a tendency of high expression in NPC.

In addition, 14 TFs regulating the DEGs of NPC were identified in the present study. Among these, FOXM1, E2F1, TFDP1 and NF-κB subunit 1 (NFKB1) have been previously studied in NPC. FOXM1 is upregulated in NPC cells and related to metastases in patients with NPC (43). Downregulation of FOXM1 by thiostrepton or small interfering RNA leads to the inhibition of NPC cell proliferation, migration, invasion, transformation ability and angiogenesis (44,45). Several signaling pathways of E2F1, including p21/CCND1/CDK6/E2F1, p15/p16/CDK4/E2F1, and p16 (INK4A)/Rb/E2F1, have been reported to be involved in carcinogenesis of NPC (46,47). Moreover, transactivation of E2F1 and/or TFDP1 might contribute to overexpression of stathmin 1 in NPC (48). Aberrant NF-κB signaling has been demonstrated to be constitutively activated by genetic alterations or EBV-encoded latent membrane protein 1 (49,50). The remaining TFs (E2F4, E2F5, E2F7, SIN3A, RFX4 and SOX8) were reported to be associated with NPC for the first time in the present study.

In summary, based on an integrated bioinformatics approach, five hub genes and 14 TFs were identified in NPC. Among them, TPX2 were reported to be associated with the pathogenesis of NPC for the first time. In addition, pathway enrichment analysis for modules suggested that CDK1 and XPO1 may be important genes in EBV-induced NPC. The pathogenesis of NPC is complex and involves the majority of human chromosomes, and EBV may be involved. Additionally, the results of the enrichment analysis revealed that, in the course of NPC, abnormalities of cell morphology, glutathione metabolism, immunity and the formation and movement of cilium, microtubules and axoneme are affected. In addition, DNA replication, cell cycle, chromosome separation and other processes were involved, which ultimately lead to malignant transformation of nasopharyngeal cells.

It must be noted that there were certain limitations to the present study. Firstly, only two datasets were integrated and analyzed, and the number of cases was small. Secondly, the distribution of DEGs on different chromosomes was studied at the level of mRNA, but chromosomal aberrations were not investigated. Thirdly, the present study only validated the expression of hub genes using datasets based on GSE13597 and lacks validation experiments at the molecular level. However, these results may promote the understanding of the molecular mechanisms underlying NPC and provide useful clues for identifying molecular targets for diagnosis and treatment of NPC in the future.

Acknowledgements

Not applicable.

Funding

The present study was supported by The Scientific and Technological Project of Zhanjiang (grant no. 2012C3105022) and The Foundation for Young Scholars of Guangdong Medical College (grant no. 2012004).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

JH, ZZ and SG conceived and designed the experiments. ZZ, SG, SL and RL performed the experiments. ZZ and SG analyzed the data. SL and RL participated in the acquisition of data, prepared figures and/or tables and interpreted the data. ZZ, SG and SL wrote the paper.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References