EBV-Associated Hub Genes as Potential Biomarkers for Predicting the Prognosis of Nasopharyngeal Carcinoma

Abstract

This study aimed to develop a model using Epstein-Barr virus (EBV)-associated hub genes in order to predict the prognosis of nasopharyngeal carcinoma (NPC). Differential expression analysis, univariate regression analysis, and machine learning were performed in three microarray datasets (GSE2371, GSE12452, and GSE102349) collected from the GEO database. Three hundred and sixty-six EBV-DEGs were identified, 25 of which were found to be significantly associated with NPC prognosis. These 25 genes were used to classify NPC into two subtypes, and six genes (C16orf54, CD27, CD53, CRIP1, RARRES3, and TBC1D10C) were found to be hub genes in NPC related to immune infiltration and cell cycle regulation. It was shown that these genes could be used to predict the prognosis of NPC, with functions related to tumor proliferation and immune infiltration, making them potential therapeutic targets. The findings of this study could aid in the development of screening and prognostic methods for NPC based on EBV-related features.

Keywords: Epstein–Barr virus; hub genes; nasopharyngeal carcinoma; prognosis.

Figure 1

Identification of genes associated with EBV infection in the GSE2371 dataset. (A,B) The gene expression profile of the raw data from GSE2371 before and after normalization. The bars represent all samples included in the analysis, with each color indicating a different sample. (C) Heatmap of the representative 20 up- and 20 downregulated EBV-associated genes. Pink and cyan indicate EBV-positive and EBV-negative NPC cell lines, respectively. Red and blue indicate up- and downregulated EBV-DEGs, respectively (p < 0.05). (D) GO term enrichment of the EBV-DEGs. BP, biological process; CC, cellular components; MF, molecular function. (E) KEGG pathway enrichment of the EBV-DEGs.

Figure 2

Screening of prognosis-associated genes in the combined datasets of GSE12452 and GSE102349. (A,B) Principal component analysis (PCA) of the combined datasets ((A) without removing batch effects; (B) removing batch effects). (C) Correlation network of 25 prognosis-related EBV-DEGs. Purple and cyan dots indicate risk and favorable factors, respectively. Pink and blue lines indicate positive and negative correlations, respectively, between two genes (p < 0.0001).

Figure 3

Kaplan–Meier (KM) survival curves of 25 prognosis-associated EBV-DEGs.

Figure 4

Unsupervised clustering of 25 EBV-associated genes. (A) Consensus matrix heatmap. The optimal number of clusters: k = 2. (B) Kaplan–Meier (KM) survival curves showing the survival probability of the two EBV-associated prognosis clusters (p = 0.012). (C) The expression level of 25 prognosis-associated EBV-DEGs in the two EBV-associated prognosis clusters (*, p < 0.05; **, p < 0.01; ***, p < 0.001; ns, no significance). (D) Heatmap showing the expression level of 25 prognosis-associated EBV-DEGs, the survival outcomes (alive or dead), and full survival times (days) in the two EBV-associated prognosis clusters.

Figure 5

Pathway enrichment analysis in the two EBV-associated prognosis clusters. (A) HALLMARK, (B) KEGG, and (C) Reactome pathway differences in NPC prognosis cluster A and B. Blue, cluster A; yellow, cluster B.

Figure 6

Immune infiltration analysis in the two EBV-associated prognosis clusters. (A) Principal component analysis (PCA) results are visualized using scatter dot plots, which depict the samples in the two clusters. (B) Differences in StromalScore, ImmuneScore, and ESTIMATEScore in the two clusters. (C) Differences in 23 types of immune-infiltrating cell enrichment in the two clusters. Blue, cluster A; yellow, cluster B; ***, p < 0.001; ns, no significance.

Figure 7

Overview of gene expression and pathway enrichment in EBV-associated low-/high-risk NPC patients. (A) Volcano plots of DEGs between cluster A and cluster B. The threshold is set at logFC > 1 and p < 0.05. (B) GO enrichment items. (C) KEGG enrichment items. (D) The correspondence between the top five KEGG enriched items and the DEGs.

Figure 8

Identification of the DEGs between NP and NPC based on DEGs in EBV-associated low-/high-risk NPC patients. (A) Forest plot of 48 prognosis-associated DEGs in NPC when the significance is set at p < 0.005. (B,C) Gene expression levels in the merged dataset without or after debatching. The bars represent all samples included in the analysis, with each color indicating a different sample. (D) Volcano plots showing DEGs between NP and NPC samples in the merged dataset. (E) Heatmap showing the top 20 up- and downregulated genes. (F) Venn diagram of 3337 upregulated genes and 48 prognosis-associated DEGs. (G) Venn diagram of 2501 downregulated genes and 48 prognosis-associated DEGs. (H) Box plot of 11 downregulated genes. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001; ns, no significance.

Figure 9

Identification of hub genes using machine learning. (A) Construction of prognosis-associated candidate genes by random forest. (B) Screening of prognosis-associated candidate genes by SVM. (C) Construction of an occurrence-associated candidate gene signature by random forest. (D) Screening of occurrence-associated candidate genes by SVM. (E) Venn diagram showing the overlap of the candidate genes in (A–D). (F) Circos plot displaying the relationship between the overlapping genes (hub genes) in E. (G) ROC curve of the NPC prognosis signature. (H) ROC curve of the NPC occurrence signature.

Figure 10

Hub gene-related immune infiltration analysis. (A) Comparison of 23 immune cell subtypes between NC and NPC patients. (B) Correlation matrix of all 23 immune cell subtype compositions in NC and NPC patients. (C) Correlation between hub genes and immune-infiltrating cells. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ns, no significance.

Figure 11

Correlation analysis between the six hub genes and all genes expressed in NPC. The top 50 positive and negative correlation genes are displayed in the heatmaps.

Figure 12

GSEA of the six hub genes. The top 20 items are displayed in the enrichment map.

Figure 13

The miRNA and transcription factor regulatory network of the six hub genes.