A prognostic signature established based on genes related to tumor microenvironment for patients with hepatocellular carcinoma

Data acquisition and preprocessing

We collected the transcriptomic data and clinical features of 374 HCC patients and 50 control samples from The Cancer Genome Atlas (TCGA) database by visiting Genomic Data Commons (https://portal.gdc.cancer.gov/) [13]. Four patients without complete clinical follow-up and microdissection data were eliminated, while the remaining samples were included in the analyses. Moreover, for external validation, the GSE76427 dataset, which used the GPL10558 platform and included 115 primary HCC tumors and 52 adjacent nontumor tissues, were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Screening and analysis of differentially expressed TRGs (DE-TRGs)

We identified 4,061 TRGs from previously published studies [14, 15]. First, “limma” package was employed to screen DE-TRGs between HCC and adjacent nontumor tissues [16]. The data in fpkm format were preprocessed by normalization and filtering under the criteria of the absolute value of FDR <0.05 and |log2Fold Change|>1.

NMF clustering

Prognosis-related DE-TRGs were screened using Univariate CRA (UCRA) under P < 0.01 as the threshold. Next, NMF clustering analysis was performed on the expression matrix of prognosis-related DE-TRGs in the “NMF” R package [17]. The correlation coefficients were calculated, the inner feature structure of gene expression matrices was predicted, and 50 iterations were recycled during the clustering.

Analysis of immune cell infiltration

Microenvironment cell populations-counter (MCP-counter) was applied to score the abundance of immune cells such as CD8+ T cells, NK and endothelial cells, myeloid DCs, neutrophils, cytotoxic and B lymphocytes, monocyte-derived cells, and fibroblasts based on the gene expression [18] using the “MCPcounter” package.

Developing and validating the RS model

We divided 370 HCC patients from TCGA into the training (n = 262) and validation (n = 108) sets at the ratio of 7:3 ratio without control (Table 1). Notably, the GSE76427 dataset contained both tumor and adjacent non-tumor samples, only the tumor samples were used for developing a RS model. Next, prognosis-related genes were mined by performing UCRA on DE-TRGs, the model was refined by building a penalized feature using the “glmnet” package [19], and multivariate CRA (MCRA) was used to screen the characteristic genes for the RS model. Finally, we calculated the coefficient value of each gene and the riskscore as follows:

$$Riskscore={\displaystyle \sum _{i=1}^n(coefi\times Xi)}$$

coef represented the Cox regression coefficient for each gene, and X referred to gene expression.

The samples were clustered by their median RS value into the low-risk group (LRG) and high-risk group (HRG). The prognostic significance of the model was analyzed by Kaplan-Meier (KM) analysis with two-sided log-rank test using the “survminer” package. Area under the receiver operating characteristic (ROC) curve (AUC) in “timeROC” [20] package was used to reflect the prediction accuracy of the model.

We also created a nomogram by integrating clinicopathological features such as gender, age, TNM stages, grade, clinical stage, and the RS model using the “restricted mean survival (rms)” package. The reliability of the nomogram was tested according to the ROC curve, calibration curve (CC, bootstrap 1000), decision curve analysis (DCA), and C-index.

Analysis of biological functions

Annotated gene sets ‘c2.cp.kegg.v7.4’ and ‘c5.go.v7.4.symbols’ obtained from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb) were used for GSEA (software version 4.1.0).

Statistical analysis

Data were analyzed using the R (version 4.1.1) package. Volcano maps and violin plots were visualized by the “ggplot2” package [21] and “ggpubr” package, respectively. CRA was performed using the “survival” package. By employing the chi-squared or Fisher’s exact test, differences between both sets and the correlation between clinical characteristics and the RS were analyzed. KM survival curves with the log-rank tests were generated using the “survminer” package. The ROC and calibration curves were plotted using the “timeROC” and “rms” packages, respectively [22]. The restricted mean survival (RMS) package was used for computing the C-index of the model. A p < 0.05 indicated a significant difference.

Data availability statement

The dataset analyzed in this study is available in (GSE76427) at (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76427).

NMF consensus clustering analysis classified three molecular HCC subtypes based on TRGs

A total of 1,717 TRGs were defined as DE-TRGs from 4061 TRGs between 370 (FDR <0.05 and |log2FC|>1, Figure 1A), and UCRA identified 704 prognosis-related TRGs. NMF consensus clustering analysis was performed based on the expression profiles of 704 prognosis-related TRGs. The cophenetic correlation coefficient was significantly reduced at k = 3, which was therefore selected as the optimal number of clusters (Figure 1B–1D). Survival analysis showed that the OS and progression-free survival (PFS) of patients in C2 and C3 was significantly worse compared to C1 (Figure 1E, 1F). Moreover, immune scores of myeloid dendritic cells, B lineage, endothelial cells, cytotoxic lymphocytes, neutrophils, fibroblasts, monocytic lineage, NK cells, and CD8+ T cells were different across C1, C2, and C3 (Figure 2A–2I).

Wound healing (Immune C1), IFN-gamma dominant (Immune C2), inflammatory (Immune C3), lymphocyte depleted (Immune C4), C5 (Immunologically quiet), and TGF-beta dominant (Immune C6) are the six immune subtypes of solid tumors [23]. Hence, the distribution of these immune subtypes in three molecular clusters was analyzed. Specifically, Immune C3 and C4 were the primary immune subtypes in C1, the proportions of Immune C4 and C3 were the highest in C2 and C3, respectively, and Immune C1 was only detected in C2 and C3 (Figure 1G).

Constructing and validating the RS model containing 13 DE-TRGs

