GDI2 is a novel diagnostic and prognostic biomarker in hepatocellular carcinoma

Clinical characteristics

The clinical characteristics of 371 hepatocellular carcinoma (HCC) patients from TCGA were collected as shown in Table 1, including gender, race, TNM stage, pathologic stage, vascular invasion and tumor status, as well as HCC-specific index as adjacent hepatic tissue inflammation, fibrosis ishak score and Child-Pugh grade. In present study, a total of 121 female patients (54.2%) and 250 male (54.2%) patients were analyzed, including 184 white patients (49.6%) and 175 non-white patients (47.2%; 158 Asian patients and 17 Black or African American patients). The tumor status involved 201 tumor free (54.2%) cases and 151 cases with tumor (40.7%), while 109 with (19.4%) and 206 without (55.6%) vascular invasion. Stage I disease was found in 171 patients (46.1%), and Stage II, III, IV in 85 (22.9%), 85 (22.9%), and 5 (1.3%) patients, respectively. Most tumors were distributed among 48.8% T1 Stage (n=181), 25.3% T2 (n=94), 21.6% T3 (n=80), and 3.5% T4 (n=13). As not all the HCC patients could offer intact clinical data, the censoring values made the total number fail to sum up to 371 cases for each clinical count.

Identification of differentially expressed genes (DEGs) in HCC

To elucidate whether GDI2 was positively correlated with HCC occurrence, the GDI2 expressions were compared between 371 HCC cases and 50 normal tissue cases via RNAseq TMP data from TCGA combined with GTEx, which showed much higher expression of GDI2 in tumor cases than in normal cases (P < 0.001, Figure 1A). Among them, 50 tumor-and-adjacent paired samples were also showed high expression of GDI2 in tumor compared with paired normal tissues (P < 0.001, Figure 1B). To further define the oncogenicity of GDI2, comparison of GDI2 expression between tumor and normal specimens was also made in pan-cancers from TCGA-GTEx database. It was shown that GDI2 significantly expressed in most of the 33 kinds of solid tumors, such as hepatocellular carcinoma (LIHC, HCC), lung adenocarcinoma (LUAD), cholangiocarcinoma (CHOL), breast cancer (BRCA), glioblastoma (GBM), endometrial carcinoma of uterus (UCEC), etc (Figure 1C). Then the differential expression of GDI2 was verified between two normal hepatic cell lines (L02, WRL-68) and seven hepatoma cell lines (Huh7, SK-HEP1, BEL-7402, PLC/PRF/5, SMMC-7721, HepG2, Hep3B) in vitro. Moreover, which transcript exerts effect in these HCC cell lines was also confirmed by designing different primers amplifying certain fragments for different transcripts of GDI2 gene. Both qRT-PCR and western blot assays confirmed that compared with normal hepatic cells L02, GDI2 expression was increased in most of hepatoma cells (Figure 1D, 1E). Although there are three transcripts GDI2 mRNA, Transcript I (NM_001115156.2) and Transcript II (NM_001494.4) are two confirmed mature transcripts, while Transcript III (XM_017016071.2) is a predictive transcript overlapped with Transcript I only with a latter Start Codon (Figure 1F). Correspondingly, PCR amplification displayed that almost both Transcript I and II of GDI2 gene contributed to GDI2 expression in HCC cell lines (Figure 1G). Finally, based on the cut-off criteria (|log2-fold change (FC)|>1, adjusted P-value<0.05) for median value of GDI2 expression, a total of 1225 Differentially Expressed Genes (DEGs) were identified after the analyses of TCGA RNA-seq data between GDI2-high and -low groups, including 654 upregulated and 571 downregulated DEGs illustrated by Volcano Map (Figure 1H), 10 of which were specifically displayed by Heat Map (Figure 1I).

Functional enrichment analysis associated with GDI2 by GO, KEGG and GSEA

