Gene expression and clinical data were extracted from multiple databases ( Supplementary Table 1 ) and RAB22A expression levels from RNA-seq data (TPM) of patients with HCC were analyzed. The Supplementary Materials and Methods (17–19) presents detailed information on the included data.

The First Affiliated Hospital of Harbin Medical University provided 30 matched sets of HCC and nearby non-tumor liver samples from patients undergoing hepatectomy between February 2020 and June 2022. This project was approved by the First Affiliated Hospital of Harbin Medical University’s Ethics Committee.

Total proteins and total RNA were extracted from HCC samples. Details of the experimental procedures are provided in the Supplementary Materials and Methods .

We used the single sample gene set enrichment analysis (ssGSEA) method (20, 21) and TIMER database (22, 23) to investigate the relationships between RAB22A expression and immune cell infiltration, as detailed in the Supplementary Materials and Methods .

Enrichment analyses of relevant functional pathways were performed using the GO and KEGG databases ( Supplementary Tables 3 , 4 ) and GSEA ( Supplementary Tables 5 , 6 ), as detailed in the Supplementary Materials and Methods (24, 25).

Multiple databases were used to predict and screen the lncRNA-miRNA-mRNA (RAB22A) ceRNA network online. Details are provided in the Supplementary Materials and Methods .

We created the protein–protein interaction (PPI) network using the Search Tool for the Retrieval of Interacting Genes (STRING) database ( Supplementary Table 6 ) (26, 27), as detailed in the Supplementary Materials and Methods .

The R package (version 3.6.3) was used for statistical analyses and plotting. RAB22A expression in unpaired and paired samples was analyzed using the Wilcoxon rank sum test, and Wilcoxon signed rank test, respectively, with the pROC (1.17.0.1) package for ROC analysis. The RAB22A expression level was analyzed by querying the GEO, TIMER, and UALCAN databases (18). Using the KM method and log-rank test, we compared the differences in 10-year OS, DSS, and PFI between patients with high RAB22A expression and those with low RAB22A expression in TCGA. Cox analysis was used to determine the correlation between RAB22A expression and clinical features. p < 0.05 was considered to indicate significance.

First, we examined the RAB22A expression levels in different malignancies by assessing TCGA databases. RAB22A was highly expressed in 33 malignant tumors, including HCC ( Figure 1A ). In addition, RAB22A was highly expressed in the GEO datasets GSE121248, GSE87630, GSE76427, GSE84005, GSE57957, and GSE39791 HCC samples (p < 0.001) ( Figures 1B–G ). Western blot analysis of human normal liver cells (L02) and HCC cells (Hep G2, SK-Hep1, Huh7, HCCLM3, and MHCC97-H) validated the high expression of RAB22A in HCC cell lines ( Figure 1H ). The same results were obtained through qRT-PCR (p < 0.001) ( Figure 1K ). Next, we extracted 30 pairs of proteins from HCC and adjacent tissues and analyzed them using western blotting, which revealed that RAB22A was highly expressed in the former ( Figure 1I ). Results of western blot analysis of the liver and adjacent tissues are shown in Supplementary Figure 1 . The high RAB22A mRNA expression levels in HCC tissues were further substantiated using qRT-PCR (p < 0.001) ( Figure 1L ). Immunohistochemistry (IHC) results also verified that RAB22A was upregulated in HCC tissues ( Figure 1J ). Finally, a receiver operating characteristic (ROC) curve was created. The ROC curve enclosed by the axes is the area below the curve (AUC). The AUC for RAB22A was 0.891, suggesting its remarkable diagnostic value for HCC ( Figure 1M ).

Using the UALCAN database to perform subgroup analysis of numerous pathological characteristics, we found that RAB22A transcript levels were elevated in patients with HCC. ( Figure 2A ). The sub-group analysis of cancer stage, ethnicity, sex, age, weight, tumor grade, and TP53 mutation showed that the expression of RAB22A in HCC patients was significantly higher than that in the normal group ( Figures 2B–H ).

Logistic regression analysis showed that the increased expression of RAB22A in HCC was significantly correlated with sex (OR = 0.627 for male vs. female, p = 0.036), weight (OR = 0.567 for weight > 70kg vs. ≤ 70kg, p = 0.009), histological grades (OR=1.611 for G3 and G4 vs G1 and G2, p = 0.028), and tumor status (OR = 1.619 for with tumors vs. tumor free, p = 0.026). Conversely, RAB22A expression was not associated with age, M stage, T stage, N stage, height, BMI, AFP, or vascular invasion ( Table 1 ).

