Prognostic analysis of tumor mutation burden and immune infiltration in hepatocellular carcinoma based on TCGA data

In order to explore the prognosis of tumor mutation burden (TMB) and the relationship with tumor infiltrating immune cells in hepatocellular carcinoma (HCC), we downloaded somatic mutation data and transcriptome profiles of 376 HCC patients from The Cancer Genome Atlas (TCGA) cohort. We divided the samples into high-TMB and low-TMB groups. A higher TMB level indicated improved overall survival (OS) and was associated with early pathological stages. One hundred and nine differentially expressed genes (DEGs) were identified in HCC. Moreover, based on four hub TMB-related signatures, we constructed a TMB Prognostic model (TMBPM) that possessed good predictive value with area under curve (AUC) of 0.701. HCC patients with higher TMBPM scores showed worse OS outcomes (p < 0.0001). Moreover, DCs subsets not only revealed higher infiltrating abundance in the high-TMB group, but also correlated with worse OS and hazard risk for high-TMB patients in HCC. Meanwhile, CD8+ T cells and B cells were associated with improved survival outcomes. In sum, high TMB indicates good prognosis for HCC and promotes HCC immune infiltration. Hence, DCs and the four hub TMB-related signatures can be used for predicting the prognosis in HCC as supplements to TMB.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the most common pathological subtype of liver cancer; it is the sixth most common type of cancer and the fourth leading cause of cancer-related death in the world [1][2][3]. Although the main treatment of early HCC is surgery, 50% of the patients are at an advanced stage at the time of diagnosis [4]. The 5-year overall survival (OS) rate for HCC is about 12.5% [3,5]. Targeted therapy has improved the OS in HCC; however, the overall efficacy is unsatisfactory. The emergence of immunotherapy has identified new therapeutic prospects for HCC [6], especially immune checkpoint inhibitors (ICIs), such as programmed death 1 (PD-1)/ programmed death ligand-1 (PD-L1) and cytotoxic Tlymphocyte associated antigen-4 (CTLA-4) [7,8]. In recent clinical studies, both nivolumab and AGING pembrolizumab have exhibited better prognosis as second-line therapy in advanced HCC after sorafenib treatment [9,10].
Although ICIs therapy has been shown to be effective in several types of malignancies [11][12][13][14][15], it has shown varying effects in various types of cancer, even in the same cancer patients [15]. Hence, the identification of accurate biomarkers is an urgent need for screening patients who can benefit from immunotherapy and to monitor the prognosis of immunotherapy. Recently, an increasing number of biomarkers have been identified for predicting the efficacy of immunotherapy, including DNA damage repair (DDR) [16], microsatellite instability (MSI) [17,18], neoantigens [19], and HLA presentation of neoantigens against tumor [20][21][22]. TMB is a novel biomarker that is calculated as genetic variations per million bases of the encoded genome [23,24]. Patients with a high TMB have a superior objective response rate (ORR) and prolonged OS [19,25,26] than those with a low TMB. ORR differences can be explained by TMB in about 55% types of tumor [20]; however, TMB is not a single biomarker for predicting the efficacy of immunotherapy that may be inconsistent with specific genetic mutations in high-TMB patients. For example, patients with EGFR mutations and ALK rearrangements (EGFR+/LK+), JAK1 mutations, or JAK2 mutations are associated with low response to immunotherapy [27][28][29][30]. Meanwhile patients with KRAS mutations (KRAS+) had a better ICIs response rate [27], and STK11 mutation was found to be the main factor for PD-1 inhibitor resistance in lung adenocarcinoma with KRAS+ [31]. Moreover, tumor microenvironment (TME) has been well identified as a molecular determinant in many cancers. TUMEH [32] found that T cell infiltration in TME is closely related to the efficacy of immunotherapy. Moreover, higher infiltrating abundance of CD8+ T cell and memory activated CD4+ T cell subsets revealed prolonged OS in the high-TMB group [33]. Therefore, combined information from TMB levels, gene mutations, and immune infiltration density can be used as a novel biomarker for predicting the efficacy of immunotherapy in HCC.
The prognostic role of TMB and the relationship between TMB and immune infiltration varied for different types of cancers [33,34], and limited studies have focused on TMB with immune infiltration in HCC. Thus, we investigated the prognostic role of TMB and the potential association between immune infiltration and hub TMB-related signature in HCC, using TCGA HCC cohort and Gene Expression Omnibus (GEO) datasets. We found that high TMB was a good prognostic predictor for HCC, and that DCs and the four hub TMB-related signatures could also be used for predicting the prognosis in HCC.