To better understand the functional implications of GDI2 in HCC, GO and KEGG functional enrichment analyses were performed based on 1225 DEGs between high- and low-GDI2 expression, which indicated that the GDI2-associated genes engaged in 342 GO terms of biological processes (BP), 62 terms of cellular components (CC) and 74 terms of molecular function (MF). Thereinto, massive enrichment of Receptor ligand activity (GO:0048018), Hormone activity (GO:0005179), Metal ion transmembrane transporter activity (GO:0046873) and Extracellular matrix structural constituent (GO:0005201) categorized by MF, Stress response to metal ion (GO:0097501), Extracellular structure organization (GO:0043062), Cell-cell adhesion via plasma-membrane adhesion molecules (GO:0098742) and Second-messenger-mediated signaling (GO:0019932) categorized by BP, Haptoglobin-hemoglobin complex (GO:0031838), Transmembrane transporter complex (GO:1902495), Protein-lipid complex (GO:0032994) and Collagen-containing extracellular matrix (GO:0062023) categorized by CC, were highly associated with the aberrant expression of GDI2 (**P < 0.01; Figure 2A).

As for advanced analysis, gene set enrichment analysis (GSEA) was applied to confirm the key signaling pathways referring to GDI2 expression data sets. Significant differences with normalized enrichment score |NES| > 1.0, adjusted P-value < 0.05, false discovery rate (FDR) < 0.25 in enrichment of the MSigDB Collection h.all.v7.0.symbols [Hallmarks] and c2.cp.v7.0.symbols [curated] were observed in a total of 39 and 690 pathways, respectively. In particular, GDI2 was related to the Fatty acid metabolism, Angiogenesis, G2M checkpoint, Epithelial mesenchymal transition of Hallmarks symbols (adj.P = 0.008, FDR = 0.002; Figure 2B–2E), and similarly, to the Fatty acid metabolism by REACTOME, Cell extracellular matrix interactions by REACTOME, P53 regulation pathway by PID, Cell adhesion molecules CAMS by KEGG of C2_curated symbols (adj. P = 0.017, FDR = 0.009; Figure 2F–2I).

Furthermore, protein-protein interaction (PPI) enrichment analysis was performed to predict co-regulatory protein network of GDI2, and the functional interactions among proteins were illustrated by STRING online database. The interactions with a combined score > 0.7 were considered statistically significant, covering 620 protein network interactions. The resultant protein network illustrated by STRING encompassed the Red Node GDI2 with interactive proteins as multicolor Bubbles, and interaction Edge indicating both functional and physical protein associations (Figure 2J).

Immune infiltration analysis associated with GDI2

The enrichment analyses indicated that the occurrence of HCC was strongly associated with tumor microenvironment (TME) like extracellular matrix (ECM) organization. Therefore, the relative tumor infiltration levels of a list of 509 genes were quantified and the abundance of a diverse set of 24 tumor-infiltrating adaptive and innate immune cell types was predicted in individual tissue sample by single-sample gene set enrichment analysis (ssGSEA) analysis (Figure 3A). Then, Spearman correlation analysis between GDI2 expression and immune cell infiltration level displayed that the GDI2 expression was positively correlated with the abundance of T helper cells (R = 0.287; P < 0.001; Figure 3B), Th2 cells (R = 0.225; P < 0.001; Figure 3C) and Tcm cells (R = 0.271; P < 0.001; Figure 3D); while negatively correlated with Cytotoxic cells (R = -0.290; P < 0.001; Figure 3E), Dendritic cells (DCs) (R = 0.-280; P < 0.001; Figure 3F), Plasmacytoid Dendritic cells (pDCs) (R = -0.291; P < 0.001; Figure 3G), Th17 cells (R = -0.196; P < 0.001; Figure 3H), B cells (R = -0.106; P = 0.041; Figure 3I), and Neutrophils (R = -0.113; P = 0.030; Figure 3J). While other immune cell subsets, including T cells, Treg cells, NK cells, and macrophages were weakly correlated with GDI2 expression (Figure 3A).

Association analyses of GDI2 expression with clinicopathologic variables

