Identification of markers for predicting prognosis and endocrine metabolism in nasopharyngeal carcinoma by miRNA-mRNA network mining and machine learning

Abstract

Background: Nasopharyngeal cancer (NPC) has a high incidence in Southern China and Asia, and its survival is extremely poor in advanced patients. MiRNAs play critical roles in regulating gene expression and serve as therapeutic targets in cancer. This study sought to disclose key miRNAs and target genes responsible for NPC prognosis and endocrine metabolism.

Materials and methods: Three datasets (GSE32960, GSE70970, and GSE102349) of NPC samples came from Gene Expression Omnibus (GEO). Limma and WGCNA were applied to identify key prognostic miRNAs. There were 12 types of miRNA tools implemented to study potential target genes (mRNAs) of miRNAs. Univariate Cox regression and stepAIC were introduced to construct risk models. Pearson analysis was conducted to analyze the correlation between endocrine metabolism and RiskScore. Single-sample gene set enrichment analysis (ssGSEA), MCP-counter, and ESTIMATE were performed for immune analysis. The response to immunotherapy was predicted by TIDE and SubMap analyses.

Results: Two key miRNAs (miR-142-3p and miR-93) were closely involved in NPC prognosis. The expression of the two miRNAs was dysregulated in NPC cell lines. A total of 125 potential target genes of the key miRNAs were screened, and they were enriched in autophagy and mitophagy pathways. Five target genes (E2F1, KCNJ8, SUCO, HECTD1, and KIF23) were identified to construct a prognostic model, which was used to divide patients into high group and low group. RiskScore was negatively correlated with most endocrine-related genes and pathways. The low-risk group manifested higher immune infiltration, anticancer response, more activated immune-related pathways, and higher response to immunotherapy than the high-risk group.

Conclusions: This study revealed two key miRNAs that were highly contributable to NPC prognosis. We delineated the specific links between key miRNAs and prognostic mRNAs with miRNA-mRNA networks. The effectiveness of the five-gene model in predicting NPC prognosis as well as endocrine metabolism provided a guidance for personalized immunotherapy in NPC patients.

Keywords: endocrine; immune checkpoint blockade; immunotherapy; miRNA-mRNA network; micro RNAs; nasopharyngeal cancer; risk model.

Figure 1

Identification of NPC-related DEmiRNAs in the GSE32960 dataset. (A) Volcano plot of 332 DEmiRNAs. (B) Clustering dendrogram of 312 samples. (C, D) Under different soft thresholds (power), we analyzed scale independence and mean connectivity. (E) Clustering dendrogram of DEmiRNAs based on TOM and dynamic cut methodology. (F) The relationships of modules with normal and tumor groups. (G) Venn plot of DEmiRNAs and miRNAs of brown and blue modules.

Figure 2

ROC analysis and survival analysis for evaluating the performance of the miRNA risk model. (A) ROC curves and survival curves in the training group. (B) ROC curves and survival curves in the test group. (C) ROC curves and survival curves of DFS, OS, MFS, and RFS in the GSE32960 dataset. (D) ROC curves and survival curves of OS and DFS in the GSE70970 dataset.

Figure 3

Kaplan–Meier survival analysis of two risk groups in the samples with different clinical characteristics (A–J).

Figure 4

Constructing a nomogram based on clinical features and risk score. (A, B) Univariate (A) and multivariate Cox regression analyses of clinical features and risk type. (C) The nomogram based on gender and risk score for predicting 1-, 3-, and 5-year survival. (D) Calibration curve of 1-, 3-, and 5-year survival. (E) DCA of nomogram, gender, and risk score. *P < 0.05, ***P < 0.001.

Figure 5

Analysis of the target genes of has-miR-142-3p and has-miR-93. (A) The mRNA–miRNA ceRNA networks. Green rhombus indicates target mRNAs, and ellipse indicates miRNAs. (B) KEGG and (C–E) GO functional analyses of potential target mRNAs. The color of dots indicate the significance of P values, and the dot size indicates the gene counts.

Figure 6

Construction of the prognostic model based on the target genes in the GSE102349 dataset. (A) The top 15 target genes associated with prognosis from 1,000-times random sampling. (B) ROC analysis and survival analysis of the five-gene prognostic model.

Figure 7

Immune characteristics of high-risk and low-risk groups in the GSE102349 dataset. (A) The ssGSEA score of 22 immune-related cells in two risk groups. (B) The ssGSEA score of adaptive and innate immune response in two risk groups. (C) ESTIMATE analysis of immune infiltration and stromal infiltration. (D) MCP-counter analysis for estimating the enrichment of 10 immune-related cells. (E) Expressions of immune checkpoint genes in two risk groups. ns, not significant. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Figure 8

Analysis of biological pathways, response to immunotherapy, and drug sensitivity of two risk groups in the GSE102349 dataset. (A) Screening differentially enriched pathways between two risk groups based on ssGSEA. (B) Correlation analysis of risk score with 13 tumor-related pathways. Orange and blue indicate positive and negative correlations, respectively. (C) SubMap analysis of expression data of GSE102349 and IMvigor210 (anti-PD-L1 treatment) datasets. (D) TIDE analysis of two risk groups for assessing immune escape and response to immunotherapy. (E) The sensitivity to chemotherapeutic drugs predicted by the pRRophetic package. *P < 0.05, **P < 0.01, ***P < 0.001, ****P<0.0001.

Figure 9

Analysis of endocrine metabolism. (A) Pearson analysis between endocrine-related genes and RiskScore. (B) Pearson analysis between endocrine pathways and RiskScore.

Figure 10

The performance of the mRNA prognostic model in three independent immunotherapy datasets. (A–C) ROC analysis, survival analysis, the distribution responder groups in two risk groups, and the risk score of responder groups in the IMvigor210 (A), GSE135222 (B), and GSE78220 (C) datasets. ns P>0.05, *P < 0.05, **P < 0.01, ***P < 0.001.