TMB was related to overall survival and clinical stage
We calculated the number of TMB per million bases for 363 samples and classified them into high-TMB and low-TMB groups using the median value as the threshold (Supplementary Table 1). Kaplan-Meier survival indicated that higher TMB was associated with better OS (p = 0.0004) ( Figure 2A). However, high TMB was not in accordance with better disease-free interval (DFI) in this research ( Figure 2B). Moreover, we found that higher TMB was also associated with tumor stage (p = 0.035; Figure 2C), pathologic stage (p = 0.020; Figure 2D), and T stage (p = 0.027; Figure 2E). The TMB levels tended to decrease with tumor progression ( Figure 2C-2E), and clinical research with a larger sample size is required to verify this result.

Analysis of differentially expressed genes between the 2 TMB groups
The heatmap and volcano plot visualized that 109 differentially expressed genes (DEGs) were identified with limma software with |logFC|> 1.5 and false discovery (FDR) < 0.05 (Figure 3A, 3B and  Supplementary Table 2). The Gene Ontology (GO) enrichment analysis demonstrated that in molecular function group, these DEGs were mainly involved in extracellular matrix structural constituent, glycosaminoglycan binding, and extracellular matrix structural constituent conferring tensile strength. In the biological process group, extracellular structure  -4C and Supplementary Table 3). Thereafter, we conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and found that TMB-related signatures were involved in several pathways, including fat digestion and absorption, cholesterol metabolism, and peroxisome proliferatorsactivated receptor (PPAR) signaling pathway ( Figure  3C and Table 2). We then selected the top gene set enrichment analysis (GSEA) results of TMB-related items, wherein dilated cardiomyopathy, ECM-receptor interaction, focal adhesion, oocyte meiosis and small cell lung cancer were associated with TMB levels ( Figure 3D). Moreover, we obtained protein-protein AGING interaction (PPI) networks with the STRING tool, and 36 proteins were related ( Figure 4D).

Constructing a risk scoring model with differentially expressed genes
Further, we identified 30 prognostic signatures associated with TMB using univariate Cox regression model from the above 109 DEGs (Supplementary Table  4). In addition, we used multivariate Cox analysis to select four independent risk signatures with p < 0.05 and acquired the coefficients (ßi) of the respective signature ( Figure 5A and Table 3). The selection process was visualized on a Venn plot ( Figure 5A). Among these four TMB-related genes, we found that the lectin galactoside-binding soluble 3 (LGALS3) AGING expression was upregulated in high-TMB groups than in low-TMB groups, and the expression of Nuclear Pore Complex Interacting Protein Family Member B15 (NPIPB15), Formimidoyltransferase Cyclodeaminase (FTCD), and Decorin (DCN) was negatively correlated with the TMB level ( Figure 5B-5F). The hazard ratio (HR) with 95% confidence interval was shown in the forest plot ( Figure 5G). Using the multivariate Cox regression model, TMBPM was constructed with the following formula: TMBPM = 0.179 x LGALS3 − 0.096 x NPIPB15 − 0.145 x FTCD − 0.126 x DCN. The multivariate model indicated that a high risk score was associated with poor survival (P < 0.001) ( Figure 6A, 6B, 6D). Moreover, we classified the HCC patients into high-risk (n = 175) and low-risk (n = 175) groups using the median value as the cutoff value ( Figure 6C, 6D). The receiver operating characteristic (ROC) curve for 5year OS prediction suggested that the model possessed predictive accuracy with AUC = 0.701 ( Figure 6C), and Kaplan-Meier plot showed that patients with high TMBPM revealed worse OS than those with low TMBPM ( Figure 6D). Furthermore, HCC prognosis model nomograms were constructed based on different clinical characteristics ( Figure 6E). We calibrated the 1year, 3-year, and 5-year OS predictions for HCC patients, and all the calibrated curves were well-fitted AGING ( Figure 6F-6H). Although TMBPM can accurately predict the 3-year and 5-year survival rate in HCC patients, whether the TMBPM maintained the independent predictive value needs to be investigated and validates on larger samples.

