For mRNA expression analysis, RNA seq data of TCGA-LIHC dataset, which was originally sourced from Broad GDAC Firehose (http://gdac.broadinstitute.org/) (Wheeler and Roberts, 2017) was extracted using UCSC XENA webserver (Goldman et al., 2020). Clinicopathological and molecular characterstics of the TCGA-LIHC dataset has been given in Table 1. Microarray gene expression data from multiple studies was accessed through the TNMplot webserver (Bartha and Győrffy, 2020). This web server hosts data from multiple HCC studies, where gene expression data has been normalized for all available studies and can be used for comparison between the collective groups of all normal samples with tumor samples. Multi-Omics dataset of hepatocellular carcinoma released by Clinical Proteomic Tumor Analysis Consortium (CPTAC) (https://cptac-data-portal.georgetown.edu/cptacPublic/) was utilized to analyze both mRNA and protein levels of RAP genes.

DNA methylation of RAP genes in TCGA cancer dataset was estimated and visualized using MEXPRESS web server (https://mexpress.be) (Koch et al., 2015; Koch et al., 2019) and TCGA Wanderer (Díez-Villanueva et al., 2015). The MEXPRESS web server uses DNA methylation data of cancer and normal tissues from TCGA datasets, which were originally developed on the Illumina Human Methylation 450 BeadChip platform. The predesignated methylation probes for each gene were taken into consideration.

Kaplan-Meier survival analysis was performed using the tool available in the KM-plotter (Nagy et al., 2018). For Kaplan-Meier analysis, patients were distributed in high and low expression groups based on median expression value as a cut-off point for each gene. For survival analysis using univariate and multivariate Cox proportionate hazard model, RAP2A gene expression was taken as a continuous variable with multiple survival parameters for the TCGA-LIHC dataset, as recommended (Liu et al., 2018).

Similarly, whole transcriptome correlations of RAP2A in the TCGA-LIHC study were downloaded from the cBioPortal website (https://www.cbioportal.org/) (Cerami et al., 2012; Gao et al., 2013). After applying a filter for a cutoff of FDR corrected p-value of 0.05 for Spearman’s r-value, 10,980 genes with Spearman’s correlation q value <0.05 were filtered and used for gene set enrichment analysis in GSEA software (Broad Institute, http://www.broad.mit.edu/gsea/). Hallmark gene set (version 7.1) (Subramanian et al., 2005) from predefined molecular signature database was used as a reference gene set for pathway enrichment (Liberzon et al., 2015).

Tumor immune estimation score (TIMER) webserver (https://cistrome.shinyapps.io/timer/), which utilizes the RNA sequencing data from TCGA for estimation of correlation between gene expression and level of immune cells, present in the tumor samples (Li et al., 2017). We utilized TIMER to calculate the association between RAP2A gene expression with infiltration of six immune cells including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells in TCGA-LIHC datasets. Default parameters were used in the TIMER database for the gene-specific analysis module.

Further, CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts) analysis data of TCGA-LIHC were extracted from a previously published study (Chen et al., 2018; Thorsson et al., 2018). This provided relative fractions of 22 different immune cells from a mixture of gene expression profiles (TCGA-LIHC study) and was used to correlate with RAP2A expression using Spearman’s correlation test. A total of 360 HCC samples were available with both gene expression data and CIBERSORT analysis estimated fractions of immune cells. Heatmap of the immune cell profiling data was generated along with hierarchical clustering using HemI (Deng et al., 2014). The default parameters of hierarchical clustering using the average linkage method and Pearson distance were used.

Data analysis was performed using Graphpad (version 6) and Stata software (version 11). Mann-Whitney U-test was used for comparison among histological subtypes, molecular subtype and grades (***p < 0.001; **p < 0.01; *p < 0.05; ns, p > 0.05). Pearson correlation analysis was used to determine the association of DNA methylation of RAP2A to its expression in the TCGA-LIHC dataset. Kaplan-Meier survival analysis was performed using the log-rank test. A p-value < 0.05 was considered statistically significant.

Further, RNA sequencing data from TCGA-LIHC study was utilized to compare RAP gene expression in tumor tissues with both tumor adjacent normal tissues from the same dataset and with non-tumor associated normal hepatic tissue. RAP1A, RAP1B, RAP2A and RAP2B exhibited significant higher expression in tumors compared to other two groups (Figures 1A,C,E,G). Although RAP2C expression was higher in tumors compared to adjacent normal tissues, but both these groups exhibited lower expression of RAP2C compared to normal tissues from GTEx (Figure 1I). Furthermore, tumor adjacent normal tissues also exhibited higher expression of RAP1A and RAP2B, while no difference was observed for RAP1B and RAP2A. Comparison between 50 paired normal and tumor tissues from TCGA-LIHC also revealed that all RAP genes exhibit higher expression in tumor tissues compared to tumor adjacent normal tissues (Figures 1B,D,F,H,J). Among all RAPs, RAP2A displayed most robust upregulation in tumor tissues in TCGA dataset (Figure 1F). Further, we utilized multiomics data of hepatocellular carcinoma developed by CPTAC study, where both mRNA and proteomic data was available. The expression analysis in CPTAC data also suggested that expression of RAP genes differ between normal and tumor tissues both at the mRNA and protein level. In CPTAC data also, RAP2A exhibited most robust upregulation of mRNA and protein levels in tumor tissues compared to normal tissues, while expression of RAP2C was found to be reduced in tumor tissues compared to normal tissues (Supplementary Figure S1).

We further performed receiver operating characteristic (ROC) curve analysis to determine potential of RAP gene expression in differentiating tumor tissues from normal liver tissues. Interestingly, among five RAP genes, RAP2A exhibited highest area under curve (AUC) of 0.8676 in TCGA-LIHC mRNA data (Figure 2A). Similarly, analysis of CPTAC mRNA data also suggested highest AUC of RAP2A (AUC: 0.9173, Figure 2B) compared to other RAP genes. Interestingly, analysis of AUC in CPTAC protein expression data revealed that RAP2C exhibit highest AUC of 0.8445 followed by RAP2A with AUC of 0.8172 (Figure 2C).

Furthermore, expression data of RAP gene family in hepatocellular carcinoma tissues and normal tissues from multiple other datasets was assessed through TNMplot web server. The analysis revealed that RAP1A, RAP1B, RAP2A, and RAP2B genes exhibit significantly higher expression in HCC tissues compared to normal tissues in comparison of both available paired (Supplementary Figure S2, left panel) and unpaired tissues (Supplementary Figure S2, right panel). However, RAP2C did not exhibit significant difference in expression in paired tissue analysis (Supplementary Figures S2J). Considering robust upregulation of RAP2A in tumors, and its established involvement in cell migration, we compared expression of RAP2A in metastatic tissues with both normal and primary tissues, which revealed highest expression of RAP2A in metastatic tissues compared to other two groups (Figure 2D).

We further assessed the association of RAP genes with clinicopathological including pathological age, gender, stage, tumor grade, blood AFP levels. Among all RAP genes, higher expression of RAP1B was associated with advanced-stage (Figure 3A). Higher expression of RAP2A and RAP2B, and low expression of RAP1A was associated with advanced grade (Figure 3B). High RAP2A expression was associated with younger age (<50 years, Supplementary Figure S3A) and female gender (Supplementary Figure S3B). Higher expression of RAP2A was also associated with increased AFP levels (Figure 4A). A history of alcohol consumption was associated with lower levels of RAP2A and RAP2C expression (Figure 4B).

To further determine whether the expression of RAP genes is associated with genetic alterations in HCC, we compared their expression in tumors with mutated or wild type TP53, CTNNB1, ALB, PCLO, and LRP1B. TP53 mutation was observed to be associated with higher expression of RAP1A, RAP1B, RAP2A, and RAP2B (Figure 5A). Further, CTNNB1 mutation was significantly associated with reduced levels of RAP1B, RAP2A, and RAP2B expression (Figure 5B). No RAP gene exhibited association with PCLO and ALB mutation status (Supplementary Figures S4A,B respectively), while higher expression of RAP1A was associated with LRPB1 mutant tumors (Supplementary Figure S5).

To determine the potential role of copy number alterations and DNA methylation in the regulation of RAP2A expression in HCC, we utilized a TCGA-LIHC study where copy number variation, DNA methylation, and gene expression data were available. The DNA methylation data in TCGA study was developed on “Illumina HumanMethylation450 Beadchip” platform, where representative CpG sites from different regions of most genes are captured. Interestingly, RAP2A gene expression was reduced in tumor tissues and exhibited a negative correlation with DNA methylation at several sites within the RAP2A promoter regions and gene body (Figures 6A,B). A similar association was also observed for normal liver tissues (Supplementary Figure S6). We observed that in both normal and tumor tissues, DNA methylation at an intragenic region represented by cg03608515 was most negatively correlated with gene expression, suggesting this region, but not promoter region is the major regulatory site for the expression (Figure 6C). Furthermore, a comparison of 47 paired normal and tumor tissue also revealed significantly reduced methylation levels of cg03608515 in tumor tissues., these results strongly suggest the role of DNA methylation in aberrant expression of RAP2A in HCC. Additionally, the expression of RAP2A was also positively correlated with its copy number (r = 0.450, p < 0.001). Further, analysis of CNV data revealed frequent alterations in RAP2A copy number in HCC tissues was associated with its higher expression with copy number gain (Figure 6D, Kruskal-Wallis test, p < 0.0001).

To determine the association of RAP gene family expression with patient prognosis, we utilized the TCGA-LIHC dataset where information for overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) was available. We performed survival analysis by constructing a Kaplan-Meier plot for all RAP genes using median expression levels for allotting patients into high and low groups. We observed that higher expression of RAP2A was significantly associated with poor OS (HR = 1.72, CI = 1.21–2.45, p = 0.0023, Figure 7A) and DSS (HR = 1.9, CI = 1.2–2.99, p = 0.005, Figure 7B), while no significant association was observed with DFI and PFI (Figures 7C,D, respectively). In light of the high positive correlation of RAP2A with other RAP genes, we also assessed their association with patient survival (Supplementary Figure S7). Among other RAPs, higher expression of RAP1A and RAP1B was also associated with poor overall survival (Supplementary Figures S7A,B). We further performed univariate and multivariate survival analysis for RAP2A expression and other clinicopathological features, such as age, gender, stage, grade, alcohol intake history, radiotherapy status, and embolization status using Cox proportionate hazard model. Interestingly, higher RAP2A expression was also associated with poor OS, DSS, and PFI in both univariate and multivariate survival analysis (Tables 2, 3, respectively). This suggested that RAP2A expression is independently associated with poor outcome in HCC patients.

To determine RAP2A associated cancer-related pathways, gene expression data of the TCGA-LIHC study was used. GSEA analysis revealed that RAP2A expression is positively correlated with cell cycle associated pathways such as mitotic spindle (Figure 8A), G2M checkpoint (Figure 8B), and E2F targets (Figure 8C) besides protein secretion (Figure 8D). Further, negatively correlated genes were enriched in metabolism associated pathways, such as oxidative phosphorylation (Figure 8E), xenobiotic metabolism (Figure 8F), fatty acid metabolism (Figure 8G), bile acid metabolism (Figure 8H), adipogenesis (Figure 8I), reactive oxygen species (Figure 8J) and others such as coagulation (Figure 8K), peroxisome (Figure 8L), interferon-alpha response (Figure 8M), DNA repair (Figure 8N) and Myc target genes (Figure 8O).

Considering the previously described role of RAP genes in immune cell functions (Carvalho et al., 2019a), we analyzed the association of RAP2A expression with the level of immune cell infiltration. Using the TIMER tool, we determine tumor purity normalized spearman correlation of RAP2A expression with infiltration level of six different immune cells. This analysis revealed a positive correlation between RAP2A expression with B cells (r = 0.3, p = 1.37e-08), CD8+ T cells (r = 0.237, p = 9.06e-06), CD4+ T cells (r = 0.474, p = 1.16e-20), macrophages (r = 0.469, p = 4.56e-20), neutrophils (r = 0.374, p = 7.19e-13), and dendritic cells (r = 0.401, p = 1.36e-14) in HCC (Figure 9A). Furthermore, we utilized CIBERSORT analysis to determine the association of RAP2A gene expression with the relative abundance of 22 different types of immune cells in the TCGA-LIHC dataset (Figure 9B, Supplementary Table S2). Among immune cells, RAP2A expression was positively correlated to CD4 Memory Resting T cells, resting dendritic cells, neutrophils, M0 type macrophages, and naïve B cells, while it exhibited negative correlations to monocytes, activated NK cells, CD4 naïve T cells, CD8 T cells.