To explore the association between GDI2 expression with 11 kinds of clinicopathologic characteristics, GDI2 expression was firstly classified into high- and low-level, and Chi-square test was applied to explore the difference between high and low expression of GDI2 in different clinicopathological characteristics of HCC patients (Table 1). Here, there was significant differences between high and low expressions of GDI2 in HCC patients with clinical index as weight (P = 0.035), BMI (P = 0.020) and AFP (ng/ml) (P = 0.019). Then Logistics regression was applied to explore the clinical risky factors in HCC patients with GDI2 expression on both TPM values and high- and low-classifications (Table 2). It was displayed that high- and low-expression of GDI2 was only associated with histologic grade (OR=1.58(1.03-2.44); P = 0.035) (Table 2; upper); while consecutive TPM values of GDI2 expression were associated with patients under more progressive stage (OR= 1.01(1.00-1.01) for T1 vs. T2/T3/T4, P < 0.001), more advanced pathologic stage (OR=1.01(1.00-1.01) for Stage I vs. Stage II/III/IV, P = 0.001), more serious histologic grade (OR=1.01(1.00-1.01) for G1/G2 vs. G3/G4, P = 0.007), and mutated TP53 status (OR=1.01(1.00-1.01) for WT vs. MUT, P < 0.001) (Table 2; lower). Though the odds ratios were critical (only 1.01(1.00-1.01)), it was indicated that the associations between GDI2 and these four clinical factors were significant (*P < 0.05) but not strongly dependent. Finally, the Wilcoxon rank sum test was further applied and demonstrated that the GDI2 expression was significantly associated with clinicopathologic features as T stage (T1-2 vs. T3-4; P = 0.004), Pathologic stage (PI-II vs. PIII-IV; P = 0.001), Histologic grade (G1-2 vs. G3-4; P = 0.004), AFP (ng/ml) level (AFP<=400 vs. >400; P = 0.009), TP53 status (WT vs. Mut; P = 0.007) and tumor status (Tumor free vs. With tumor; P = 0.047) (Figure 4A–4F); and consistently associated with personal physical index as Age (age<=60 vs. >60; P = 0.040), Weight (weight<=70 vs. >70; P = 0.006) and BMI (BMI<=25 vs. >25; P = 0.002) (Figure 4G–4I). Those association analyses suggested that GDI2 expression was probably an independent risky factor for HCC patients, and significantly associated with clinicopathologic characteristics as T stage, pathologic stage, histologic grade and TP53 status.

Prognostic analysis and model construction for GDI2 expression in HCC patients

To explore the risky factors associated with patients’ prognosis, the Univariate Cox regression analysis was performed to figure out the hazard clinicopathologic indicators for HCC survival based on high- and low-GDI2 expression, from which the significant variables (P < 0.1) as T stage (P < 0.001, hazard ratio (95% confidence interval) [HR (95%CI)] = 2.109 (1.469-3.028)), M stage (P = 0.018, HR (95%CI) = 4.032 (1.267-12.831)), pathologic stage (P < 0.001, HR (95%CI) = 2.074 (1.418-3.032)), tumor status (P < 0.001, HR (95%CI) = 2.361 (1.620-3.441)), and GDI2 expression (P < 0.001, HR (95%CI) = 1.844 (1.298-2.620)), were put into the Multivariate Cox regression analysis for further investigation. The results confirmed that tumor status (P = 0.003, HR (95%CI) = 2.145 (1.300-3.539)) and GDI2 (P = 0.014, HR (95%CI) = 1.836 (1.130-2.983)) were independent prognostic factors for HCC patients’ survival (*P < 0.05) (Table 3).

To validate whether the model construction was effective, receiver operating characteristic (ROC) curve was displayed to measure the discrimination value of GDI2. The calculated area under curve (AUC) of GDI2 was 0.748, indicating that GDI2 owned an efficient ability to discriminate hepatic carcinoma from normal liver and might be a potential diagnostic biomarker (Figure 5A). Then to provide clinicians with a quantitative approach to predicting the prognosis of HCC patients, a nomogram integrating GDI2 and tumor status was constructed based on multivariate Cox analysis. We compared the predictive accuracy of this nomogram with that of GDI2 and tumor status, obtaining the nomogram performance (C-index: 0.599 (0.573-0.626)) between predicted values 0.5 and 1.0. In addition, HCC patients with tumor (100 points) and high GDI2 level (77 points) could receive a total point score of 177 in this nomogram. Then the probabilities of 1-, 3-, 5-year survival were determined by drawing a vertical line from the total point axis at a value of 177 straight downward to the outcome axis, respectively (Figure 5B).