Comparison of immune cell abundance between high-TMB and low-TMB groups
Based on the newly developed CIBERSORT software [35], we intended to compare the differential profiles of immune fractions between the high-TMB and low-TMB groups. After filtering out patients with P > 0.05 with the "CIBERSORT" package, we obtained fractions of 22 immune cells in 56 HCC patients, and the results were displayed in the box ( Figure 7A), where different colors represented various cell subsets.
Meanwhile, we revealed the differential abundance of immune cells between low-TMB and high-TMB groups with heatmap plot ( Figure 7B), wherein we could intuitively find that CD8+ T cell, M0 macrophage, and M2 macrophage formed the majority of the components. Moreover, the Wilcoxon rank-sum test indicated that the infiltration levels of dendritic cells resting (P = 0.001), eosinophils (p < 0.001), and T cells regulatory (Tregs) (p = 0.02) were higher in the high-TMB group than in the low-TMB group ( Figure 7C). In addition, the density of neutrophils (p = 0.006) showed a lower infiltrating level in the high-TMB group. In accordance with previous mutation analysis and Kaplan-Meier analysis, lower TMB commonly inhibited the immune infiltration levels in HCC patients, contrary to clear cell renal cell carcinoma (ccRCC) [36]. Macrophages were identified using CIBERSORT in HCC. Macrophages were significantly enriched in tumors and displayed similar proportions in the immune fractions of high-and low-TMB groups. Although M2 was a dominant Macophages, there were not clearly distinguished in M0, M1 and M2 macrophages, consistent with other previous reports [37,38]. AGING Table 2. Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis for the differential genes. We intended to compare the differential profiles of immune fractions between the high TMB and low TMB groups; therefore, we further intended to investigate the potential prognosis of immune cells based on the Tumor Immune Estimation Resource (TIMER) database.
Furthermore we investigate 28 subpopulations of immune cells based on the GSEA database (https://tcia.at/home) ( Figure 9A and Table 5

DCs infiltrates in HCC indicated survival outcomes
The  Table 5) which are significantly associated with bad prognosis in other cancer [39] were not hazard risk (Table 5) for HCC in our study, and we also don't observed separation of the MDSCs and Tregs subpopulations related to immune suppression from the subpopulations related to the effector function (activated T cells, Tcm, Tem CD4+ and CD8+ cells) ( Figure 9B, 9C). The separation was found in other cancers including lung squamous cell carcinoma, ovarian cancer, pancreatic cancer, and melanoma [39].  The progression was also characterized by distinct immune cell patterns between high-and low-TMB group based on GSEA database. Within each TMB group, the composition of the immune cells was divergent in HCC. For example, activated DCs, EOS and Neu were enriched in stage I and stage II tumors and only depleted in stage III in high-TMB patients ( Figure 9D), and that activated CD4 and Mem_B were only depleted in late stage tumors and enriched in early stage tumors ( Figure 9D) in low-TMB samples. There were not any pots in stage IV for high-TMB, because there were not enough samples for analysis ( Figure 9D). Meanwhile, Univariate Cox analysis showed that higher levels of mDCs (HR = 21.823, P = 0.013) and Imm_DC (HR = 10.151, p = 0.022) with poor OS rates (Tables 4, 5 and Figure 8C, 8I) were hazard factors in High TMB group, while higher density of pDCs (HR = 0.302, p = 0.021) with improved OS rates ( Figure 8J, P < 0.01) was marginally protective infiltrating cells in low-TMB groups. These results might explain why some patients with high TMB are not responsive to therapy with checkpoint blockers and also why some patients with low TMB are responders.