Next, we collected data from TCGA database to determine the clinicopathological parameters of RAB22A in different patients with HCC. Detailed information on the clinical data is provided in Table 2 . After excluding cases without the necessary clinical data, 374 cases with a median age of 61.5 (range: 49.25−70.00) years and male preponderance of 67% were included. High expression of RAB22A in HCC was positively associated with tumor status (tumor-free vs. with tumor, p = 0.033), sex (female vs. male, p = 0.047), weight (≤ 70 vs. > 70, p = 0.012), and histological grade (grades 3 and 4 vs. grades 1 and 2, p = 0.031). These results indicate that the overexpression of RAB22A in HCC is closely related to the clinicopathological characteristics.

Kaplan–Meier survival curves were analyzed to determine the connection between RAB22A expression and overall survival (OS), disease-free survival (DSS), and progression-free interval (PFI) in the prognosis of patients with HCC. Increased levels of RAB22A expression were inversely related to prognosis ( Figures 3A–C ). Additionally, subgroup analysis was performed on patients with low RAB22A expression and AFP < 400, and these patients had better OS, DSS, and PFI prognosis ( Figures 3D–F ). However, the groups with AFP (ng/mL) > 400 showed no significant differences ( Supplementary Figures 2A–C ). The high expression of RAB22A in stage M0 liver cancer was associated with poor OS, DSS, and PFI in a subgroup of patients ( Figures 3G–I ). The subgroups of T3 versus T4, stages III vs. IV, and tumor versus tumor-free status had significantly worse OS ( Supplementary Figures 1D–F ). Finally, we compared predictive variables in patients with HCC obtained by univariate regression analysis to those obtained via multivariate survival analysis (OS) ( Supplementary Table 7 ). Pathologic stage (stages I and II compared with stages III and IV; p <0.001), tumor size (T stages 1 and 2 versus T stages 3 and 4; p < 0.001), metastatic spread (M stages 0 and 1; p = 0.017), and tumor status (without or with tumor; p < 0.001) were highly significant in the univariate analysis. The multivariate analysis showed that with tumor (p = 0.014) was significant, suggesting that it is an independent risk factor.

GO and KEGG pathway co-expression analyses of RAB22A-related genes in liver cancer mRNA sequencing data with 371 patients from the TCGA were performed using the functional module of Linkedomics. The top 50 marker genes and their connections with RAB22A expression are displayed on the heat map ( Figures 4A, B ; Supplementary Table 8 ). These findings revealed a widespread effect of RAB22A on the transcriptome.

Next, we conducted an enrichment analysis using the GO and KEGG databases to support the concept that RAB22A-related DEGs play a biological role in HCC ( Figures 4C, D ). The results of GO analysis showed that these DEGs were related to biological processes (BP), cellular components (CC), and molecular functions (MF). In the GO analysis, DEGs were enriched in diverse biological pathways, including proteasomal protein catabolic process, ncRNA processing, ribosome ribosomal subunit, protein serine/threonine kinase activity, and protein serine kinase activity. In the KEGG analysis, DEGs were highly concentrated in endocytosis and non-alcoholic fatty liver disease. GSEA was used to analyze the biological functions related to RAB22A expression.

The later criteria were enrichment score | NSE | > 1 (p < 0.05), according to which the five most positively relevant signal pathways were selected. GO analysis revealed that RAB22A expression was strongly positively correlated with the processes of homophilic cell adhesion via plasma membrane adhesion, immunoglobulin, T cell receptor, and plasma membrane signaling receptor complex ( Figure 4E ). The expression of RAB22A was inversely linked to that of co-translational proteins that bind to the membrane, the cytosolic ribosome, the structural components of the ribosome, the ribosomal subunit, and the nonsense-mediated decay of nuclear-transcribed mRNA catabolic processes ( Figure 4H ). KEGG analysis revealed that RAB22A expression was most strongly negatively connected with axon guidance, extracellular matrix receptor interaction, focal adhesion, FCγR-mediated phagocytosis, and the interaction with neuroactive ligand receptors ( Figure 4F ). The ribosome, Parkinson’s disease, retinol metabolism, oxidative phosphorylation, and complement and coagulation cascades were the top five most negatively correlated pathways ( Figure 4I ). REACTOME pathway analysis determined that phospholipids play a role in phagocytosis, Fc gamma receptor FCGR-dependent phagocytosis, CD22-mediated B cell receptor (BCR) regulation, and FCGR activation, and that second messenger is activated by BCR antigen that all positively correlated with RAB22A expression ( Figure 4G ). Eukaryotic translation initiation, eukaryotic translation elongation, translocation, the response of eukaryotic initiation factor 2 alpha subunit kappa B cyclin N2 to amino acid deprivation, co-translational protein of SRP-dependent targeting to the membrane, and nonsense-mediated decay were all negatively correlated with RAB22A expression ( Figure 4J ).

