Skip to main content

ORIGINAL RESEARCH article

Front. Endocrinol., 19 July 2023
Sec. Systems Endocrinology
This article is part of the Research Topic Machine Learning-Assisted Diagnosis and Treatment of Endocrine-Related Diseases View all 13 articles

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

Xixia ZhangXixia Zhang1Xiao LiXiao Li2Caixia WangCaixia Wang2Shuang WangShuang Wang2Yuan ZhuangYuan Zhuang2Bing LiuBing Liu1Xin Lian*Xin Lian1*
  • 1Department of Otolaryngology Head and Neck Surgery, Shengjing Hospital of China Medical University, Shenyang, China
  • 2Department Obstetrics and Gynecology, Shengjing Hospital of China Medical University, Shenyang, China

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.

Introduction

Nasopharyngeal cancer (NPC) is an epithelial malignancy, which has discrepant occurrence in different regions and countries. The etiology of NPC is multiple and remains incompletely understood, but most cases are closely linked to Epstein–Barr virus (EBV) infection (1). In Western countries, the incidence rate is relatively low with a ratio of 1:100,000 each year. However, in the regions of Southern China and Asia, the annual incidence elevates to 25–50 cases per 100,000 (2), which accounts for approximately 70% new cases worldwide (3). The incidence also strikingly increases in the recent decades in China. The age-standardized incidence rate (ASIR) was 3.3/100,000 in 1990, whereas it reached 5.7/100,000 in 2019. The rising incidence rate was especially startling in men, with ASIR of 4.3/100,000 to 8.6/100,000 from 1990 to 2019 (4). The distant metastasis contributes to an extremely poor prognosis in NPC patients, and the patients of stages III and IV have a 5-year survival rate less than 10%.

With intensive usage of modulated treatment including chemotherapy and radiotherapy, the control of NPC metastasis reaches a satisfactory outcome and the 5-year survival rate drastically improves (57). Nevertheless, the prognosis and treatment efficiency were awfully unfavorable in the NPC patients of late stages. It is still a great challenge for clinicians to control and lessen the metastasis and recurrence of advanced NPC patients. In the recent years, immunotherapy reaches a milestone in clinical cancer therapy, not excluding in NPC. Immune checkpoint blockade (ICB) therapy is one of the promising strategies to increase antitumor activity in NPC. Studies have shown that NPC patients expresses high expression levels of programmed death protein 1 (PD-1) and programmed death ligand 1 (PD-L1) that are associated with poor outcomes and recurrence (8, 9). The blockade of PD-1/PD-L1 expression can recover the ability of cytotoxic lymphocytes exerting anticancer response. Lines of clinical trials of ICB therapy have demonstrated the positive response to anti-PD-1/PD-L1 drugs such as pembrolizumab in phase 1 and 2 studies (10, 11). In the phase 2 study, of 44 enrolled NPC patients, eight patients reached a partial response and one patient reached a complete response, showing an objective response rate (ORR) of 20.5% (11). As we know, in the tumor microenvironment, the expression of PD-1 on the surface of tumor cells and the combination of PD-L1 on the surface of tumor-infiltrating lymphocytes can inhibit the activity of T cells, lead to the loss of function of effector factors TNF-a, IFN-gamma, and IL-2, and inhibit cytolysis function through granulozyme B and perforin, and ultimately promote immune escape (12, 13). Obviously, still a high proportion of NPC patients showed a negative response to ICB therapy, which may result from their disadvantageous tumor microenvironment. Therefore, in order to raise the treatment accuracy, understanding the mechanism of immune evasion and developing molecular biomarkers is critically needed.

Over the last few decades, it has become clear that endocrines are also involved in regulating tumor cells and that cancer cells themselves abnormally express and respond to many hormones (14). Both leptin and its receptor have been reported in cancer biopsy specimens, indicating autocrine and/or paracrine roles in tumorigenesis (15). Several hypothalamic hormones have been implicated in a variety of human cancers, including growth hormone-releasing hormone, luteinizing hormone-releasing hormone, somatostatin, and bombotin (16, 17). It is becoming increasingly clear that many human cancer cells are sensitive to a variety of hormones and that they themselves express many hormones that play an important role in the development and progression of cancer.

Non-coding RNAs play essential roles in gene modulation and pathway regulation. Some microRNAs (miRNAs) were unveiled to participate in NPC invasion and metastasis, immune escape, and resistance to chemotherapy and radiotherapy (18), endowing miRNAs possible to serve as potential therapeutic targets (19). For example, miR-26a was found to have an anticancer effect and overexpressing miR-26a could inhibit the metastatic feature in NPC cells (20). MiR-663 targeting WAF1/CIP1 promotes the proliferation and tumorigenesis in NPC cells (21). The crosstalk between extracellular microRNAs and the tumor microenvironment has also been profoundly parsed by previous studies (22, 23). Consequently, our study tried to mine the key miRNAs that had a prognostic value in NPC. We identified two potentially key miRNAs (has-miR-142-3p and has-miR-93) closely involved in NPC prognosis. By building competing endogenous RNA (ceRNA) networks using multiple miRNA tools, we determined 125 potential target genes (mRNAs) and screened five key genes contributable for the prognostic model. The five-gene prognostic model manifested a satisfactory performance in prognosis prediction and estimating the response to immunotherapy and chemotherapy.

Materials and methods

Data source and preprocessing

GSE32960, GSE70970, and GSE102349 datasets containing the expression profiles of NPC samples were accessed from the Gene Expression Omnibus (GEO) database (24), where GSE32960 and GSE70970 include miRNA expression data and GSE102349 includes mRNA expression data. We screened the samples of GEO datasets according to the following criteria: 1) remove samples without survival time and survival status; 2) convert probes into gene symbols; 3) remove a probe matching to multiple genes; 4) select the averaged expression of a gene with multiple probes; 5) for miRNA data, only human samples were remained. After preprocessing, 312 NPC and 18 normal samples were remained in the GSE32960 dataset; 253 NPC and 10 normal samples were remained in the GSE70970 dataset; and 88 NPC samples were remained in the GSE102349 dataset.

Identification of NPC-associated differentially expressed miRNAs