AGING
In summary, these analyses showed that both the genomic profiles and the specific tissue context contribute to the cellular composition of the immune infiltrates. Furthermore, whether using CIBERSORT, TIMER database or GSEA database, it is noteworthy that inflated levels of DCs and prognosis of DCs were different between high-TMB and low-TMB groups. mDCs and immature DCs subsets showing higher infiltrating abundance in the high TMB group was associated with lowered OS and higher risk factor, indicating that these high TMB HCC patients with higher infiltrates of mDCs and immature DCs had poor survival outcomes, meanwhile the higher fraction of pDCs was lower risk factor and was associated with better OS in low TMB samples. Accordingly, Kaplan-Meier analysis and Wilcoxon rank-sum test showed that dendritic cells can also be used as secondary prognosis indicators for immunotherapy in HCC.

DISCUSSION
Liver cancer (LC) remains the most common cause of cancer-related death and the second major cause of cancer-related death worldwide, despite considerable improvements in its treatment [2]. Since 2017, immune checkpoint inhibitors (ICIs) have showed promising results in the immunotherapy for advanced HCC [9]. However, few HCC patients are able to obtain benefits from this treatment. Thus, many researches have focused on the identification of predictive biomarkers for immunotherapy.
In the current research, we found that mutations in TP53 (28%), CTNNB1 (24%), and TTN (25%) are frequently found in HCC (Figure 1). Wang found that HCC patients with TP53 mutation were significantly correlated with high TMB (P = 0.0005) and exhibited poor prognosis (OS: HR = 1.58, P = 0.0109) [40], indicating that HCC patients with TP53 mutation were more likely to benefit from immune treatment. Furthermore, CTNNB1 mutation was positively correlated with TMB-H and TP53 in HCC [41]. A metaanalysis found that NSCLC patients with EGFR mutation could not benefit from immunotherapy, and EGFR-positive patients had low TMB [42,43].
TMB is a novel prognosis biomarker for ICIs therapy in breast cancer [44] and other tumors [45,46]. Robert found that high somatic TMB patients treated with ICIs exhibited better OS [45]. Hellmann demonstrated that small-cell lung cancer patients with high TMB treated with either nivolumab plus ipilimumab or nivolumab monotherapy had better prognosis, irrespective of the level of PD-L1 expression [25]. In this study (Figure 2), HCC patients with high TMB had significantly better OS (p = 0.0014) than those with other malignancies [33]; moreover, higher TMB may induce immune recognition and prognosis improvement. Further, TMB was negatively correlated with the tumor stages (P = 0.035) and AJCC-T stages (P = 0.027) ( Figure 2E). These results suggest that the TMB level declined with AGING tumor invasion and progression; this prompted us to investigate the potential relationship among TMB, DEGs, and immune infiltrates for identifying more prognostic biomarkers for HCC.
Subsequently, four hub TMB-related signatures were identified ( Figure 5 and Table 3) on univariate and multivariate Cox analysis (positive correlation: LGALS3, negative correlation: NPIPB15, FTCD and DCN). Furthermore, a prognostic model (TMBPM) was developed using four hub TMB-related signatures that can be very useful for survival prediction ( Figure 6 and Table 3). To our knowledge, this is the first TMB prognostic model to predict survival outcomes in HCC.

