We downloaded HCC mRNA-seq and clinical information from the TCGA database, including 371 tumor and 50 normal samples. The expression of genes is all logarithmic. At the same time, the expression data of GSE76427 and the corresponding clinical data from the GEO database (https://www.ncbi.nlm.nih.gov/geo) served as validation data set. GSE76427 contains 115 tumor tissues and 52 adjacent tissues. The non-zero survival time was retained for clinical analysis in the TCGA and GSE76427 data sets.

Clinical evaluation of RAC1 integrates the expression data of RAC1 and the survival data of the patients to evaluate the prognostic value. We applied t test to infer the expression changes of RAC1 in various stages of HCC. Then, we analyzed the correlation between RAC1 and the overall survival (OS) through the Kaplan-Meier (K-M) curve. The log-rank test was used to compare the prognostic differences between the groups. The diagnostic capability of RAC1 is measured by the area under the time-dependent receiver operating characteristics (tROC) curve (AUC).

The Cox regression model establishes the relationship between the survival rate of patients and several variables and evaluates the effects of these factors on survival time. When p < 0.05, the impact on survival time is statistically significant, identifying the expression profiles of 338 immune-related genes. The survival analysis aims to study the relationship between the variable X and the survival function (cumulative survival rate) S (t, X). X = (X_1,…, X_m) is a vector, and S (t, X) is affected by many factors. The traditional method considers the regression equation, that is, the influence of variables X_i on S (t, X). However, the data in the survival analysis contains censored data, and the time t usually does not meet the requirements of normal distribution and homogeneity of variance. These reasons make it difficult to study the above-mentioned relationship with the general regression method. Therefore, we use the Cox regression model as a specific tool (Goerdten et al., 2020). The Cox regression is implemented through the “survival” package in R.

Support vector machine (SVM) is a class of generalized linear classifiers that classify data binary in a supervised learning manner. Its decision boundary is the maximum margin hyperplane that is solved for learning samples (Huang et al., 2018). CIBERSORT performed deconvolution based on ν’s support vector machine (ν-SVR) method, a support vector machine (SVM) optimization method for binary classification problems. Based on deconvolution and linear support vector regression principles, CIBERSORT evaluates 22 immune cell molecular subtypes and uses v-SVR to estimate immune scores based on the expression matrix of IRGs. Therefore, the SVM algorithm is used in this paper through CIBERSORT.

CIBERSORT is a deconvolution tool that links immune cell infiltration to gene expression. It includes the LM22 gene signature file, which contains 547 genes used to identify 22 immune cell subtypes. We uploaded the mRNA expression data to the CIBERSORT portal (https://CIBERSORT.stanford.edu/). The gene expression data can then be standardized by running 1,000 permutations through the default feature matrix. Finally, we applied WilCoxon test to determine the differentially infiltrated immune cells in tumors. The deconvolution algorithm CIBERSORT (Newman et al., 2015) is used to determine the immune-related features from 365 labeled immune genes, and to quantify the relative score of each immune cell type. The method is used to evaluate the relative proportions of 22 tumor-infiltrating immune cell profiles based on expression files, including B cells, T cells, natural killer cells, macrophages, dendritic cells, and bone marrow subpopulations. And we used gene expression data to estimate the abundance of member cell types in the mixed cell population. Monte Carlo sampling is used to obtain the p-value for each deconvolution sample through CIBERSORT algorithm. The number of permutations is set to 1,000, and p < 0.05 is considered significant. The immune cell matrix of each sample was obtained in the transcriptome data through the CIBERSORT algorithm.

The R platform’s “limma” software package is used to find genes that are differentially expressed in the high expression community vs. the low expression progenitor of RAC1 (p 0.05 and (log 2 FC) > 1). Then we built a co-expression network of differentially expressed genes using the R package “WGCNA.” WCGNA clusters IRGs with similar expression patterns to construct a scale-free gene co-expression network and analyzes the correlation between modules and specific phenotypes. The important IRG modules are identified according to the module-IRG, module-module correlation. WGCNA is a robust algorithm highlighted by the modular clustering of genes and the association analysis between the modules and clinical traits.

First, we used a hierarchical cluster analysis on the expression profile. S = [ S X Y ] = [ | c o r ( X , Y ) | ]

To create the adjacency matrix A, calculate the correlation index an XY of any gene pair using the square of the correlation coefficient S XY. A = [ a X Y ] = [ p o w e r ( S X Y , β ) ] = [ | S X Y | β ]

Then, a topological overlap matrix (TOM) is constructed, T O M = [ ω X Y ] = [ l X Y + a X Y min { k X , k Y } + 1 − a X Y ]

The R software package “survival v3.1-8” was used to perform univariate Cox regression on the immune module genes in the TCGA cohort. p < 0.05 is used to screen for genes significantly related to the overall survival rate of patients. We used the R software package “glmnet v3.0-2” to perform LASSO analysis to optimize the results further. LASSO appeared as a constraint on the objective function. min ω ∑ j = 1 m ( y j − ∑ i = 1 n x j i ω i ) 2 , s . t . ∑ i = 1 n | ω i | ≤ λ

The risk score (RS) referred to R S = ∑ i = 1 n β i × E x p i where n denoted the prognostic signature genes. The LASSO coefficient of a gene is represented by i. The expression of a gene is represented by Exp _i .

LASSO COX performs collinearity processing on the filtered genes. Thus, the algorithm has advantages in processing high-dimensional data. The screening performed 10-fold cross-validation, and the survival status and time were used as dependent variables to analyze the influence of multiple IRGs on the dependent variables. Finally, the LASSO coefficient is determined when the lambda is the minimum, forming the IRG signatures.

The “R v3.6.1″ framework is used for all statistical analysis R packages. Statistical significance was described as a p-value of less than 0.05.

The TCGA cohort collected the pathological characteristics of 365 patients to evaluate the clinical manifestations. We studied the gene expression profile and clinical information of 115 patients with a non-zero survival time of GSE76427. RAC1 expression was substantially higher in tumor tissues than in normal ones (p = 1.9e-12; Figure 1A). Kaplan-Meier results showed that the increase in RAC1 expression was closely linked to the poor patient survival rate (p < 0.001, Figure 1B). In addition, RAC1 has excellent capabilities in HCC (1-, 3-, 5-year AUC = 0.690, 0.635, 0.661; Figure 1C). The GSE76427 cohort also confirmed that RAC1 has a good prognostic ability (Figure 1). Interestingly, RAC1 expression was significantly correlated with the grade, T, and stage (Figures 1D–F).

In the TCGA cohort, T (HR = 1.0960, 95% CI = 0.8442–2.2774, p < 0.001), Stage (HR = 1.0045, 95% CI = 0.8703\−2.3584, p < 0.001) and RAC1 (HR = 1.7852, 95% CI = 1.3775–2.3136, p < 0.001) was significantly correlated with the overall survival of the patient. Multivariate Cox regression found RAC1 (HR = 1.590, 95% CI = 1.209–2.09, p < 0.001) was an independent progostic marker of HCC (Figure 2A). In the GSE76427 cohort, we found that RAC1 (HR = 3.1108, 95% CI = 1.1863–8.1572, p = 0.0210) significantly affected the survival of patients. In addition, multivariate Cox regression found that Stage III (HR = 3.32, 95% CI = 1.083–10.21, p = 0.0358) and RAC1 (HR = 3.590, 95% CI = 1.347–9.54, p = 0.0106) are risk factors for HCC (Figure 2B).

The CIBERSORT deconvolution algorithm obtained the infiltration level of 22 immune cells. Grouping 365 tumor samples with high and low expression of RAC1, we found differences in the infiltration of nine immune cells in the TCGA cohort (Figure 3A). In the GSE76427 cohort, we found seven differences in immune cell fluctuations (Figure 3B). Macrophages M0, M1, and M2 were identified by Wayne analysis (Figure 3C). We speculate that RAC1 mediates macrophages to play a critical role in the development of HCC. The expression of immune checkpoint inhibitors (PD-1 and PD-L1) was significantly correlated with macrophage M1 (Figure 3D).

We obtained 1,417 differentially expressed genes (DEGs), of which 80 were upregulated and 1,337 were downregulated (Figure 4A). WGCNA analysis was performed based on the mRNA expression data of 1,417 DEGs. The soft threshold is selected as 5 to conform to the scale-free network rule through hierarchical clustering (Figure 4B). A total of four gene modules were identified using the clustering criteria of minModuleSize = 30 and mergeCutH8 = 0.25 (Figure 4C). TOM is used to describe the pair-wise relationship between genes. In the topological overlap matrix (TOM), each row and each column corresponds to a gene. Light colors indicate low topological overlap, and gradually darker ones indicate higher topological overlap. The dark squares along the diagonal correspond to the modules. The gene tree diagram and module allocation are shown on the left and top. The four gene modules are gray, turquoise, yellow, blue, and brown modules. The correlation analysis between the module and the clinical phenotype showed that the turquoise module is related to PD-1 (cor = 0.52, p = 6e-24), PD-L1 (cor = 0.3, p = 2e-08), and macrophages M0 (cor = 0.19, p = 6e-04) is significantly correlated (Figure 4D). Therefore, the turquoise module (n = 838) was identified as immune-related genes (IRGs) for subsequent analysis.

We obtained 388 IRGs that significantly correlated with the survival through univariate Cox analysis. Subsequently, LASSO Cox analysis was used to optimize the results further. We finally obtained 11 genes (CCDC136, CTSE, CXCL5, CXCL6, FAM188B, GALNT6, HOMER3, KCNF1, KCNQ3, NRCAM, PTP4A3), which are significantly related to the survival rate of patients (Figures 5A–C). Figure 5 shows the identification of immune-related genes. LASSO COX determines the number of factors by introducing shrinkage penalties and limiting the coefficients. With the continuous selection and simulation of the number of features, the best model and the simplest model are finally obtained. Subsequently, we further calculate the final score of the model. The final risk score is obtained by multiplying the expression of each gene with its corresponding coefficient and adding them together. And, with the median risk score of all patients as the critical value, the patients were divided into high-risk and low-risk two groups.

The risk scoring formula is as follows: RS  = CCDC 136 *0 . 0228 + CTSE*0 . 0176 + CXCL 5 *0 . 0244 + CXCL 6 *0 . 006 + FAM 188 B*0 . 0016 + GALNT 6 *0 . 0326 + HOMER 3 *0 . 001 + KCNF 1 *0 . 0056 + KCNQ 3 *0 . 0916 + NRCAM*0 . 0298 + PTP 4 A 3 *0 . 0247

Patients were divided into low- and high-risk groups (152 vs. 153) according to the risk score in the TCGA cohort. Compared with high risk, low risk can significantly improve the prognosis of patients (p = 0.00012). The 1-, 3-, and 5-year AUC (0.706, 0.638, and 0.642) confirmed the signature’s good prognostic performance (Figure 6A). In addition, the K-M analysis (p = 0.029) and ROC analysis also confirmed the above results in the validation cohort. The 1-, 3-, and 5-year survival rates of the risk prognostic feature prediction validation cohort were 0.700, 0.701, and 0.750, respectively (Figure 6B).

We perform multivariate Cox analysis on the immune and risk characteristics obtained by the algorithm above. M0 macrophage (HR = 6.050, 95% CI = 1.3113–27.91, p = 0.021) and risk score (HR = 3.591, 95% CI = 2.038–6.33, p < 0.001) were significantly correlated with the survival (Figure 7A). In addition, M0 macrophage (HR = 0.19, 95% CI = 0.061–0.61, p = 0.005), macrophage M2 (HR = 0.11, 95% CI = 0.040–0.31, p < 0.001) and risk score (HR = 0.32, 95% CI = 0.229–0.45, p < 0.001) was significantly correlated with RAC1 (Figure 7B). Therefore, the M0 macrophage and risk score were determined as prognostic factors.