Next, the Kaplan-Meier survival analysis was performed in 184 HCC patients with GDI2-high expression and 186 cases with GDI2-low expression to evaluate the relationship between GDI2 expression and survival status of HCC patients in TCGA cohort. It showed that high GDI2 expression was strongly associated with a worse overall survival (OS) (P < 0.001, HR (95%CI) = 1.84 (1.30-2.62); Figure 5C), a worse progression-free interval (PFI) (P=0.027, HR (95%CI) = 1.40 (1.04-1.87); Figure 5D), as well as a worse disease-specific survival (DSS) (P = 0.036, HR (95%CI) = 1.61 (1.03-2.51); Figure 5E) than that of low GDI2 expression. In survival prediction nomogram for OS, the calibration curve conformed well to observations in all patients, integrating GDI2 and tumor status index, as well as T stage, pathologic stage and histologic grade with a Hosmer-Lemeshow test, implied no departure from perfect fit (Figure 5F). These validations suggested that the models based on GDI2 expression were effective for predicting short-term or long-term survival in HCC patients.

As for the subtypes of clinicopathological features associated with prognosis, the K-M survival analysis for sub-population displayed that GDI2 expression significantly impacted the OS rate in HCC patients with certain clinical index as: Age>60 (HR=2.08(1.31-3.28); P = 0.002; Figure 6A), BMI>25 (HR=2.10(1.21-3.66); P = 0.009; Figure 6B), and AFP(ng/ml)<=400 (HR=2.01(1.21-3.35); P = 0.007; Figure 6C). Meanwhile, HCC patients with low GDI2 expression could improve their prognosis when they were with clinicopathological subtypes as: T1 Stage (N=49, P = 0.010; HR (95%CI) = 2.18 (1.21-3.94); Figure 6D), negative lymph nodes (N=252, P < 0.001; HR (95%CI) = 2.252 (1.427-3.555); Figure 6E), without metastasis (N=266, P = 0.004; HR (95%CI) = 1.919 (1.232-2.990); Figure 6F), with fibrosis ishak score-1/2&3/4&5/6 (HR=2.56(1.24-5.25); P=0.011; Figure 6G); early pathologic stage I&II (HR=1.66(1.03-2.68); P=0.037; Figure 6H), and no vascular invasion (HR=2.23(1.32-3.76); P=0.003; Figure 6I). Then significant subtypes with prognostic impact were together intuitively illustrated by Forest Map (Figure 6J). These results suggested that GDI2 could be applied as a prognostic indicator in early stage of HCC.

Data collection for bioinformatics analysis from TCGA data repository

We downloaded level 3 HTSeq - FPKM formatted RNAseq data and clinical information from LIHC-hepatocellular carcinoma project from the website: https://portal.gdc.cancer.gov/ of The Cancer Genome Atlas (TCGA). A total of 371 cases with gene expression data and clinical information were collected by discarding those RNAseq data without clinical information [39]. Meanwhile, 50 normal control cases from GTEx database were downloaded from UCSC XENA (https://xenabrowser.net/datapages/) and unified handled by Toil procedure [40]. Level 3 RNAseq data formatted as HTSeq-FPKM (Fregments Per Kilobase per Million) from TCGA and GTEx were transformed into TPM (transcripts per million reads) for subsequent analyses. This study was in accordance with the publication guidelines provided by TCGA (https://cancergenome.nih.gov/publications/publicationguidelines). All data used in this study were obtained from TCGA without containing any human participants or animals performed by any of the authors, and hence ethics approval and informed consent were not required.

Analysis of immune infiltration characteristics by ssGSEA

The ssGSEA (single-sample Gene Set Enrichment Analysis) method classifies marker gene sets in a single sample with common biological functions, chromosomal localization, and physiological regulation [41]. Normalized HCC gene expression profiles from formatted TPM data of single sample were compared with the immunocyte signatures using GSVA (R package) [42]. We quantified the relative tumor infiltration levels of immunocyte signature genes, including a total of 509 genes predicting the abundance of 24 tumor-infiltrating adaptive and innate immune cell types in individual tissue sample [43].

The following 24 types of immune cells were obtained: B cells; T cells, Helper T cells (Th), Cytotoxic T cells (Tc), CD4+ T cells, CD8+ T cells, type-1 T helper cells (Th1), type-2 T helper cells (Th2), type-17 T helper cells (Th17), Regulatory T cells (Treg), gamma delta T cells (γδT), central memory T cells (Tcm), effector memory T cells(Tem), follicular helper T cells(Tfh); Dendritic cells (DCs), activated Dendritic cells (aDCs), immature Dendritic cells (iDCs), plasmacytoid Dendritic cells (pDCs); natural killer (NK) cells, CD56 bright natural killer cellsr (CD56+NK), CD56 dim natural killer cells (CD56-NK); eosinophils, mast cells, neutrophils and macrophages. The correlation between GDI2 and these immune cells was analyzed by Spearman correlation, and Wilcoxon rank sum test was adopted to explore the association of the infiltration levels of immune cells between the high- and low- expression groups of GDI2 gene.