The PPI network of co-expressed genes conforming to the STRING conditions was assembled and visualized using Cytoscape, and analysis of the interactions among 108 DEGs in the HCC group was conducted. A total of 51 proteins and 534 edges were screened ( Figure 5A ; Supplementary Table 9 ).

After screening 12 nodes and 212 edges, a primary gene cluster with a total score ≥ 14,000 was discovered ( Figure 5B ). Finally, we screened the top 10 central genes, namely RAB22A, RABGEF1, VPS45, VPS18, VPS11, MON1A, VPS39, VPS16, ZFYV20, and VPS8 ( Figure 5C ).

M6A methylation affects the development of HCC (24–27). The expression of RAB22A was compared with that of the 23 M6A methylation genes reported in the literature to verify this conclusion ( Figure 6A ). RAB22A expression was closely connected with that of the 23 m6A-related genes in HCC ( Figures 6B–X ). Moreover, groups were formed according to RAB22A median expression. By analyzing the differences in the 23 m6A methylation genes in RAB22A between the high- and low-expression groups of patients with HCC, we observed that the expression levels of all genes in the RAB22A high-expression group were upregulated ( Figure 6Y ). Overall, we observed an obvious relationship between m6A methylation and RAB22A expression levels in HCC.

Increasing evidence has demonstrated the regulatory effect on the lncRNA-miRNA-mRNA ceRNA network in HCC. The Venn diagram showed 41 overlapping miRNAs in the Targerscan, starBase, and MiRDB databases ( Figure 7A ). Five human-derived miRNAs (miR-328-3p, miR-3163, miR-2114-5p, miR-664b-3p, and miR-204-5p) were verified to negatively correlate with RAB22A expression ( Figure 7B ). The expression of RAB22A and target microRNAs is displayed as a scatter plot ( Figures 7C–G ). We consulted the Rnalnter and starBase databases to predict lncRNAs that can have a mutual effect on target miRNAs (miR-204-5p and miR-328-3p) ( Figures 7H–I ). The expression levels of lncRNAs and miRNAs were inversely correlated, which accounted for the mutual influence between the two. We used the starBase database to filter and identify lncRNAs that were adversely associated with the two target miRNAs in HCC. Nine HCC-related ceRNA regulatory networks were constructed ( Figure 7J ).

Using the ssGSEA method, we verified the strong connection between RAB22A and immune cells ( Figure 8A ). The expression of RAB22A was positively connected with T helper cells, Tcm cells, and Th2 cells (p < 0.001) but negatively with cytotoxic cells, DCs, and pDCs (p < 0.001) ( Figures 8B–G ). RAB22A may be heavily involved in the T-cell immune response to HCC. Moreover, RAB22A expression in HCC correlated with various immune cell markers ( Table 3 ). In the M2 macrophages in HCC, we found that RAB22A expression was substantially relevant to the expression of the immunological markers CD163, VSIG4, and MS4A4A. These results indicate that RAB22A caused the macrophages in HCC to adopt an M2 phenotype. The expression of RAB22A was substantially linked to 66 immunological markers, including CD8A, CD3D, and T-bet, in an analysis of functional T-cell immunity indicators. Furthermore, RAB22A expression was linked to immunological markers for B cells, T cells, TAMs, and neutrophils ( Table 3 ). The TIMER database was utilized to determine whether RAB22A expression in HCC was connected with immune cell invasion levels. The results indicated that the CNV of RAB22A was related to the level of neutrophil infiltration. ( Figure 8H ). Subsequently, the infiltration of macrophages, T helper cells, Tcm, and Th2 cells increased (p < 0.001) in the RAB22A high-expression group; however, cytotoxic cells, DCs, and pDCs decreased (p < 0.001) ( Figure 8I ). These results verified that the increased expression of RAB22A in HCC is tightly linked with the infiltration of immune cells.