AGING
Patients with a high TMBPM score had worse survival outcomes. Moreover, the AUC of the ROC curve was 0.701, indicating the prediction accuracy of the TMBPM. Although we calibrated the curves of 1-year, 3-year, and 5-year prognosis prediction models, further research on larger samples is needed before clinical application. In bladder cancer (BLCA) [33], tumor mutation burden related signature (TMBRS) mode were constructed by eight hub TMB-related signatures for prognosis prediction. Although this method can  accurately predict the 3-year and 5-year survival rate in HCC and BLCA, whether the method maintained the independent predictive value needs to be investigated and validates on larger samples in different cancers.

AGING
LGALS3 (Galectin-3) is mainly involved in cell growth, cell adhesion, cell differentiation, and tumor progression and metastasis owing to its action of binding to glycoproteins. Kada found that decreased expression of Galectin-3 in gastric cancer indicated poor prognosis [47]. The core proteoglycan (DCN), a main component of ECM, negatively regulates Tregs [48][49][50][51][52]. DCN can indirectly inhibit the formation of foxp3+ Tregs via the inhibition of the TGF-β signaling pathway [53,54]. Many studies have demonstrated that DCN overexpression inhibits the progress of various tumors, such as breast cancer and colon cancer [35,55,56].
Low expression of FTCD is correlated with poor prognosis (P < 0.001) in HCC as per the TCGA data [57]; meanwhile, FTCD overexpression suppressed the proliferation of BEL-7402 and SNU499 cells, resulting in increased PTEN protein and decreased PI3K, total Akt, and phosphorylated Akt protein in BEL-7402 and SNU499 cells. As per a recent study, the PTEN-PI3K/AKT signal transduction pathway was involved in tumor immune escape via the regulation of PD-L1 expression [58]. PD-L1 played an important role in inducing specific T cell apoptosis and tumor immune escape. Thus, FTCD overexpression inhibits tumor progression and tumor immune escape in HCC via the suppression of PI3K/AKT pathway activation. Therefore, FTCD could also be a promising biomarker and a potential target for HCC treatment.
TIMER is the first method for performing integrative analysis of tumor immune cell, clinical, and genomics data [59]. TIMER database was used for assessing the relationships of 9 TMB-related signature with immune infiltration levels in clear cell renal cell carcinoma (ccRCC) [36], the 9 signatures were associated with lower immune infiltrates. After constructing the TMB prognostic model, we also compared the abundance of immune cells utilizing TIMER database, CIBERSORT software, and GSEA database, and found that the prognosis of DCs infiltration was different between high-and low-TMB group. At present, DCs are considered as prognostic indicators in cancers, because the higher infiltrated DCs are associated with better prognosis [60,61]. Cai et al. found that high infiltration of DC in hepatoma indicates a higher disease-free survival time [60]. Single-cell sequencing of 16,498 HCC cells found that, compared with primary hepatocellular carcinoma, more DCs and CD8+ T cells were infiltrated in early recurrent tumors, while fewer regulatory T cells played an immunosuppressive role [37] in early recurrent tumors. Due to the high affinity of PD-L1 binding to CD80 on the DC surface, CD80 was preferentially bound to PD-L1. Thereby the competitive inhibition that CD80-CD28 mediated COstimulation of DC on CD8+ T cells blocked antigen presentation and inhibited the activation of CD8+ T cells. This suggests that the mechanism of immune escape in early recurrent tumors is different from that in primary hepatocellular carcinoma. This will contribute to the development of more effective therapeutic targets and biomarkers for immunotherapy in HCC patients [37]. According to the origin, DCs are divided into DCs derived from myeloid DC (mDC) and lymphatic dendritic cells (LDC) or plasmacytoid dendritic cells (pDC) [62]. In this study, the prognosis of different DC subgroups was different between high-TMB and low-TMB groups. For example, higher fraction of pDCs associated with improved OS rates (Table 5 and Figure  8J) was marginally protective infiltrating cells and immune response in low-TMB group, while higher fraction of pDCs was a hazard risk for high-TMB patients (Table 5) AGING immature pDCs in the tumor microenvironment promoted the tumor [63,64]. These previous research results might explain that higher infiltration of mature pDCs was favorable factor while immature pDCs was risk factor in HCC. On the other hand, although the fraction of pDCs was not different between the two TMB groups (Figure 9A), pDCs is might mainly mature pDCS in low TMB while immature pDCS in high TMB in HCC. However, this hypothesis needs to be verified by flow cytometry, and/or immunohistochemistry.  In this study, we found that high-TMB patients had higher infiltration levels of resting dendritic cells (DCs) (Figure 7), poor OS (Figure 8), and that higher infiltration levels of mDCs (Table 4), especially the immature DCs (Table 5) was a hazard factor in HCC based on TIMER database and GSEA database. Myeloid-DCs (mDCs) can induce both primary immunosuppression and tumorigenesis [65]. mDCs are considered to be a negative factor in anti-tumor immunity due to their decreased expression and function [66]. Further, in the present study, the infiltration level density of T cells regulatory (Tregs) was higher in the high-TMB group. Previous studies have shown that Tregs could inhibit the proliferation of CD4+ T cells and the secretion of IL-2 [67,68]. SATO found that the improvement in tumor-specific CD8+ T cells and reduction in Tregs cells can effectively improve the prognosis in cancer patients [69]. Empirical researches have shown that mDCs are related to the selection of thymic Tregs, differentiation, proliferation, and functional regulation of peripheral lymphoid tissues [70,71].
Although higher macrophage level was a hazard risk factor associated with poor survival outcomes, which is in agreement with that of breast cancer [72], no significantly different macrophage levels were observed in the two TMB groups. Higher MDSCs infiltrations was associated with improved OS, but we couldn't cleanly distinguish the difference of MDSCs infiltrations levels and risk factor of MDSCs between AGING high-and low -TMB groups in HCC. Therefore, we could not easily determine the potential relationship between TMB, infiltration of macrophage and MDSCs.
In our study, infiltrating levels of DCs, TMB, and TMB-related hub signatures can be used as prognostic indicator in HCC. However, TMB should be analyzed together with DCs infiltrating levels, especially in stage III and stage IV patients, wherein TMB is generally lower. Moreover, the four TMB-related signatures can also be prognosis index in HCC. While there are some limitations should be taken into consideration: (1) although the role of TMB in the prognosis of HCC and its potential relationship with immune cell infiltration had been identified by CIBERSORT software, TIMER database and ssGSEA, different analysis methods also have conflicting conclusions, which need to be verified through a large number of clinical data or clinical trials in the future; (2) The four TMB-related genes and immune cell infiltrates are also needed to be validated based on basic experiment in the future.