UCRA and LASSO regression analysis was conducted using the DE-TRGs identified from the training set. Independent variables gradually increased to zero with lower λ (Figure 3A). In addition, 32 genes were mined using the partial likelihood deviance method (Figure 3B). Finally, the RS model was constructed with the 13 TMGs identified by MCRA as follows:

Risk score = (−0.5974 × CAPN3 expression level) + (0.3376 × HAVCR1 expression level) + (0.4161 × GRIN2D expression level) + (0.6536 × BMI1 expression level) + (0.4733 × SLC30A3 expression level) + (0.2293 × GLP1R expression level) + (0.2211 × CSAG3 expression level) + (−0.5401 × NR4A3 expression level) + (−0.8607 × INPP4B expression level) + (0.1339 × ETV4 expression level) + (−0.2406 × MYH4 expression level) + (−0.4394 × ICA1 expression level) + (0.4485 × ST6GALNAC4 expression level.

Next, the AUC value was calculated to assess the RS model in all cohorts. The AUC values for 1-, 3-, and 5-year OS of patients in the training set were 0.835, 0.854, and 0.911, respectively (Figure 3C). In the entire TCGA-LIHC cohort, the AUC values for 1, 3, and 5-year OS of patients were 0.790, 0.785, and 0.830, respectively (Figure 3D). In the TCGA validation cohort (the validation set), the AUC values for 1-, 3-, and 5-year OS of patients were 0.704, 0.614, and 0.552, respectively (Figure 3E). Though the AUC values of the 1- (0.561), 3- (0.614), and 5-year (0.692) OS of patients in the GSE76427 cohort were <0.7, the RS model was still accurate to some extent, especially in predicting the 5-year OS. Such a result might be related to the short median duration of follow-up of patients in this cohort (Figure 3F).

Patients from all HCC cohorts were assigned with a RS and divided into HRG and LRG based on the median RS. The KM survival curves showed that in the training set, the prognosis in the LRG was significantly better than in HRG (Figure 4A). The results from survival analyses on patients in the validation set, the entire TCGA, and the GSE76427 cohort were the same results. In the HRG, the patient’s OS was significantly poor than in LRG (Figure 4B–4D). Additionally, in the TCGA-LIHC cohort, the RS varied in patients with different clinicopathological characteristics and it was positively correlated with high pathological grades and advanced T, TNM stages (Figure 4I).

Construction of a nomogram by integrating RS and clinicopathological features

We performed UCRA and MCRA to determine the correlation between the OS and clinicopathological features such as age, tumor grade, and stage, gender, T and M stage, as well as RS. The results demonstrated that only the RS (HR: 1.026, 95% CI: 1.014–1.037, P < 0.001) could independently predict risk (Table 2). A nomogram integrating clinicopathological features and RS was developed (Figure 4E) to predict 1-, 3-, and 5-year survival. For example, the clinical features of a 50-year-old female with LIHC were T2NxM0 stage II, grade 1, and high RS. Then, final total score of the variables was 399. The AUC values of patients’ 1-, 3-, and 5-year survival rates were 0.767, 0.561, and 0.322, respectively. The CC of the nomogram showed a high consistency between the observed 1-, 3- and 5-year survival rates and the predicted values in the training sets (Figure 4F). In addition, compared to the other clinicopathological characteristics, the AUC value of the RS was higher (Figure 4G). DCA showed that the RS and nomogram could effectively predict the OS compared to all other clinicopathological features (Figure 4H).

Comparing the TRG-RS model with other models

We compared our 13-gene signature RS model to other models, including 8-gene signature related to inflammation response (a TME signature) [24], 3 immune-related prognostic genes [14], and 3 tumor doubling time-related immune genes [15]. Firstly, we calculated the RS for all patients using MCRA and plotted the ROC curve using their corresponding genes. Then, the patients were classified by the median RS value into HRG and LRG. A significant survival difference of patients in the two groups was observed. In addition, the prognosis of patients in LRH was better than in HRG (Figure 5A–5D). As shown by the ROC curve, the AUC value of our model in predicting average 1- year prognosis was the highest. The AUC values of our model for predicting 1-, 3-, and 5-year survival were higher compared to the other three models (Figure 5E–5H). The C-index of our model was 0.734, which was higher than other models (Figure 5I). Finally, we used the RMS time to assess the prediction accuracy of our model at different time points, and observed that our model outperformed the other three models when predicting survival time longer than 60 months. These results indicated that our model was more effective and accurate in predicting 5-year OS and survival timer longer than 5 years (Figure 5J).

Functions of our RS model

Based on the normalized enrichment score and adjusted p-value (q-value), GSEA was applied to identify enriched pathways in the two groups. The complement and coagulation pathways, fatty acid and cytochrome p450 drug metabolism was enriched by genes in the LRG (Figure 6A). Moreover, GO terms such as the process of cellular amino acid and alpha-amino acid catabolism and the T cell receptor complex were enriched genes in the LRG (Figure 6B). The correlation between the RS and target genes was analyzed according to the expression of genes related to classical cancer-related pathways including fatty acid metabolism, DNA replication, and cell cycle. DNA replication-related genes (MCM6, MSH2, and MSH6) were positively correlated with the RS (Figure 6C). Additionally, immune cell infiltration in all patients from TCGA-LIHC was analyzed. The RS showed a positive relationship with CD8+ T cells, monocytic lineage, and myeloid DCs (Figure 6D). The correlation between the RS/tumor mutation burden (TMB) and infiltrating immune cells was analyzed using an MCP counter, and the results are presented in Figure 6E.

Survival prediction for patients with different clinicopathological characteristics using the RS

We predicted the survival of patients with different clinical traits using the RS model. Compared to those in HRG, the KM survival showed better prognostic outcomes of patients with age >65 and ≤65, pathological grade (G1-2 and G3-4), gender, stage I–II, and III–IV in the LRG (Figure 7).