Expression of HNRNPR in Dataset
We analyzed the expression of HNRNPR in different tumors in the Oncomine (www. oncomine. org) online database, which allowed us to compare the discrepancies in HNRNPR expression levels between cancer and normal groups using the Student’s test, and we obtained the expression of HNRNPR in various pan-cancers and ESCA from the TCAG database. The difference in HNRNPR expression between esophageal squamous cell carcinoma and normal esophageal epithelium was studied and compared using the GSE45670 and GSE20347 data set, between tumor and normal adjacent tissues of patients with ESCC.
We analyzed the ESCA data sets in the TCAG database to compare the relationship between HNRNPR expression levels and clinical manifestations in ESCA patients, and then constructed a receiver operating characteristic (ROC) curve to determine the diagnostic prediction accuracy of HNRNPR in ESCA patients. The Kaplan–Meier Plotter (http://kmplot. com/analysis/) was used to evaluate the effect of gene expression on cancer patient survival.
Immunohistochemistry (IHC)
IHC staining was used to stain paraffin-embedded 4-µm-thick tissue sections. After dewaxing, fixed sections were microwaved for 10 min in Tris-EDTA (pH 9. 0) (ab93684; Abcam, USA). Normal serum was sealed at room temperature for 30 min. The sections were incubated with HNRNPR rabbit polyclonal antibody (1:50, 15018-AP, Proteintech, USA) at 4°C. The next day, the slices were rewarmed and incubated with secondary antibody HRP-labelled anti-rabbit (1:500; ab6802; Abcam, USA) antibodies at room temperature for 1h. The sections were then stained using the Diaminobenzaldehyde (DAB) reagent.
The expression of HNRNPR was evaluated using ImageJ software. When the software analysis result produced a negative result, the score was 0; when the result was low positive, positive, or high positive, the scores were 1, 2, and 3, respectively. At least five fields were selected for each slice and evaluated to average the results to arrive at a final score.
Quantitative RT-PCR (qRT-PCR)
The RNA was isolated from cancerous and par-cancerous tissues using the Trizol reagent (Ambion, USA). PrimeScript™ RT Master Mix (TakaRa, Japan) was used for real-time polymerase chain reaction (RT-PCR). SYBR Green Real-time PCR Master Mix (Takara, Japan) was used to detect mRNA expression. The 2−ΔΔCT formula was used to calculate the relative expression multiple of candidate genes in unpaired samples, while HNRNPR/ACTH calculate the expression in paired samples. (HNRNPR: forward primer (5′-3′), GGAGGCAAGAGAAAGGCAGATGG, and reverse primer (5′-3′), GCTGAGCGATGGGTTGGGAAC. ACTB: forward primer (5′-3′), GCACAGAGCCTCGCCTT, and reverse primer (5′-3′), GTTGTCGACGACGAGCG.)
18 F-FDG PET/CT Imaging and Analysis of Related Parameters
PET/CT images of ESCA in patients who fasted for 6 hours, were obtained 50–60 min after 18F-FDG (3.7 MBq/kg) injection, as described in previous experiments [5, 21], using a 64-detector PET/CT scanner (Biograph mCT-S, Siemens, USA). After iteratively reconstructing PET images and computing semi-quantitative parameters, regions of interest (ROIs) were drawn around the tumors. The SUVmax, mean standardised uptake value (SUVmean), total lesion glycolysis (TLG) and metabolic tumour volume (MTV) of each tumor were automatically calculated and recorded.
Gene Set Enrichment Analysis (GSEA) and Bioinformatics Analysis
To determine whether HNRNPR is involved in a diverse and important biological process in ESCA, we used GSEA to evaluate the gene map and associated gene correlation information for HNRNPR samples from the TCGA ESCA data set. The reference gene set was set to c2. cp. v7. 2. symbols. gmt [Curated]. When the FDR (q-value) < 0. 25 and P < 0. 05, the data was considered to be statistically significant.
Enrichment Analysis of HNRNPR Gene Co-Expression Network in ESCA
We investigated the co-expressed genes associated with HNRNPR expression using the stat packet of R software to analyze the TCGA ESCA data sets using Pearson’s correlation coefficient to evaluate for statistical correlation. A volcano and heat map was drawn using the ggplot2 package in R software. Following Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, we performed a visual analysis of the data.
Gene-Gene Interaction and Protein-Protein Interaction(GGI&PPI)
We used the GeneMANIA database (www. genemania. org) to gain a better understanding of R function and genome to generate gene lists and construct interactive networks. When calculating and analyzing the PPI of HNRNPR using the STRING database (www. string-db. org), we selected the minimum required interaction score as the highest confidence level (0. 900).
Correlations between HNRNPR Expression and m6A in ESCA
There were 20 m6A-related genes used to examine the correlation between HNRNPR expression in the TCGA database and the GSE69925 data set, as well as the differential expression between groups with high and low HNRNPR expression [22, 23]. Data were analyzed using the Ggplot2 software package.
Correlations between HNRNPR Expression and Glycolysis in ESCA
Glycolysis-related genes were used to analyze the correlation between HNRNPR expression in the TCGA database and the GSE69925 data set to further investigate the role of HNRNPR in glycolysis [24, 25].
RESULT
Differences in HNRNPR mRNA Expression in Pan-cancer
We compared the transcription levels of HNRNPR in a variety of pan-cancer and normal individuals using the Oncomine online database and TCGA dataset. When we analyzed data from the Oncomine database with a fold change > 1. 5 and a p-value < 0. 001, we discovered that HNRNPR was expressed at a higher level in many cancers than in normal tissue (Fig. 1A).
Furthermore, the difference in HNRNPR expression between cancer and normal tissues was evaluated using the TCGA data set, which revealed that the mRNA level of HNRNPR was increased in a variety of cancers, including ESCA (Fig. 1B).
HNRNPR Expression in ESCA Patients and the Prediction of Poor Patient Survival
The analysis of ESCA patient data from the TCGA and GEO database revealed that HNRNPR expression is significantly higher in cancer tissue than in normal tissue or para-cancerous tissue (Fig. 1C-E).
Then, we used IHC and qPCR to confirm data accuracy and to further investigate the role of HNRNPR in ESCA. The results indicated that, whether at the protein level or mRNA level, HNRNPR expression was significantly higher in para-cancerous tissues, particularly when paired samples of qPCR were analyzed (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.
HNRNPR Co-Expression Networks Indicate the Potential Function of HNRNPR in ESCA
There were 20140 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 of the first 400 co-expressed genes positively correlated with HNRNPR expression were analyzed at p < 0.05 and q < 0.25. 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 an alysis 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.
Correlations 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 data set, these 20 M6-related genes are positively correlated with the expression of HNRNPR (P < 0.05), whereas in the GEO data set, 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.
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 (Vir Like M6A Methyltransferase Associated). Additionally, except for the GSE69925 data set, G6PD, HK1, and SLC2A1 exhibited a negative correlation, whereas the rest showed a positive correlation.
Figure 8B-G shows that HNRNPR is involved in a variety of glycolysis-related pathways, such as influencing FDG absorption.