In this study, The Cancer Genome Atlas (TCGA, http://cancerge.nome.nih.gov/) (Tomczak et al., 2015) and the International Cancer Genome Consortium (ICGC, www.icgc.org) (Zhang et al., 2019) were used to download somatic mutation data (MAF files) of TCGA-LIHC cohort and the LIRI-JP cohort. The primary objective of the ICGC database was to provide a comprehensive elucidation of genome changes in multiple cancers that result in human disease burden. Among the ICGC database, the tumor data from different cancer types (or subtypes) were collected, including abnormal gene expression, somatic mutation, epigenetic modification, and clinical data. The ICGC database contains 25,000 tumor genomes. Meanwhile, the clinicopathologic characteristics and the prognostic information of the patients in TCGA-LIHC cohort, such as gender, age, and clinical stage, were obtained from the UCSC Xena website (http://xena.ucsc.edu/) (Goldman et al., 2019). Moreover, RNA sequencing data (count value), containing mutation data and survival data of 353 patient samples (TCGA-LIHC), were downloaded from TCGA database for subsequent analysis and were annotated by the annotation file of the GRCh38 version from the Ensembl database (http://ftp.ensembl.org/pub/current_gtf) (Howe et al., 2021). In addition, the copy number variations data were obtained from TCGA database. The clinical characteristics of patients are listed in Table 1.

With the development of tumor genomics, the mutation annotation format (MAF) is being widely accepted and used to store detected somatic mutations. In this study, the maftools package (Mayakonda et al., 2018) and the GenVisR package (Skidmore et al., 2016) were utilized to visualize the somatic mutation data downloaded from TCGA; meanwhile, the GenVisR package was also used to visualize the somatic mutation data obtained from ICGC. Moreover, the mutation of DOCK2 was revealed by the G3viz package (Guo et al., 2019). Additionally, to evaluate whether the genes have copy number variation in liver cancer, GISTIC2.0 in the Genepattern (https://cloud.genepattern.org/) cloud analysis platform was used to analyze the copy number variation data of liver cancer in TCGA database (M. et al., 2006).

The liver cancer patients were divided into mutation group (DOCK2-MUT) and wild group (DOCK2-WT) according to the DOCK2 mutation status of the gene expression data downloaded from TCGA. A survival analysis was performed based on the DOCK2 mutation and prognosis information of liver cancer patients, thus investigating the prognostic difference between the DOCK2 mutation group and the wild group. Moreover, the patients’ data obtained from TCGA was randomly divided into a training set (N = 264) and a testing set (N = 89) at a ratio of 7:3. The DOCK2 mutation prediction model was conducted using the random forest (RF) method (Yperman et al., 2019) in the training set, and the model performance was quantified via the receiver operating characteristic (ROC) curve.

The prognostic model was built in light of the gene expression data of 28 DOCK2 mutant liver cancer patients with clinical information. First, univariate Cox proportional hazard regression analysis was performed to initially identify overall survival (OS)-related genes (p-value<0.05). Next, RF and multivariate Cox regression analyses were conducted to construct a prognostic model. The formula for calculating the risk score is risk score = exp gene 1 × β gene 1 + exp gene 2 × β gene 2 + exp gene 3 × β gene 3 + … exp gene n × β gene n (exp gene n indicates the expression level of gene n; β gene n indicates the regression coefficient of gene n calculated by multivariate Cox regression). Moreover, correlation analysis was performed between the DOCK2 mRNA expression and the risk score, as well as between the DOCK2 mRNA expression and the characteristic genes mRNA expression in the model.

The liver cancer patients in the DOCK2 mutant group with clinical information were divided into high-risk groups and low-risk groups in light of the median risk score. The OS analysis was performed using the Kaplan–Meier (KM) survival curve and time-dependent ROC, thus evaluating the prediction accuracy of the model. Then, the univariate Cox regression analysis and the multivariate Cox regression analysis were conducted in light of the age, gender, clinical stage, tumor stage, and risk score in DOCK2 mutant liver cancer patients. Meanwhile, the risk score and clinical characteristics were analyzed using correlation analysis.

Given that different DOCK2 mutation types may have different effects on tumorigenesis, the expression data of liver cancer patients were divided into two subgroups: inactivated mutation subgroup and other non-silent mutation subgroups. The two subgroups were assessed via the KM survival curve and time-dependent ROC.

Tumor mutation burden (TMB) refers to the total number of somatic mutations in the exon coding region of the genome that have substitutions, insertions, or deletions per Mb base in a tumor sample. The TMB score of each liver cancer sample is the total number of somatic mutations (including non-synonymous point mutations, insertions, and deletions in the coding region of exons)/target region size, and the unit is mutations/Mb (Chan et al., 2019). Microsatellite (MS) is defined as a short tandem repeat (STR) in the human genome including single-nucleotide repeats, dinucleotide repeats, and even more nucleotide repeats; microsatellite instability (MSI) refers to the change of any length of microsatellites due to the insertion or deletion of repeat units in tumor tissues compared to normal tissues (Hile et al., 2013). MSI is calculated as the number of insertions or deletions in gene repeats. In this study, the relationship between the risk score and TMB and the correlation between the risk score and MSI were analyzed, respectively.

In order to assess the impact of gene expression value on the DOCK2 mutant type compared to the DOCK2 wild type, the limma R package (Ritchie et al., 2015) was used to conduct the discrepant analysis between the DOCK2 mutant group and DOCK2 wild-type group. The absolute value of log fold change (logFC) > 0.5 and p-value <0.05 were set as the threshold for differentially expressed genes. Among them, the genes with logFC > 0.5 and p-value <0.05 were considered upregulated differential genes, while the genes with logFC < −0.5 and p-value <0.05 were regarded as downregulated differential genes, and the aforementioned results were displayed with heat maps and volcano maps.

Gene Ontology (GO) enrichment analysis is a common method for large-scale functional enrichment studies of genes in different dimensions and levels and generally includes three aspects: biological process (BP), molecular function (MF), and cellular component (CC) (Ashburner et al., 2000). The Kyoto Encyclopedia of Genes and Genomes (KEGG) (Ogata et al., 1999) is a widely applied database that stores numerous data about genomes, biological pathways, diseases, and drugs. Additionally, the clusterProfiler R package (Yu et al., 2012) was applied to identify significantly enriched biological processes and pathways by GO functional annotation and KEGG biological pathway enrichment analysis. A p-value <0.05 was considered statistically significant.

Gene Set Enrichment Analysis (GSEA) is a calculation method to assure whether a set of predefined genes show statistical differences between two biological states, generally applied to estimate changes in the pathway and bioprocess activity in sample expression datasets (Subramanian et al., 2005). Based on the gene expression profile data of DOCK2 mutant group and DOCK2 wild-type group patients in TCGA-LIHC dataset and the reference gene sets “c5.go.v7.4.entrez.gmt” and “c2.cp.kegg.v7.4.entrez.gmt” downloaded from the MSigDB database (Liberzon et al., 2015), the GSEA method included in the clusterProfiler R package was used to conduct enrichment analysis of TCGA-LIHC gene expression profile data, thus studying the differences in the biological processes of genes between the DOCK2 mutant group and DOCK2 wild group. A p-value <0.05 was considered statistically significant.

Gene Set Variation Analysis (GSVA) (Hnzelmann et al., 2013; Liberzon et al., 2015), a nonparametric unsupervised analysis method, is widely utilized in the evaluation of metabolic pathways enriched in different samples by converting the expression matrix of genes between different samples into the expression matrix of gene sets between samples. To study the biological process variation of the DOCK2 mutant group compared with the DOCK2 wild group, the “GSVA” R package (Hnzelmann et al., 2013; Liberzon et al., 2015) was used to perform gene set variation analysis, and the enrichment scores of each sample in each pathway in the reference gene set “h.all.v7.4.symbols.gmt” were downloaded from the MSigDB database. Moreover, the GSVA results were also analyzed for correlation with the risk scores.

The immune microenvironment is a comprehensive LoAD system, which is mainly composed of immune cells, inflammatory cells, fibroblasts, interstitial tissues, and various cytokines and chemokines. The infiltration analysis of immune cells in tissues has an important guiding role in disease research and treatment prognosis.

ESTIMATE analysis, an algorithm that quantifies the immune activity (immune infiltration level) in tumor samples on the basis of gene expression data, can reflect the richness of the gene characteristics of the matrix and immune cells. The content of stromal cells and immune cells in TCGA-LIHC was calculated by an ESTIMATE R package (Yoshihara et al., 2013). The correlation between the ESTIMATE score and the expression level of characteristic genes and DOCK2 in the prognostic model was also evaluated.

CIBERSORT is an algorithm that deconvolves the expression matrix of immune cell subtypes in light of the principle of linear support vector regression, making use of RNA-Seq data to assess the abundance of immune cells in the tissue. In this study, the proportion of 22 immune cell subtypes in TCGA-LIHC immune microenvironment was calculated by the CIBERSORT algorithm (Newman et al., 2019) in the R package. The number of permutations was 1,000, and a p-value <0.05 was considered accurate for calculating the content of immune cells. Based on Pearson correlation analysis, the correlation between the expression levels of characteristic genes and DOCK2 in the prognostic model and 22 types of immune cells in liver cancer was calculated.

To test the biological processes and cell signal transduction pathways that the characteristic genes of the prognostic model may be involved in, an immune gene set was obtained from the ImmPort database (Bhattacharya et al., 2014) (https://www.immport.org), and the relationship between the characteristic genes and DOCK2 in the prognostic model and immune genes was determined. The correlation between the expression of the HLA family and the risk score of the prognostic model was also conducted.

All data calculations and statistical analysis were performed using R programming (https://www.r-project.org/, version 3.6.3). Multiple testing corrections were determined using the Benjamini–Hochberg (BH) method, and FDR correction was conducted using multiple tests to reduce the false-positive rate. For the comparison of two groups of continuous variables, the statistical significance of normally distributed variables was estimated by independent Student’s t test, and the differences between non-normally distributed variables were analyzed by using the Mann–Whitney U test (the Wilcoxon rank-sum test). The survival R package (Durisová and Dedík, 1993) was utilized in the survival analysis. The Kaplan–Meier survival curve was used to reveal differences in survival, and the significance of the difference in survival time between the two groups of patients was assessed via a log-rank test. Univariate and multivariate Cox analyses were applied to determine independent prognostic factors. As for assessing the accuracy of the risk score to estimate the prognosis, the receiver operator characteristic (ROC) curve was drawn by the pROC package and ROCR package, and the area under the curve (AUC) was calculated (Sing et al., 2005; Robin et al., 2011). All p-values reported from statistical tests were two-sided, and a p-value <0.05 was considered statistically significant.

First, 54 genes with mutation frequencies greater than 6% in TCGA-LIHC patients obtained from TCGA were identified (Figure 1A). Moreover, the 54 genes were further confirmed using the data downloaded from the ICGC database (Figure 1B). Among them, the mutation frequency of DOCK2 was relatively high, and the DOCK2 mutation was visualized (Figures 1C, D). GISTIC 2.0 was utilized to analyze copy number variation data in TCGA, identifying obviously amplified or deleted genes, and the results showed that DOCK2 had no significant amplification or deletion (Figures 1E, F).

Survival analysis was conducted based on the DOCK2 mutation data and prognostic information of liver cancer patients, and the results revealed that the mutation of DOCK2 had an essential impact on the prognosis and survival of patients (Figure 2A). In the training set, the RF method was used to construct a DOCK2 mutation prediction model in the mRNA data (Figures 2B, C). The ROC curve was used to evaluate the performance of the model, and AUC scores close to 1 indicated that the model had high sensitivity under a very low false-positive rate. The model AUC value in the training cohort was 1.00 and that in the validation cohort was 80.4% (Figure 2D), which demonstrated that the performance of this model was sufficient to effectively predict DOCK2 mutation in other transcription cohorts.

Univariate Cox proportional hazard regression analysis was carried out in the gene expression data of 28 DOCK2 mutant LIHC patients with clinical information and 641 genes related to OS were discovered (p-value <0.05) (Figure 3A). Then, we conducted the RF method to find out the most important features connected with prognosis, and 15 genes were screened out (Figure 3B). Finally, a multivariate Cox proportional hazard regression analysis identified the five genes associated with OS, which are secretin receptor (SCTR), tetratricopeptide repeat, ankyrin repeat and coiled-coil domain-containing 1 (TANC1), Alkb homolog 7 (ALKBH7), FRAS1-related extracellular matrix 2 (FREM2), and G protein subunit gamma 4 (GNG4). Cox regression coefficients of the characteristic genes were calculated, and the risk score of each sample was defined as the sum of the expression of each characteristic gene multiplied by its regression coefficient. To assess the predictive power of the prognostic model, the risk scores of DOCK2 mutant and DOCK2 wild-type patients were calculated and ranked, the survival status of each patient was displayed on the dot chart, and the expression of characteristic genes was shown on the heat map (Figures 3C, D). Meanwhile, the correlation between DOCK2 expression and risk score and characteristic gene expression was analyzed, respectively. The expression of DOCK2 was dramatically negatively correlated with the risk score (Figure 3E). DOCK2 expression was significantly positively correlated with SCTR (r = 0.293, p-value = 9.6e-10), TANC1 (r = 0.607, p-value = 1.8e-43), and FREM2 (r = 0.252, p-value = 1.8e-07), whereas the expression of DOCK2 had a significant negative correlation with ALKBH7(r = −0.162, p-value = 0.0009) (Figure 3F).

According to the median risk score, DOCK2 mutant liver cancer patients with clinical information were divided into the high-risk group and low-risk group. The results of survival analysis showed that there was a significant difference in OS between the two risk groups of 28 DOCK2 mutant samples (Figure 4A). The 1- and 3-year AUCs on the basis of the risk score obtained by the prognostic model were 0.791 and 0.822, respectively (Figure 4B). Additionally, the correlation analysis results of the risk score and the clinical characteristics of 28 DOCK2 mutant samples revealed that there were no significant differences in risk scores, different ages, genders, clinical stages, and tumor stages (Figures 4C–F). Then, univariate Cox analysis and multivariate Cox analysis were performed based on the age, gender, clinical stage, tumor stage, and risk score of DOCK2 mutant liver cancer patients, thus building a clinical prediction model, the efficacy of which in 28 DOCK2 mutant samples was 85.8% (Figure 4G). Meanwhile, the calibration curve showed both good discrimination ability and calibration (Figure 4H).

Given that different DOCK2 mutation types may have different effects on the occurrence of liver cancer, this study further divided the gene expression data of 28 DOCK2 mutant LIHC patients into two subgroups: the inactivated mutation subgroup (n = 8, containing nonsense mutation and silent mutation) and other non-silent mutation subgroups (n = 20). Survival analysis showed that significant differences in OS were observed between the two risk groups of samples in other non-silent mutation subgroups (Figures 5A, B). The time-dependent ROC analysis showed that in the subgroup of inactivated mutations, the 1- and 3-year AUCs of the risk score were both 0.833 (Figure 5C); moreover, in other subgroups of non-silent mutations, the 1- and 3-year AUCs of risk scores were 0.651 and 0.665, respectively (Figure 5D), suggesting that the risk score could still maintain good predictive performance in subgroups with different mutation types. After acquiring the total number of mutations to obtain TMB and assessing the relationship between the risk score and the TMB, we found that there were obvious differences in TMB between samples with different risk scores (p-value < 0.05) (Figure 5E). In addition, MSI between samples with different risk scores also had a significant difference (p-value < 0.05) (Figure 5F).

To analyze the effect of gene expression values on the DOCK2 mutant samples compared with the DOCK2 wild-type samples, we conducted a limma discrepant analysis to obtain differentially expressed genes. The gene expression profile data of 28 DOCK2 mutant samples and 325 DOCK2 wild-type samples were included in TCGA-LIHC, from which 12 upregulated differential genes (p-value <0.05, logFC > 0.5) and 4 downregulated differential genes (p-value <0.05, logFC < −0.5) were screened out, and the volcanic map and heat map of the differential genes were shown in Figures 6A,B. To determine the value of the differential genes, the biological processes, the cellular components, and the molecular functions were performed. GO functional enrichment analysis was first assessed on the 16 differential genes (Figure 6C and Table 2), and the results showed that these genes were mainly enriched in biological processes such as antimicrobial humoral response, antimicrobial humoral immune response mediated by antimicrobial peptides, regulation of cardiac muscle contraction, humoral immune response, regulation of striated muscle contraction, regulation of membrane potential, cardiac muscle contraction, and skeletal muscle tissue development (Figure 6D); in cellular components including fascia adherens, transport vesicle membrane, GABA-A receptor complex, GABA receptor complex, mast cell granule, integral component of synaptic vesicle membrane, postsynaptic membrane, and dendrite membrane (Figure 6E); and in molecular functions including benzodiazepine receptor activity, secondary active monocarboxylate transmembrane transporter activity, GABA-gated chloride ion channel activity, amino acid:sodium symporter activity, oligosaccharide binding, inhibitory extracellular ligand-gated ion channel activity, peptidoglycan binding, and amino acid:cation symporter activity (Figure 6F). Then, pathways significantly affected by 16 differential genes were also performed (Figure 6G and Table 3), and the data revealed that the 16 differential genes were involved in GABAergic synapse, nicotine addiction, endometrial cancer, adherens junction, and bacterial invasion of epithelial cells (Figure 6H).

GSEA biological function enrichment analysis of DOCK2-MUT and DOCK2-WT genes was performed, and the results showed that the genes in DOCK2-MUT and DOCK2-WT were enriched in biological processes including coagulation and regulation of cytosolic calcium ion concentration (Figures 7A, B and Table 4).

Next, the results of GSEA biological pathway enrichment analysis suggested that biological pathways such as complement and coagulation cascades, glycosphingolipid biosynthesis, and ganglio series were identified among the targets in DOCK2-MUT and DOCK2-WT (Figures 7C, D and Table 4).

Furthermore, in order to comprehensively evaluate the roles of the targets in DOCK2-MUT and DOCK2-WT in liver cancer, the GSVA was conducted. The data showed three hallmarks: reactive_oxygen_species_pathway, spermatogenesis, and uv_response_dn (Figure 7E). Among them, spermatogenesis was significantly negatively correlated with risk score (p-value <0.05); uv_response_dn was obviously positively related with risk score (p-value <0.05); however, reactive_oxygen_species_pathway had no significant correlation with risk score (Figures 7F–H).

As liver cancer is considered an immunogenic tumor, the relationship between the expression of DOCK2, SCTR, TANC1, ALKBH7, FREM2, and GNG4 and the levels of immune cells and stromal cells was assessed (Figures 8A, B). The data showed a positive correlation between stromal cells and DOCK2, SCTR, TANC1, and FREM2, and a negative correlation between stromal cells and ALKBH7 and GNG4. Moreover, immune cells had a positive correlation with DOCK2 and TANC1 and a negative correlation with ALKBH7 (p-value <0.05).

In addition, the six target genes, DOCK2, SCTR, TANC1, ALKBH7, FREM2, and GNG4, were significantly correlated with specific immune-related genes. For example, DOCK2 was significantly related to the immune gene of SEMA3F; SCTR was correlated with SEMA3F and FGR; and the same situation occurred between GNG4 and SEMA3F and NFYA, FREM2 and FGR, NFYA and MPO, ALKBH7 and SEMA3F, FGR, NFYA and CALCR (p-value <0.05) (Figure 8C).

More importantly, the six target genes were obviously interrelated with the infiltration of numerous immune cells. DOCK2 gene expression was distinctly related to the infiltration of 11 immune cells; SCTR gene expression was dramatically correlated with the infiltration of 7 immune cells; TANC1 gene expression was obviously interrelated with the infiltration of 10 immune cells; FREM2 gene expression was distinctly related to the infiltration of one immune cell; ALKBH7 gene expression was dramatically correlated with the infiltration of 4 immune cells; GNG4 gene expression was markedly interrelated with the infiltration of 7 immune cells (p-value <0.05) (Figure 8D). Furthermore, the expression value of HLA-DOA was statistically significant in different risk groups (Figure 8E).