In the GSE32960 dataset, differentially expressed miRNAs (DEmiRNAs) were screened by Limma R package (25) from NPC and normal samples under conditions of P < 0.05 and |log2(fold change, FC)| > 1.2. Gene modules were identified by weighted correlation network analysis (WGCNA) (26). Firstly, samples were clustered and the co-expression network was constructed. A scale-free network was ensured under scale-free R2 = 0.85, and soft threshold (power) = 3 was determined. The co-expression network was then converted to the adjacent matrix and was further converted to the topological overlap matrix (TOM). Subsequently, we clustered genes by the dynamic cutting method and average-linkage hierarchical clustering based on TOM. Eigengenes were calculated for each gene, and gene modules were clustered and merged under parameters of deepSplit = 2, and minModuleSize = 60, height = 0.25. The module–trait relationships were assessed by Pearson correlation analysis. By overlapping the miRNAs in the NPC-associated gene modules and DEmiRNAs, NPC-associated miRNAs were determined.

Construction of a miRNA-based risk model

Random grouping of samples from the GSE32960 dataset into training group and test group at a ratio of 3:2 was performed. Two-group differences were assessed using Student’s t test. Univariate Cox regression analysis screened NPC-associated miRNAs from the training group. MiRNAs with P < 0.01 were selected as prognostic miRNAs. Multivariate analysis (stepAIC) was introduced to measure the coefficients of prognostic miRNAs. Then, the miRNA-based risk model was defined as: risk score = Σ(coef i*expression i), where coef indicates the coefficients of miRNAs and i represents miRNAs.

Evaluation and optimization of the risk model

Each sample obtained a risk score calculated by the risk model. The median risk score was employed in dividing samples into low-risk and high-risk groups. Receiver operating characteristic (ROC) curve analysis was used to predict the efficiency of the risk model in predicting overall survival through timeROC R package (27). The prognosis difference of two risk groups was studied by Kaplan–Meier survival analysis. Univariate and multivariate Cox regression models were used to analyze the hazard ratio of risk type. A nomogram was used to optimize the clinical use of the risk model with the rms package.

Construction of a mRNA-based prognostic model

First of all, the potential target genes of miRNAs were predicted by different online tools and software including microT (28), miRanda (29), mircode (30), miRDB (31), miRmap (32), miRtarbase (33), PicTar (34), PITA (35), TargetMiner (36), TargetScan (37), RNA22 (38), and starbase (39). The target genes predicted by at least six tools were remained as key potential target genes. The WebGestaltR package (40) was utilized to annotate the enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms. Then, we used the key target genes to establish the mRNA-based prognostic model. Prognostic target genes in 70% samples of the GSE102349 dataset were screened by univariate Cox regression analysis. Random sampling with 70% samples in the GSE102349 dataset was performed for 1,000 times corresponding with univariate analysis. The top five frequent target genes from the results of 1,000 times of univariate analysis were selected as the final target genes for constructing the mRNA-based prognostic model (defined by the same formula with the miRNA risk model).

Endocrine metabolism analysis