DEGs analysis between high and low GDI2 expression groups

We firstly used Wilcoxon Rank Sum Test to compare the expression of GDI2 gene between tumor and normal specimens from TCGA combined with GTEx database in pan-cancers, referring to 33 kinds of solid tumors. Then in LIHC, GDI2 expression was compared in both paired and unpaired normal-versus-tumor tissues. Finally, the GDI2expression profiles (HTSeq-counts) were compared between the divided high and low GDI2 expression groups to identify differentially expressed genes (DEGs) by using DESeq2 R package [44]. |log2-fold change (FC)|>1 and adjusted P-value<0.05 were considered as threshold values for the DEGs.

GO, KEGG and GSEA enrichment analysis

Gene Ontology (GO) functional analysis, including cellular component (CC), molecular function (MF), and biological process (BP), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were performed on the DEGs based on high and low GDI2 expression levels by using the ClusterProfiler package (http://www.bioconductor.org/) [45] with P value adjusted by Benjamini and Hochberg method, from which adjusted P-value of <0.05 was considered as significance for GO function and KEGG pathway.

Gene Set Enrichment Analysis (GSEA; http://software.broadinstitute.org/gsea/index.jsp) is a computational method that determines whether prior defined functions or pathway sets of genes show statistical significance, concordant differences between two biological states [46]. Thus in this study, GSEA was applied to generate an ordered list of all genes according to their correlation with GDI2 expression, then carried out by R package to elucidate the significant differences observed between high- and low- GDI2 expression. The datasets h.all.v7.0.symbols.gmt [Hallmarks] (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#H) and c2.cp.v7.0.symbols.gmt [Curated] (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C2) from MSigDB Collections were chosen as reference gene sets, and the expression level of GDI2 was regarded as a phenotype. Gene set permutations were performed 1000 times for each analysis. A function or pathway term with adjusted P-value <0.05 and false discovery rate (FDR) <0.25 was considered to be statistically significant enrichment.

Protein-protein interaction (PPI) network

Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 11.0) online database was used to predict PPI network co-regulated by GDI2 gene and exploring the functional interactions between co-regulatory proteins [47]. An interaction with a combined score > 0.7 was considered statistically significant.

Clinical association and prognostic analysis

Wilcoxon rank sum test and t. test were utilized to evaluate the clinical index (age, gender, weight, BMI, AFP(ng/ml), etc.) in non-paired samples and paired samples with different expression levels of GDI2 (the Fisher exact test was used when needed) [48]. The association of GDI2 expression with clinicopathological features, such as TNM stage (TNM stage is a way of staging tumors, in which T represents the range of primary tumors, N represents the presence and extent of regional lymph node metastasis and its scope, M represents the presence or absence of a distant transfer), pathologic stage, histologic grade, and TP53 status, was evaluated by Chi-square test, and Logistics regression was performed on both continuous TPM values and high- and low-classifications of GDI2 expression.

Clinical prognosis correlating GDI2 to Overall Survival (OS), Progression-Free Interval (PFI) and Disease-Specific Survival (DSS) of HCC patients, was assessed using Survminer package [49]. Survival curves were constructed using the Kaplan-Meier method, and the differences between the survival curves were examined by the log-rank test [50]. Clinicopathological risky factors were evaluated by Cox regression analyses. Univariate Cox proportional hazards regressions were applied to estimate the individual hazard ratio (HR) with GDI2 gene expression for HCC survival. The significant variables from the univariate analyses (P < 0.1) were then put into the multivariate analysis [51]. Multivariate Cox analysis was used to compare the influence of GDI2 expression on survival along with other clinical characteristics in order to find independent variables. The HR with 95% confidence interval (CI) was measured to estimate the hazard risk of individual factors. Prognostic analyses were performed in both main-group and sub-group of clinicopathological indicators with GDI2 expression, respectively.

Nomogram evaluation

Calibration and discrimination are the most commonly used methods for evaluating the performance of models [52]. In this study, the Calibration curves were graphically assessed by mapping the nomogram-predicted probabilities against the observed rates, and the 45° line represented the best predictive values. A concordance index (C-index) was used to determine the discrimination of the nomogram, and it was calculated by a bootstrap approach with 1000 resamples [53]. The predictive accuracies of the nomogram and separate prognostic factors were compared using both the C-index and ROC (receiver operating characteristic) analysis. As the frequently-used method for binary assessment, ROC analysis was performed by pROC package [54] to assess the effectiveness of the transcriptional expression of GDI2. The computed area under the curve (AUC) value ranging from 0.5 to 1.0 indicates the discrimination ability from 50 to 100%.

Cell lines and cell culture

The human normal hepatic cell line L02 and hepatic embryonic cell line WRL-68 were the kind gifts from Infection Department of First People’s Hospital of Yunnan Province. Seven human hepatoma cell lines: Huh7, SK-HEP1, BEL-7402, PLC/PRF/5, SMMC-7721, HepG2, Hep3B were purchased from Tongpai Biotechnology Company (Shanghai, China). All the cell lines were cultured in DMEM/high-glucose medium (Gibco, USA) supplemented with 10% fetal bovine serum (Gibco, USA) and 1% Penicillin-Streptomycin, and performed STR authentication [55].

RNA isolation and qRT-PCR analysis

Total RNA from HCC cell lines was isolated using RNAprep FastPure Animal/Cell Total RNA Extraction Kit (TSINGKE, Beijing, China) according to manufacturer's protocols. Then approximately 2μg of RNA were reverse transcribed using the Goldenstar™ RT6 cDNA Synthesis kit (TSINGKE, Beijing, China). For GDI2 detection, qRT-PCR labeled by SYBR green master mix (TSINGKE, Beijing, China) was performed on a LightCycler 480 system (Roche, USA). To confirm the expression from different transcripts of GDI2, primers indicating different transcripts and common primers of GDI2 for qRT-PCR are listed as below: GDI2: Forward: 5’ATTCCACAGAACCAAGTCAATCGA 3’; Reverse: 5’CTTCTCAGGCTCCTTGGTTTCC 3’; GDI2-Transcript I: Forward: 5'GGGCACCGGCCTGACGGAATGTA 3'; Reverse: 5'TGCCAGGGCTTCTGCTTCAGTGG 3'; GDI2-Transcript II: Forward: 5'CCACCCGAGTCAATGGGGAGAGG 3'; Reverse: 5'CTTGATGGGGTGGCTGAGGATGC 3'. ACTB was used as an internal reference: Forward: 5’CACCATTGGCAATGAGCGGTTCA 3’; Reverse: 5’AGGTCTTTGCGGATGTCCACGT 3’.

Western blot

Cultured cells were collected and washed twice with ice-cold PBS, then lysed in RIPA lysis mixed with PSMF (Beyotime, China) buffer on ice for 30 minutes. Cell lysates were collected immediately according to the manufacturer’s instructions and subjected to bicinchoninic acid (BCA) protein assay (Beyotime, China) for concentration determination. Equal amounts of proteins in each lane were separated by SDS-PAGE and subsequently transferred onto polyvinylidene fluoride (PVDF) membranes (Millipore, USA). The membranes were blocked with 3% (w/v) skim milk and then incubated with specific primary antibody at dilution of 1:5000 (anti-GDI2; Invitrogen, USA) overnight followed by incubation with horseradish peroxidase (HRP)-conjugated secondary antibody at dilution of 1:10000 (anti-rabbit-HRP; Cell Signaling Technology, USA). β-actin was used as an internal control at dilution of 1:10000 (anti-beta Actin; Cell Signaling Technology, USA). An enhanced chemiluminescence (ECL) chemiluminescence kit (ABclonal, China) was used to detect immunoreactive protein bands using a Gel Doc XR imaging system (Bio-Rad, USA).

Statistical analysis

All statistical analysis and plots referring to phenotype and expression profiles in 371 HCC patients from TCGA-LIHC were conducted using R package (v.3.5.1) (http://cran.r-project.org/web/packages/rms/index.html). In gene expression analysis, the median GDI2 expression was regarded as the cut-off value. All hypothetical tests were two-tailed and all reported P values < 0.05 were considered significant, marked as *P< 0.05; **P < 0.01. The statistical analysis was carried out using Graphpad Prism 8.0 and SPSS version 22.0.