CONCLUSIONS
In conclusion, higher TMB correlated with improved survival outcomes and might induce local immune recognition in HCC. Patients were divided in two risk groups using TMBPM based on 4 hub TMB-related signature, and TMBPM index may serve as a promising prognostic biomarker for HCC in the future. Different DC subgroups (mDCS, pDCs and immDCs) and the infiltrating levels of these DC subgroups were different risk factors and also were associated with different survival outcomes between high-and low-TMB samples, while may play an important role in metastasis, suppression and progress of HCC.

Acquisition of multi-omics data resource
All the data used in this article were obtained from the TCGA database on the GDC website (https://portal.gdc.cancer.gov/), including HCC mutation data; transcriptome profiles; and clinical information, such as age, sex, AJCC-TNM stages, pathological stages, tumor stages, and survival outcomes. The statistical results of somatic mutation were visualized with the maftools software. The transcriptome profiles with HTSeq-FPKM Format of 371 patients were also downloaded from the TCGA database with the GDC tool.

Correlation analysis between TMB and prognosis
First, we filtered out the germline mutations annotated in the dbsnp and ExAC databases. Then, we defined and calculated the TMB of each sample as the total amount of coding variants/the length of exons (38 million), where the detected variants were considered as frameshift deletion mutation, in-frame deletion, frameshift insertion mutation, in-frame insertion, missense mutation, nonsense mutation, nonstop mutation and silent. As per the median TMB value, the patients were divided into the high-TMB and low-TMB groups. Kaplan-Meier survival analysis with log-rank test was performed to assess the differential survival rate between the high-TMB and low-TMB groups. In order to improve the accuracy of Kaplan-Meier survival analysis, we excluded these samples with survival duration of <10 d. Moreover, the association of TMB distribution with several clinical variables (AJCC-TNM stage, tumor stage, and case stage) was evaluated using Kruskal-Wallis test.

Differential expressed genes and functional analysis in the 2 TMB groups
For the high-TMB and low-TMB groups, we used the limma software to identify DEGs in the two groups with |logFC| > 1.5 and FDR < 0.05. The heatmap was displayed with the "heatmap" package to visualize the differences. ClusterProfiler package was implemented for GO and KEGG analysis with q value < 0.05. Furthermore, GSEA analysis was performed with FDR < 0.3 based on JAVA8 platform using the TMB level as the phenotype. Moreover, STRINGdb software was used for protein-protein interaction analysis.

Identification of prognostic genes
Thirty TMB-related genes were obtained from 109 DEGs using univariate Cox survival analysis with P < 0.05. Subsequently, four hub TMB-related signatures were identified using the stepwise regression screening method; the process has been shown in the form of a Venn diagram. TMBPM was used to calculate TMBPM = Ʃ (βi × EXPi) based on multivariate Cox analysis. ROC curve was utilized to assess the predictive score of four TMB-related signatures in HCC.
The nomogram of the HCC prognosis model was established via univariate and multivariate analyses, combined with clinical characteristics. Then, calibration curves were constructed for the prediction of 1-year, 3year, and 5-year survival in HCC.

Assessment of immune-infiltrating cells and prognostic analysis
CIBERSORT software (https://cibersort.stanford.edu) was used to evaluate the compositions of the immune cells in HCC samples, with P > 0.05 to improve accuracy of the estimated results, Heatmap has been AGING used to show the distribution of immune cell fractions, and violin plot visualizes the differential distributions of cells with two TMB levels.
Based on the TIMER database, multivariate Cox analysis that was fitted by function coxph (βi) from R package "survival" was used to identify the immune infiltration cells. Subsequently, we calculated the HR with 95% confidence interval (95%CI). Further, Kaplan-Meier survival analysis was performed with a P value < 0.05 of log-rank test to show the differential survival outcomes between different levels of immune infiltrates in HCC.
Subsequently, we identified 28 kinds of immune cells that are over-represented in the tumor microenvironment using single sample gene set enrichment analysis (ssGSEA) [39]. All immune cell types were considered enriched in a patient or group of patients when FDR (qvalue) ≤ 0.1 and normalized enrichment score (NES) >0. The similarity of the enrichment of immune infiltrates (averaged NES) were calculated using multidimensional scaling (MDS). The distribution of selected cell types for individual patients were analyzed with t-distributed stochastic neighbor embedding (t-SNE) using the Matlab toolbox t-SNE.

Statistical analyses
All the statistical analyses were conducted using R software (Version 3.5.2). The p-values were adjusted for multiple testing based on FDR according to the Benjamini-Hochberg approach, P value < 0.05 was regarded statistically significant. Differential analysis and normalization were mainly conducted by "limma" package of the R software (version 3.5.2). The Kaplan-Meier analysis with log-rank test or Cox regression model was performed by "survival" package. Student's t test was used for continuous variables, while categorical variables were compared by χ2 test. The non-parametric two-sided Wilcoxon-rank sum test was utilized for comparing two groups and the Kruskal-Wallis test was suitable when it comes to two or more groups. Differences and correlations among immune cells were analyzed with the vioplot (https://cran.rproject.org/web/packages/vioplot/index.html).

Editorial note
& This corresponding author has a verified history of publications using a personal email address for correspondence.