There were 31 SECRETORY_PATHWAYS screened from the KEGG PATHWAY Database (https://www.genome.jp/kegg/pathway.html). SECRETORY_PATHWAYS scores in the GSE102349 dataset were calculated by ssGSEA (41). There were 26 SECRETORY-related genes obtained from KEGG (https://www.genome.jp/dbget-bin/www_bfind_sub?mode=bfind&max_hit=1000&locale=en&serv=kegg&dbkey=genes&keywords=Secretory&page=1). Pearson analysis was conducted to determine the correlation between SECRETORY_PATHWAYS scores/SECRETORY-related genes and RiskScore.

Analysis of immune characteristics

We used ssGSEA (41) to estimate the immune cell proportion in two risk groups. The gene sets of 22 immune-related cells, innate immune response, and adaptive immune response were obtained from Charoentong et al. (42). The ESTIMATE algorithm (43) was employed to measure immune cell infiltration and stromal cell infiltration. The ESTIMATE score represents the combined score of immune score and stromal score. MCP-counter (44) was used to evaluate 10 immune-related cells including monocytic lineage, CD3+ T cells, NK cells, CD8+ T cells, endothelial cells, cytotoxic lymphocytes, B lymphocytes, neutrophils, fibroblasts, and myeloid dendritic cells based on the expression matrix. The immune checkpoint genes were downloaded from a previous study (42).

Assessment of biological pathways

Biological pathways were accessed from “h.all.v7.4.symbols.gmt” downloaded from the Molecular Signatures Database (MSigDB) (45). There were 13 tumor-related genes obtained from a previous research (46), which are classic cancer pathways, involved in the development and progression of cancer, including DNA damage repair (DDR), epithelial–mesenchymal transition (EMT), cell cycle, mismatch repair, CD8 T effector, FGFR3, nucleotide excision repair, base excision repair, DNA replication, homologous recombination, and WNT target. ssGSEA was performed to determine the enrichment score of biological pathways. The relation of risk score with pathways was inspected with Pearson correlation analysis.

Predicting the response to immunotherapy and chemotherapy

We conducted SubMap analysis (47) to compare the expression profiles between GSE102349 and IMvigor210. IMvigor210 contains the expression profiles of patients with metastatic urothelial carcinoma treated by PD-L1 inhibitors (48). The higher similarity of samples in GSE102349 with complete response (CR) and partial response (PR) groups in IMvigor210 suggested that the samples were more sensitive to anti-PD-L1 treatment. The TIDE (http://tide.dfci.harvard.edu/) algorithm (49) was employed to predict the escape and response to immune checkpoint inhibitors, according to the score of T-cell dysfunction and exclusion, the enrichment of immunosuppressive cells. The sensitivity of two risk groups to chemotherapeutic drugs was estimated using the pRRophetic R package (50).

The performance of the prognostic model was further examined in immunotherapy datasets including IMvigor210, GSE135222, and GSE78220. GSE135222 contains non-small cell lung cancer patients treated with immune checkpoint inhibitors (51). GSE78220 includes patients suffering from metastatic melanoma treated by anti-programmed cell death protein 1 (anti-PD-1) therapy.

Cell culture

After resuscitating NPC cells and normal cells, they were placed in 1,640 and 5A cell medium containing 10% fetal bovine serum, respectively, at 37°C, 5% CO2 concentration, and constant temperature to around 80%–90%, then passed and spread on a plate.

RT-qPCR

Total RNA was extracted from cells by RT-qPCR, and cDNA was synthesized by reverse transcription. After reverse transcription, samples were added according to the experimental instructions. The reaction conditions of RT-qPCR were 95°C for 30 s. 95°C, 5 s; 60°C, 30 s, 40 cycles. The mRNA expressions of miR-142-3p and miR-93 in the experimental group and control group were analyzed.

Statistical analysis

The statistical methods used in this study were performed in R software (v4.2.0). The Sangerbox platform (52) was used to provide an assistant in data analysis. The log-rank test was applied in survival analysis. Difference between two groups was examined by the Wilcoxon test. The Kruskal–Walls test was used to test the difference among over two groups. P < 0.05 was determined as statistically significant.

Results

Identification of DEmiRNAs related to NPC based on WGCNA

To begin with, we assessed the expression difference of miRNAs between normal and tumor samples in the GSE32960 dataset. The results showed that 332 miRNAs were differentially expressed between normal and tumor groups, including 168 upregulated and 164 downregulated miRNAs in tumor groups (Figure 1A; Table S1). Then, we applied WGCNA to cluster samples and dig out key gene modules based on DEmiRNAs. To ensure the co-expression network to be a scale-free network, the Pearson coefficient was selected as 0.85 and the power of the soft threshold was determined as 3 to construct an adjacency matrix (Figures 1B–D). Next, a topology overlap matrix (TOM) was generated based on the adjacency matrix and the genes were divided into different modules using the Dynamic Tree Cut algorithm (Figure 1E). Four modules were finally determined after merging the adjacent modules. Then, we analyzed the correlation of four modules with different sample groups. As a result, brown and blue modules were found to be evidently associated with sample groups and they manifested opposite correlations with two groups (Figure 1F). Specifically, blue and brown modules were negatively correlated with the tumor group (R = -0.35 and -0.53, respectively). The Venn plot was constructed to describe the intersection between DEmiRNAs and miRNAs of blue and brown modules (Figure 1G). There were 99 DEmiRNAs including 42 upregulated and 55 downregulated found in the blue module. There were 169 DEmiRNAs including 79 upregulated and 90 downregulated found in the brown module. The above total 268 DEmiRNAs were determined as potential miRNAs associated with NPC.

FIGURE 1
www.frontiersin.org

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.

Establishment and verification of a miRNA-related risk model

We randomly divided 312 NPC samples of the GSE32960 dataset at a ratio of 3:2 into the training and test groups (Table S2). We used the training group to screen prognostic miRNAs from 268 DEmiRNAs according to the univariate Cox regression model. The analysis identified two miRNAs (hsa-miR-142-3p and hsa-miR-93) that were significantly associated with the overall survival (P < 0.01). Based on the multivariate coefficients of two miRNAs by stepAIC analysis, we established the risk model defined as follows (Figure S1A).

Risk score = -0.835*(hsa-miR-142-3p) + 0.85*(hsa-miR-93)

Moreover, the RT-qPCR analysis results showed that miR-142-3p expression was downregulated and miR-93 was upregulated in five NPC cell lines in comparison with normal NP69 cells (Figures S1B, C).

To validate the performance of the risk model, we calculated the risk score for each tumor sample in the training and test groups. The AUC for 1-, 3-, and 5-year survival derived from the ROC curve analysis was 0.56, 0.70, and 0.70 in the training group and 0.96, 0.71, and 0.71 in the test group, respectively (Figures 2A, B). Furthermore, based on the median risk score, tumor samples were classified into low-risk and high-risk groups. The two risk groups had distinct overall survival (OS) in the training and test groups, as shown by Kaplan–Meier survival analysis (P = 0.00035 and P = 0.013, respectively, Figures 2A, B).

FIGURE 2
www.frontiersin.org

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.

In the total GSE32960 dataset, the risk model still showed good performance in predicting overall survival (Figure 2C). In addition, we evaluated the effectiveness of the risk model in different survival times including recurrence-free survival (RFS), disease-free survival (DFS), and metastasis-free survival (MFS). High-risk and low-risk groups exhibited differential prognosis of DFS and MFS (P < 0.0001) but showed no significant difference in RFS classification (P = 0.063, Figure 2C). In another independent dataset (GSE70970), the risk model also showed an effect classification for both OS and DFS (P = 0.017 and P = 0.022, respectively) (Figure 2D). Moreover, we verified the effectiveness of the risk model in the groups of different clinical characteristics. In addition to the female group and N0 group, the risk model was sufficient to distinguish samples into different risk levels in the groups of age >45, age ≤45, male, T1–T2, T3–T4, N1–N3, I–II, and III–IV (Figure 3).

FIGURE 3
www.frontiersin.org

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

Boosting the prediction efficiency of the miRNA-related risk model by constructing a nomogram

Using univariate and multivariate Cox regression analyses, we evaluated the relation of clinical features and risk type with prognosis. The result displayed that only gender and risk type were independent risk factors in both univariate and multivariate analyses (Figures 4A, B). Gender had a hazard ratio (HR) of 2.2 in both univariate and multivariate analyses (P = 0.019 and P = 0.02, respectively). Risk type had HR of 3.7 (P = 1.1e-6) and 3.6 (P = 2.1e-6) in univariate and multivariate analyses, respectively. Therefore, we included gender and risk score to construct the nomogram (Figure 4C). Compared with gender, risk score contributed the more total points in the nomogram. Calibration curve analysis suggested that the predicted 1-, 3-, and 5-year survival was highly accordant with the observed survival (Figure 4D). In addition, decision curve analysis was implemented to evaluate the benefit that patients may obtain from gender, risk score, and nomogram. As a result, nomogram performed a more favorable performance than gender and risk score (Figure 4E). Consequently, the nomogram based on gender and risk score was more efficient to predict the prognosis of NPC patients.

FIGURE 4
www.frontiersin.org

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.

Construction of miRNA–mRNA competing endogenous RNA networks

In the above sections, we identified two key miRNAs (hsa-miR-142-3p and hsa-miR-93) that had a close relation with NPC prognosis. To understand the potential molecular mechanism of the miRNAs, we applied 12 tools (starbase, PITA, TargetMiner, miRanda, microT, miRmap, miRtarbase, mircode, TargetScan, RNA22, miRDB, PicTar) to predict the potential targets of two miRNAs. The results output 39 target genes of hsa-miR-142-3p and 86 target genes of hsa-miR-93, which were visualized in ceRNA networks (Figure 5A). KEGG analysis demonstrated the involvement of these target genes in mitophagy and autophagy (Figure 5B). The top 10 enriched terms of cellular component and molecular function, as well as biological process, were visualized using GO function analysis (Figures 5C–E). For example, biological process terms of positive regulation of amyloid−beta metabolic process, amyloid−beta formation, and negative regulation of phosphatase activity were enriched (Figure 5C). Cellular component terms of clathrin-coated pit and trans-Golgi network membrane were enriched (Figure 5D). Molecular function terms of clathrin adaptor activity and clathrin heavy chain binding were enriched (Figure 5E).

FIGURE 5
www.frontiersin.org

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.

Establishing a mRNA prognostic model based on key target genes of hsa-miR-142-3p and hsa-miR-93

We predicted a total of 125 potential target genes of hsa-miR-142-3p and hsa-miR-93. Then, we used the univariate Cox regression model to analyze the relation between target genes and overall survival. Random sampling for 1,000 times from the samples of the GSE102349 dataset was performed. The target genes closely related to overall survival were remained (P < 0.05) and were ranked by the occurring frequency from 1,000-times analysis. The top five frequent target genes were selected as key target genes for establishing the mRNA prognostic model (Figure 6A). The coefficients of five target genes were calculated by multivariate analysis. Finally, the prognostic model was defined as the following: mRNA risk score = 1.333* E2F1 - 1.766*KCNJ8 + 1.075*SUCO - 1.030* HECTD1 - 0.340*KIF23.

FIGURE 6
www.frontiersin.org

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.

The risk score for each sample in the GSE102349 dataset was determined with the mRNA prognostic model. ROC curve analysis showed that the model was efficient in predicting 1-, 3-, and 5-year survival, with AUCs of 0.87, 0.81, and 0.79, respectively (Figure 6B). The median risk score value was used in classifying NPC samples into two groups, high-risk and low-risk groups. Kaplan–Meier survival curves of two risk groups showed that they had an apparently different prognosis (P = 0.0018, Figure 6B). Hence, the five target genes could be validated by the above results to be closely involved in the prognosis.

Immune characteristics of high- and low-risk groups

The tumor microenvironment plays a central role in antitumor response and immunotherapeutic response. We used several tools to assess the immune cell component, as well as immune and stromal infiltration, and analyzed immune checkpoint genes for their expression levels. Two risk groups showed a significantly different immune microenvironment. We estimated an ssGSEA enrichment score of 28 immune-related cells and found that 25 of them had a differential enrichment score, such as myeloid-derived suppressor cells (MDSCs), activated CD8 T cells, regulatory T cells, natural killer cells, activated B cells, and macrophages (Figure 7A). Most immune cells showed a higher abundance in the group with a low risk. In addition, ssGSEA also revealed higher enrichment of both adaptive and innate immune response scores in the low-risk group (P < 0.0001, Figure 7B). The ESTIMATE algorithm was used to evaluate stromal and immune infiltration of two groups, and not surprisingly, the low-risk group showed both higher immune score and stromal score than the high-risk group (Figure 7C, P < 0.0001). Moreover, MCP-counter was employed to dig out similar results with ssGSEA. A total of 10 types of immune-related cells all showed a higher enrichment score in the low-risk group (Figure 7D). Immune checkpoint expression levels also had an extreme difference in two risk groups that most of immune checkpoints, such as CTLA4, TIGIT, LAG3, and PCDC1, were more highly expressed in the low-risk group than in the high-risk group (Figure 7E). Immune checkpoints’ differential expression may contribute to different antitumor immune responses. The distinct immune microenvironment suggested that the five prognostic genes may play critical roles in immune modulation.

FIGURE 7
www.frontiersin.org

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.

Analysis of biological pathways, immunotherapeutic response, and drug sensitivity in two risk groups

To further understand the molecular mechanisms contributing different prognoses of two risk groups, the enrichment score of pathways from the “h.all.v7.4.symbols.gmt” file was calculated using ssGSEA. There were 30 pathways differentially activated in the two risk groups (Figure 8A). Immune-related pathways, for instance, interferon alpha response, interferon gamma response, inflammation response, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, and complement, were significantly more activated in the low-risk group, which was consistent with the result of immune analysis. Cell cycle-related pathways such as E2F targets, MYC target V2, MYC target V1, and G2M checkpoint were less enriched in the low-risk group. In addition, some oncogenic pathways were more activated in the high-risk group, such as Wnt signaling and Hedgehog signaling. Notably, apoptosis was more enriched in the low-risk group, which was associated with good prognosis. Risk score also manifested significant correlations with the above pathways (Figure S2).

FIGURE 8
www.frontiersin.org

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.

Gene sets of 13 tumor-related pathways were obtained from a previous study, and their enrichment scores were calculated using ssGSEA. Pearson correlation analysis uncovered a negative association of risk score with DNA repair-related pathways including DDR (R = -0.41), base excision repair (R = -0.25), nucleotide excision repair (R = -0.33), homologous recombination (R = -0.33), and mismatch repair (R = -0.33) (Figure 8B).

Next, we evaluated the response of two risk groups to chemotherapy and immunotherapy. The similarity of expression profiles between GSE102349 and IMvigor210 (treated by PD-L1 inhibitors) datasets was shown by SubMap analysis. Higher similarity between two datasets indicates higher sensitivity to PD-L1 inhibitors. The results presented that the low-risk group had a higher similarity with CR and PR groups (P < 0.05, Figure 8C), implying that the low-risk group could obtain more benefit from immunotherapy. Furthermore, the TIDE algorithm was implemented to predict immune escape to immune checkpoint inhibitors. Two risk groups showed no significant difference of TIDE score (Figure 8D). However, the low-risk group exhibited more severe T-cell dysfunction than the high-risk group in which enrichment of MDSCs and M2 tumor-associated macrophages (TAM) was higher, contributing to its higher T-cell exclusion. In the response to chemotherapeutic drugs, we screened a total of 103 drugs including 46 drugs such as YM155 and vinorelbine sensitive to high-risk groups and 57 drugs such as sunitinib and temsirolimus sensitive to the low-risk group (Figure 8E).

The analysis of endocrine metabolism

Abnormal endocrine metabolism is one of the complications of tumor treatment. As shown in Figure 9A, endocrine-related genes, such as MON1B, SCAMP2, and FAM20A, were negatively associated with RiskScore; SCAMP3 was positively correlated with RiskScore. Pearson analysis between endocrine pathways and RiskScore showed that most endocrine pathways were negatively related to RiskScore (Figure 9B).

FIGURE 9
www.frontiersin.org

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

The performance of the mRNA prognostic model in immunotherapy datasets

To further examine the exhibition of the five-gene prognostic model, we introduced three independent datasets containing patients administrated by immunotherapy (IMvigor210, GSE135222, and GSE78220). Using the same calculation method, we measured the risk score of each patient in three datasets. In the IMvigor210 dataset, the high-risk group displayed apparently worse prognosis than the low-risk group (P < 0.0001, Figure 10A). The prediction efficiency was evaluated by ROC, with AUCs of 0.61, 0.66, and 0.64 at 0.5-, 1-, and 1.5-year survival, respectively (Figure 10A). In addition, we analyzed the proportion of different response groups in two risk groups. The low-risk group had a higher proportion of CR and PR groups (11% and 20%, respectively), compared with the high-risk group (6% and 9%, respectively). The CR and PR groups also showed a lower risk score than PD and SD groups. In the GSE135222 and GSE78220 datasets, we observed correspondent results (Figures 10B, C). The GSE135222 dataset showed AUCs of 0.82 and 0.85 at 0.5- and 1-year survival, respectively. The low-risk group exhibited a markedly higher percentage of CR/PR patients compared with the high-risk group (50% versus 8% in low- versus high-risk groups). The GSE78220 dataset showed AUCs of 0.74 and 0.71 at 1- and 2-year survival, respectively. Expectedly, the CR and PR groups were more accumulated in the low-risk group (21% and 43%) compared with the high-risk group (8% and 31%). Moreover, the CR/PR group showed a lower risk score than the PD/SD group in both GSE135222 and GSE78220 datasets, but the difference was not significant in the GSE78220 dataset. The above observation suggested that patients with a low risk score were more sensitive to immunotherapy and could attain longer survival.

FIGURE 10
www.frontiersin.org

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.

Discussion

MiRNAs are considered as potential therapeutic targets for NPC treatment, and till now abundant miRNAs have been unveiled to be dysregulated in NPC patients (19). MiRNAs serve as importantly regulatory roles in gene expression and pathway modulation. In this study, we dug out the potential key miRNAs and mRNAs that were probably responsible for NPC development and progression. We interpreted the association between key miRNAs and mRNAs and decoded the relationships of prognostic miRNA-related mRNAs with tumor microenvironment, immunotherapy, and functional pathways.

To begin with, we deciphered the miRNA expression data and screened 332 aberrantly expressed miRNAs in NPC samples compared with normal samples. By using WGCNA, we further identified two key miRNAs (hsa-miR-142-3p and hsa-miR-93) strongly related to NPC prognosis and phenotype. We constructed a miRNA risk model based on hsa-miR-142-3p and hsa-miR-93. The risk model exhibited an intensive relation with NPC prognosis in two independent datasets (GSE32960 and GSE70970). Patients with a high risk showed significantly worse OS and DFS than the low-risk group. In addition, the risk model was also effective to predict OS in NPC samples with different clinical characteristics including ages, T stage, N1–N3 stage, and AJCC stage I–IV. These results demonstrated that hsa-miR-142-3p and hsa-miR-93 were highly responsible for NPC development.

Moreover, RiskScore was negatively correlated with endocrine genes, especially FAM20A, which had been important for endocrine-related tumors (53). The FAM20 family of kinases is a newly discovered class of secreted kinases that are capable of phosphorylating secreted proteins and proteoglycans (54). FAM20A may play a more complex role in gliomas, as correlations between FAM20A genes and low-grade gliomas have been found (55). Combining the known literature and the results of this study, it is suggested that endocrine gene FAM20A may be closely related to NPC.

The two key miRNAs have been reported in the contribution of other cancer types. For example, in esophageal squamous cell carcinoma (ESCA), hsa-miR-142-3p was identified as a prognostic biomarker (56). Non-small cell lung cancer (NSCLC) cells could be promoted by overexpression of miR-142-3p via interfering TGFβR1 expression (57). However, Dong et al. revealed that miR-142-3p suppressed the growth of human cervical cancer cells by attenuating HMGB1 expression levels (58). Moreover, Sharma et al. excavated that miR-142-3p functioned as a tumor-suppressive miRNA through modulating the expression of HMGA1, A2, B1, and B3 in human cervical cancer (59). miR-142-3p in different cancer types showed a discord in expression that miR-142-3p was suggested to play complicated roles (both promotive and suppressive) by interacting with specific pathways and genes in different cancers. In our study, the miR-142-3p level was significantly decreased in NPC samples, implying that high-expressed miR-142-3p may facilitate the progression of NPC.

The role of hsa-miR-93 has also been uncovered in different cancer types. For instance, a miRNA microarray result showed that miR-93 was downregulated in human colon cancer stem cells and overexpressing miR-93 strikingly inhibited cell proliferation and colony formation (60). In triple-negative breast cancer cells, cell migratory capability and invasive potential could be weakened by overexpression of mature miR-93-5p possibly by targeting WNK1 (61). In uterine cancer, the high miRNA-93 expression group had an evidently higher survival rate than the low miRNA-93 expression group (62). The results of the above studies were accordant with our study that miRNA-93 expression was elevated in the NPC group compared with the normal group.

To clarify the potential mechanisms of the two miRNAs in NPC, their potential target genes (mRNAs) were predicted by utilizing 12 different tools to build ceRNA networks. As a result, we confirmed 39 target genes of miR-142-3p and 86 target genes of miR-93. KEGG analysis revealed that these target genes were significantly involved in autophagy and mitophagy. Autophagy occurs under stressful situations such as the presence of abnormal proteins and nutrient deprivation, which degrades cellular proteins and organelles to provide precursors for recycling (63). Some clinical trials have presented that inhibiting autophagy has feasible benefits in multiple cancer types such as glioblastoma, melanoma, and pancreatic cancer (64). Autophagy and mitophagy are demonstrated to contribute to the reprogramming of cancer metabolism that is a major challenge for anticancer therapy (65). Therefore, we supposed that maybe one of the mechanisms of miR-142-3p and miR-93 in NPC development was their participation in autophagy and mitophagy process responsible for cancer metabolism.

In order to distinguish key target genes of the two miRNAs, we applied random sampling and univariate Cox regression on 125 potential target genes. As a consequence, we confirmed five target genes that had prognostic effects on NPC, namely, E2F1, KCNJ8, SUCO, HECTD1, and KIF23, where SUCO and HECTD1 are the targets of miR-142-3p and E2F1, KCNJ8, and KIF23 are the targets of miR-93. E2F1 has been widely reported to regulate cell cycle and cell death and has a significant role in multiple cancer types. E2F1 target pathways are considered as important targets for cancer treatment (66). Limited studies of cancer have been found related to the other four genes. Based on these five target genes, we further established a prognostic model. The five-gene prognostic model manifested substantial performance in predicting 1-, 2-, and 3-year survival with AUCs of 0.87, 0.81, and 0.79 respectively. According to the model, high-risk and low-risk groups were defined with disparate prognosis.

We further investigated the differences of the tumor microenvironment and functional pathways between two risk groups for excavating the biological influence of the five target genes in NPC. Two risk groups had distinct immune cell infiltration. Anticancer immune cells, for instance, activated B cells, NK cells, dendritic cells, and CD8 T cells, were evidently higher in the low-risk group compared with the high-risk group, which led to the stronger anticancer response and clearance of cancer cells. The high-risk group showed both higher innate and adaptive immune response than the low-risk group. Analysis of biological pathways unveiled that the low-risk group displayed higher activation of immune-related pathways, for instance, IL2-STAT5 signaling, interferon alpha response, interferon gamma response, IL6-JAK-STAT3 signaling, inflammation response, and complement. Importantly, DNA repair pathways, cell cycle-related pathways, and apoptosis were also more enriched in the low-risk group, supporting that the two key miRNAs may have an interaction with these target genes. However, the specific association between the miRNAs and five target genes should be further validated in future experiments.

The different proportion of tumor-infiltrating immune cells has a profound effect on both cancer prognosis and immunotherapy (67). SubMap analysis predicted that the low-risk group was more sensitive to ICB therapy than the high-risk group, which may result from the higher expression of critical immune checkpoints such as CTLA4, IDO1, LAG3, and PDCD1 (PD-1) in the low-risk group. Anti-PD-1/PD-L1 therapy has shown some favorable outcomes in clinical trials of NPC (68). Our five-gene model can help clinicians to better select the patients sensitive to ICB therapy and raise the efficiency of immunotherapy.

Conclusions

In conclusion, this study identified two key miRNAs (miR-142-3p and miR-93) and predicted their potential key target genes. MiR-142-3p and miR-93 contributed to NPC survival possibly through regulating autophagy pathways. In addition, we confirmed five prognostic target genes (E2F1, KCNJ8, SUCO, HECTD1, and KIF23) and constructed a five-gene prognostic model. The model was effective to predicting NPC prognosis and could provide a guidance for personalized immunotherapy in NPC patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

All authors contributed to this present work. XXZ, XiaL and CW designed the study, SW acquired the data. YZ and BL drafted the manuscript. XinL revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by Key Projects of Intergovernmental Cooperation in National Key R&D Programs (Grant no. 2022YFE0131800), Natural Science Foundation of Liaoning Province (Grant no. 2022-MS-235), and Project of Shenyang Science and Technology Bureau (Grant no. RC210316).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1174911/full#supplementary-material

Abbreviations

ASIR, age-standardized incidence rate; ceRNA, competing endogenous RNA; CR, complete response; DEmiRNAs, differentially expressed miRNAs; DFS, disease-free survival; DDR, DNA damage repair; EBV, Epstein–Barr virus; GEO, Gene Expression Omnibus; GO, Gene Ontology; HR, hazard ratio; ICB, immune checkpoint blockade; KEGG, Kyoto Encyclopedia of Genes and Genomes; MFS, metastasis-free survival; miRNAs, micro RNAs; MDSCs, myeloid-derived suppressor cells; NPC, nasopharyngeal cancer; PR, partial response; PD-L1, programmed death ligand 1; PD-1, programmed death protein 1; ROC, receiver operating characteristic; ssGSEA, single sample gene set enrichment analysis; TOM, topological overlap matrix; WGCNA, weighted correlation network analysis.

References

1. Chan KCA, Woo JKS, King A, Zee BCY, Lam WKJ, Chan SL, et al. Analysis of plasma Epstein-Barr virus DNA to screen for nasopharyngeal cancer. New Engl J Med (2017) 377(6):513–22. doi: 10.1056/NEJMoa1701717

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Jain A, Chia WK, Toh HC. Immunotherapy for nasopharyngeal cancer-a review. Chin Clin Oncol (2016) 5(2):22. doi: 10.21037/cco.2016.03.08

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Chen YP, Chan ATC, Le QT, Blanchard P, Sun Y, Ma J. Nasopharyngeal carcinoma. Lancet (London England). (2019) 394(10192):64–80. doi: 10.1016/S0140-6736(19)30956-0

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Bai R, Sun J, Xu Y, Sun Z, Zhao X. Incidence and mortality trends of nasopharynx cancer from 1990 to 2019 in China: an age-period-cohort analysis. BMC Public Health (2022) 22(1):1351. doi: 10.1186/s12889-022-13688-7

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Bhattacharyya T, Babu G, Kainickal CT. Current role of chemotherapy in nonmetastatic nasopharyngeal cancer. J Oncol (2018) 2018:3725837. doi: 10.1155/2018/3725837

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Verma V, Allen PK, Simone CB 2nd, Gay HA, Lin SH. Addition of definitive radiotherapy to chemotherapy in patients with newly diagnosed metastatic nasopharyngeal cancer. J Natl Compr Cancer Network JNCCN. (2017) 15(11):1383–91. doi: 10.6004/jnccn.2017.7001

CrossRef Full Text | Google Scholar

7. Tseng M, Ho F, Leong YH, Wong LC, Tham IW, Cheo T, et al. Emerging radiotherapy technologies and trends in nasopharyngeal cancer. Cancer Commun (London England). (2020) 40(9):395–405. doi: 10.1002/cac2.12082

CrossRef Full Text | Google Scholar

8. Zhang J, Fang W, Qin T, Yang Y, Hong S, Liang W, et al. Co-Expression of PD-1 and PD-L1 predicts poor outcome in nasopharyngeal carcinoma. Med Oncol (Northwood London England). (2015) 32(3):86. doi: 10.1007/s12032-015-0501-6

CrossRef Full Text | Google Scholar

9. Zhou Y, Miao J, Wu H, Tang H, Kuang J, Zhou X, et al. PD-1 and PD-L1 expression in 132 recurrent nasopharyngeal carcinoma: the correlation with anemia and outcomes. Oncotarget (2017) 8(31):51210–23. doi: 10.18632/oncotarget.17214

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Hsu C, Lee SH, Ejadi S, Even C, Cohen RB, Le Tourneau C, et al. Safety and antitumor activity of pembrolizumab in patients with programmed death-ligand 1-positive nasopharyngeal carcinoma: results of the KEYNOTE-028 study. J Clin Oncol (2017) 35(36):4050–6. doi: 10.1200/JCO.2017.73.3675

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ma BBY, Lim WT, Goh BC, Hui EP, Lo KW, Pettinger A, et al. Antitumor activity of nivolumab in recurrent and metastatic nasopharyngeal carcinoma: an international, multicenter study of the Mayo clinic phase 2 consortium (NCI-9742). J Clin Oncol (2018) 36(14):1412–8. doi: 10.1200/JCO.2017.77.0388

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Chen G, Huang AC, Zhang W, Zhang G, Wu M, Xu W, et al. Exosomal PD-L1 contributes to immunosuppression and is associated with anti-PD-1 response. Nature (2018) 560(7718):382–6. doi: 10.1038/s41586-018-0392-8

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Garcia-Diaz A, Shin DS, Moreno BH, Saco J, Escuin-Ordinas H, Rodriguez GA, et al. Interferon receptor signaling pathways regulating PD-L1 and PD-L2 expression. Cell Rep (2017) 19(6):1189–201. doi: 10.1016/j.celrep.2017.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Holly JM, Perks CM. Cancer as an endocrine problem. Best Pract Res Clin Endocrinol Metab (2008) 22(4):539–50. doi: 10.1016/j.beem.2008.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Cascio S, Bartella V, Auriemma A, Johannes GJ, Russo A, Giordano A, et al. Mechanism of leptin expression in breast cancer cells: role of hypoxia-inducible factor-1alpha. Oncogene (2008) 27(4):540–7. doi: 10.1038/sj.onc.1210660

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Schally AV, Varga JL, Engel JB. Antagonists of growth-hormone-releasing hormone: an emerging new therapy for cancer. Nat Clin Pract Endocrinol Metab (2008) 4(1):33–43. doi: 10.1038/ncpendmet0677

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Schally AV. New approaches to the therapy of various tumors based on peptide analogues. Hormone Metab Res = Hormon- und Stoffwechselforschung = Hormones metabolisme. (2008) 40(5):315–22. doi: 10.1055/s-2008-1073142

CrossRef Full Text | Google Scholar

18. Eichmüller SB, Osen W, Mandelboim O, Seliger B. Immune modulatory microRNAs involved in tumor attack and tumor immune escape. J Natl Cancer Institute (2017) 109(10). doi: 10.1093/jnci/djx034

CrossRef Full Text | Google Scholar

19. Wang S, Claret FX, Wu W. MicroRNAs as therapeutic targets in nasopharyngeal carcinoma. Front Oncol (2019) 9:756. doi: 10.3389/fonc.2019.00756

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yu L, Lu J, Zhang B, Liu X, Wang L, Li SY, et al. miR-26a inhibits invasion and metastasis of nasopharyngeal cancer by targeting EZH2. Oncol Letters. (2013) 5(4):1223–8. doi: 10.3892/ol.2013.1173

CrossRef Full Text | Google Scholar

21. Yi C, Wang Q, Wang L, Huang Y, Li L, Liu L, et al. MiR-663, a microRNA targeting p21(WAF1/CIP1), promotes the proliferation and tumorigenesis of nasopharyngeal carcinoma. Oncogene (2012) 31(41):4421–33. doi: 10.1038/onc.2011.629

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bell E, Taylor MA. Functional roles for exosomal MicroRNAs in the tumour microenvironment. Comput Struct Biotechnol J (2017) 15:8–13. doi: 10.1016/j.csbj.2016.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Liao C, Liu H, Luo X. The emerging roles of exosomal miRNAs in nasopharyngeal carcinoma. Am J Cancer Res (2021) 11(6):2508–20.

PubMed Abstract | Google Scholar

24. Clough E, Barrett T. The gene expression omnibus database. Methods Mol Biol (Clifton NJ). (2016) 1418:93–110. doi: 10.1007/978-1-4939-3578-9_5

CrossRef Full Text | Google Scholar

25. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Langfelder P, Horvath S. WGCNA: an r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

27. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med (2013) 32(30):5381–97. doi: 10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Paraskevopoulou MD, Georgakilas G, Kostoulas N, Vlachos IS, Vergoulis T, Reczko M, et al. DIANA-microT web server v5.0: service integration into miRNA functional analysis workflows. Nucleic Acids Res (2013) 41(Web Server issue):W169–73. doi: 10.1093/nar/gkt393

PubMed Abstract | CrossRef Full Text | Google Scholar

29. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets. PloS Biol (2004) 2(11):e363. doi: 10.1371/journal.pbio.0020363

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Jeggari A, Marks DS, Larsson E. miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinf (Oxford England). (2012) 28(15):2062–3. doi: 10.1093/bioinformatics/bts344

CrossRef Full Text | Google Scholar

31. Chen Y, Wang X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res (2020) 48(D1):D127–d31. doi: 10.1093/nar/gkz757

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Vejnar CE, Blum M, Zdobnov EM. miRmap web: comprehensive microRNA target prediction online. Nucleic Acids Res (2013) 41(Web Server issue):W165–8. doi: 10.1093/nar/gkt430

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Huang HY, Lin YC, Li J, Huang KY, Shrestha S, Hong HC, et al. miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res (2020) 48(D1):D148–d54. doi: 10.1093/nar/gkz896

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Krek A, Grün D, Poy MN, Wolf R, Rosenberg L, Epstein EJ, et al. Combinatorial microRNA target predictions. Nat Genet (2005) 37(5):495–500. doi: 10.1038/ng1536

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. The role of site accessibility in microRNA target recognition. Nat Genet (2007) 39(10):1278–84. doi: 10.1038/ng2135

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Bandyopadhyay S, Mitra R. TargetMiner: microRNA target prediction with systematic identification of tissue-specific negative examples. Bioinf (Oxford England). (2009) 25(20):2625–31. doi: 10.1093/bioinformatics/btp503

CrossRef Full Text | Google Scholar

37. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife (2015) 4:e05005. doi: 10.7554/eLife.05005

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Loher P, Rigoutsos I. Interactive exploration of RNA22 microRNA target predictions. Bioinf (Oxford England). (2012) 28(24):3322–3. doi: 10.1093/bioinformatics/bts615

CrossRef Full Text | Google Scholar

39. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-seq data. Nucleic Acids Res (2014) 42(Database issue):D92–7. doi: 10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

40. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

42. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17(1):218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Systems. (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Roh W, Chen PL, Reuben A, Spencer CN, Prieto PA, Miller JP, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Trans Med (2017) 9(379):eaah3560. doi: 10.1126/scitranslmed.aah3560

CrossRef Full Text | Google Scholar

48. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, et al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet (London England). (2017) 389(10064):67–76. doi: 10.1016/S0140-6736(16)32455-2

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Kim JY, Choi JK, Jung H. Genome-wide methylation patterns predict clinical benefit of immunotherapy in lung cancer. Clin Epigenetics. (2020) 12(1):119. doi: 10.1186/s13148-020-00907-4

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Shen W, Song Z, Xiao Z, Huang M, Shen D, Gao P, et al. Sangerbox: a comprehensive, interaction-friendly clinical bioinformatics analysis platform. iMeta (2022) 1(3):e36. doi: 10.1002/imt2.36

CrossRef Full Text | Google Scholar

53. Wells SA Jr. Progress in endocrine neoplasia. Clin Cancer Res (2016) 22(20):4981–8. doi: 10.1158/1078-0432.CCR-16-0384

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zhang H, Zhu Q, Cui J, Wang Y, Chen MJ, Guo X, et al. Structure and evolution of the Fam20 kinases. Nat Commun (2018) 9(1):1218. doi: 10.1038/s41467-018-03615-z

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Feng J, Zhou J, Zhao L, Wang X, Ma D, Xu B, et al. Fam20C overexpression predicts poor outcomes and is a diagnostic biomarker in lower-grade glioma. Front Genet (2021) 12:757014. doi: 10.3389/fgene.2021.757014

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Lin RJ, Xiao DW, Liao LD, Chen T, Xie ZF, Huang WZ, et al. MiR-142-3p as a potential prognostic biomarker for esophageal squamous cell carcinoma. J Surg Oncol (2012) 105(2):175–82. doi: 10.1002/jso.22066

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lei Z, Xu G, Wang L, Yang H, Liu X, Zhao J, et al. MiR-142-3p represses TGF-β-induced growth inhibition through repression of TGFβR1 in non-small cell lung cancer. FASEB J (2014) 28(6):2696–704. doi: 10.1096/fj.13-247288

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Dong H, Song J. miR-142-3p reduces the viability of human cervical cancer cells by negatively regulating the cytoplasmic localization of HMGB1. Exp Ther Med (2021) 21(3):212. doi: 10.3892/etm.2021.9644

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Sharma P, Yadav P, Jain RP, Bera AK, Karunagaran D. miR-142-3p simultaneously targets HMGA1, HMGA2, HMGB1, and HMGB3 and inhibits tumorigenic properties and in-vivo metastatic potential of human cervical cancer cells. Life Sci (2022) 291:120268. doi: 10.1016/j.lfs.2021.120268

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Yu XF, Zou J, Bao ZJ, Dong J. miR-93 suppresses proliferation and colony formation of human colon cancer stem cells. World J Gastroenterology. (2011) 17(42):4711–7. doi: 10.3748/wjg.v17.i42.4711

CrossRef Full Text | Google Scholar

61. Shyamasundar S, Lim JP, Bay BH. miR-93 inhibits the invasive potential of triple-negative breast cancer cells in vitro via protein kinase WNK1. Int J Oncol (2016) 49(6):2629–36. doi: 10.3892/ijo.2016.3761

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Fang S, Gao M, Xiong S, Chen Q, Zhang H. Expression of serum hsa-miR-93 in uterine cancer and its clinical significance. Oncol Letters. (2018) 15(6):9896–900. doi: 10.3892/ol.2018.8553

CrossRef Full Text | Google Scholar

63. Yun CW, Lee SH. The roles of autophagy in cancer. Int J Mol Sci (2018) 19(11):3466. doi: 10.3390/ijms19113466

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Levy JMM, Towers CG, Thorburn A. Targeting autophagy in cancer. Nat Rev Cancer. (2017) 17(9):528–42. doi: 10.1038/nrc.2017.53

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Ferro F, Servais S, Besson P, Roger S, Dumas JF, Brisson L. Autophagy and mitophagy in cancer metabolic remodelling. Semin Cell Dev Biol (2020) 98:129–38. doi: 10.1016/j.semcdb.2019.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Huang Y, Chen R, Zhou J. E2F1 and NF-κB: key mediators of inflammation-associated cancers and potential therapeutic targets. Curr Cancer Drug Targets. (2016) 16(9):765–72. doi: 10.2174/1568009616666160216130755

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Lu J, Chen XM, Huang HR, Zhao FP, Wang F, Liu X, et al. Detailed analysis of inflammatory cell infiltration and the prognostic impact on nasopharyngeal carcinoma. Head Neck. (2018) 40(6):1245–53. doi: 10.1002/hed.25104

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Adkins DR, Haddad RI. Clinical trial data of anti-PD-1/PD-L1 therapy for recurrent or metastatic nasopharyngeal carcinoma: a review. Cancer Treat Rev (2022) 109:102428. doi: 10.1016/j.ctrv.2022.102428

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Citation: Zhang X, Li X, Wang C, Wang S, Zhuang Y, Liu B and Lian X (2023) Identification of markers for predicting prognosis and endocrine metabolism in nasopharyngeal carcinoma by miRNA–mRNA network mining and machine learning. Front. Endocrinol. 14:1174911. doi: 10.3389/fendo.2023.1174911

Received: 27 February 2023; Accepted: 21 June 2023;
Published: 19 July 2023.

Edited by:

Wenjie Shi, Otto von Guericke University Magdeburg, Germany

Reviewed by:

Guangying Cui, Zhengzhou University, China
Liang Yu, Second Affiliated Hospital of Harbin Medical University, China

Copyright © 2023 Zhang, Li, Wang, Wang, Zhuang, Liu and Lian. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xin Lian, amVzc2ljYXNoZW5namluZ0Bob3RtYWlsLmNvbQ==

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.