Xixia Zhang1
Xiao Li
Xin Lian- 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 (5–7). 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 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 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 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 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 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 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 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 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 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 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
53. Wells SA Jr. Progress in endocrine neoplasia. Clin Cancer Res (2016) 22(20):4981–8. doi: 10.1158/1078-0432.CCR-16-0384
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
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
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
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
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
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
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
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
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
63. Yun CW, Lee SH. The roles of autophagy in cancer. Int J Mol Sci (2018) 19(11):3466. doi: 10.3390/ijms19113466
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
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
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
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
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, GermanyReviewed by:
Guangying Cui, Zhengzhou University, ChinaLiang 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, jessicashengjing@hotmail.com