Differences in HNRNPR mRNA expression in pan-cancer
We analyzed the transcription level of HNRNPR in various pan-cancer and normal individuals using the Oncomine online database and TCGA dataset. Data from the Oncomine database with fold change > 1. 5 and a p-value < 0. 001 were analyzed. Results showed that HNRNPR expression level was higher level in many cancers than in normal tissue (Fig. 1A).
Further analysis revealed that the mRNA level of HNRNPR was increased in many cancers, including ESCA (Fig. 1B).
HNRNPR expression predicted poor survival in ESCA patients
The analysis of ESCA patient data from the TCGA and GEO database revealed that HNRNPR expression was significantly elevated in cancer tissue than in normal tissue or para-cancerous tissue (Fig. 1C–E).
Next, IHC and qPCR experiments were conducted to confirm the expression and prediction accuracy of HNRNPR in ESCA. The results indicated that HNRNPR mRNA and protein level was significantly higher in para-cancerous tissues than in normal tissues, particularly in the IHC test (Fig. 1F–H).
Additionally, Kaplan–Meier survival analysis suggested that patients whose tumors exhibited higher expression of HNRNPR had a decreased overall survival (Fig. 2A). The ROC curve analysis revealed that HNRNPR had high predictive accuracy, with a ROC curve of 0. 897 (95% CI:0. 801–0. 993) (Fig. 2B).
Additionally, we examined the significant clinicopathological role of HNRNPR in ESCA using the TCGA database, which included histological type, weight, Body mass index (BMI), reflux history, smoking, and alcohol intake. The results indicated that HNRNPR expression was significantly higher in ESCC than in EAC (Fig. 2C). HNRNPR mRNA expression was decreased in people with higher body weight and BMI (Fig. 2D, E). HNRNPR also affected whether the patient developed reflux or not (Fig. 2F). There was, however, no significant change in HNRNPR expression between patients who smoked and those who did not (Fig. 2G, H).
Relationship between HNRNPR expression and PET metabolic parameters
By comparing the PET parameter values of patients to corresponding IHC results, we discovered that FDG uptake rate also had an effect on HNRNPR expression (Fig. 3A, B). Notably, HNRNPR was expressed only in the cytoplasm or in both the nucleus and cytoplasm, which may be related to the state of the cells at the time of expression and the functions activated by HNRNPR [19, 26].
To further investigate the relationship between 18F-FDG PET/CT metabolic parameters and HNRNPR expression in patients with ESCA, the IHC staining intensity of HNRNPR in tumor samples was determined, as well as the correlation with SUVmax, SUVmean, TLG, and MTV. As shown in the figure (Fig. 4A–D), the HNRNPR score was positively correlated with SUVmax, SUVmean, and TLG (rho = 0. 369, 0. 411, 0. 379, respectively, p < 0. 05), but not with MTV. We hypothesized that SUVmax, SUVmean, and TLG may be used to predict HNRNPR expression levels in ESCA based on the above data.
Logistic analysis was used to determine the significance of each significant parameter in the HNRNPR expression. The area under the SUVmax, SUVmean, and TLG ROC curves was 0. 730, 0. 733, 0. 703; and the combined index of the ROC curves for these three indexes was 0. 753, implying that the accuracy of the answer obtained by considering all three parameters comprehensively is more accurate (Fig. 4E, F).
HNRNPR Co-expression networks indicate the potential function of HNRNPR in ESCA
There were 20,140 HNRNPR-related co-expressed genes in TCGA ESCA, and we performed gene enrichment analysis using the LinkedOmics database. As shown in Fig. 5A, 11827 genes were positively correlated with HNRNPR (red dot) while 8313 genes were negatively correlated (green dot). KDM1A (r = 0.720, p = 1.04E−30), KHDRBS1(r = 0.687, p = 4.88E−27), CENPF (r = 0.680, p = 2.30E−26), C1orf135(r = 0.674, p = 1.07E-25), DEPDC1(r = 0.673, p = 1.19E−25) had the highest positive correlation with HNRNPR. Notably, MGLL (r = − 0.571, p = 2.71E−17) and SDCBP2 (r = − 0.548, p = 8.36E−16) showed the strongest negative correlation with HNRNPR. A heat map showing the top 50 significant genes that were positively and negatively correlated with HNRNPR expression is shown in Fig. 5B. Among them, KDM1A and KHDRBS1 were previously found to be significantly related to the occurrence and prognosis of ESCA [27, 28].
The GO function and KEGG pathway enrichment analysis of the first 400 co-expressed genes positively correlated with HNRNPR expression were analyzed in line with the p < 0.05 and q < 0.25 as cut-off values. HNRNPR co-expressed genes were implicated in 18 KEGG and multiple GO functions, and the GO functions involved were classified as follows: 549 biological processes (GO-BP), 102 cell components (GO-CC), and 70 molecular functions (GO-MF). The bubble graph demonstrates the top 5 messages, respectively (Fig. 5C–F).
Functions of HNRNPR and its related genes
20 HNRNPR-related genes were analyzed in terms of function by the GeneMANIA database. They surrounded HNRNPR and formed the GGI network. Genes in the GGI network were intertwined to form many "intersection points", which represented physical interactions (77.64%), co-expression (8.01%), predicted (5.37%), co-localization (3.63%), genetic interactions (2.87%), pathway (1.88%), and Shared protein domains (0.60%). Five genes, including SRSF9 (serine and arginine-rich splicing factor 9), U2AF2 (U2 small nuclear RNA auxiliary factor 2), AGTPBP1 (ATP/GTP binding protein 1), HNRNPA1 (heterogeneous nuclear ribonucleoprotein A1), and SYNCRIP (synaptotagmin binding cytoplasmic RNA interacting protein), had the strongest correlation with HNRNPR. The strongest interaction among the 21 genes was physical contact. Most of these 21 genes were related to splicing, metabolism, and, mRNA transport (Fig. 6A).
STRING was used to further analyze the PPI network of HNRNPR. The genes associated with HNRNPR were HNRNPA1, HNRNPM, HNRNPL, HNRNPH1, HNRNPK, HNRNPC, HNRNPA2B1, PTBP1, RBMX, and ILF2. Their combined scores are shown in Fig. 6B.
We deduced from gene enrichment analysis that HNRNPR participates in and regulates the cell cycle in a variety of pathways (Fig. 6C–F). Figure 6G–J shows the list of several other pathways in which HNRNPR is involved, including the Fanconi pathway, which is a risk factor for ESCA.
Correlation between HNRNPR Expression and m6A modification in ESCA
m6A methylation is a critical RNA modification in eukaryotic cells, affecting numerous aspects of RNA metabolism. It not only contributes to the occurrence and development of cancer but also provides new ideas for its early diagnosis and treatment. The PI3K–AKT–mTOR, RAS–MAPK, and/or MYC signaling pathways are the main signal pathways of m6A [23].
According to the analysis of the TCGA dataset, 20 m6A-related genes were positively correlated with the expression of HNRNPR (P < 0.05), whereas in the GEO dataset, the HNRNPR expression is positively correlated with Methyltransferase 3 (METTL3), Methyltransferase 14 (METTL14), WT1 Associated Protein (WTAP), Vir Like m6A Methyltransferase Associated (VIRMA), RNA Binding Motif Protein 15 (RBM15), YTH Domain Containing 1 (YTHDC1), YTH N6-Methyladenosine RNA Binding Protein 1 (YTHDF1), YTH N6-Methyladenosine RNA Binding Protein 2 (YTHDF2), HNRNPC, RBMX, HNRNPA2B1; and negatively correlated with RBM15B (Fig. 7A, B). As shown in Fig. 7C, HNRNPR is involved in a number of the pathways outlined previously. There is differential expression in many m6A-related genes between high and low groups with HNRNPR expression in ESCA (Fig. 7D).
Correlations between HNRNPR expression and glycolysis in ESCA
The correlation between HNRNPR and glycolysis-related genes was analyzed using the TCGA and GEO databases. As shown in Fig. 8A, HNRNPR had a strong correlation with ADH5 (Alcohol Dehydrogenase 5), ALDH5A1 (Aldehyde Dehydrogenase 5 Family Member A1), DLD (Dihydrolipoamide Dehydrogenase), ENO2 (Enolase 2), G6PD (Glucose-6-Phosphate Dehydrogenase), GAPDH (Glyceraldehyde-3-Phosphate Dehydrogenase), HK1 (Hexokinase 1), LDHB (Lactate Dehydrogenase B), RBMX (RNA Binding Motif Protein X-Linked), SLC2A1 (Solute Carrier Family 2 Member 1), SLC2A3 (Solute Carrier Family 2 Member 3), and VIRMA. Additionally, except for the GSE69925 data set, G6PD, HK1, and SLC2A1 exhibited a negative correlation, whereas the rest showed a positive correlation.
Fig. 8B–G shows that HNRNPR was involved in many glycolysis-related pathways.