Prognostic Analysis and Biomarkers Identification of Immune Infiltration in Early and Late Stage Hepatocellular Carcinoma Based on TCGA Data

Background

Primary liver cancer is the most common cancer in the world and a common cause of cancer death.¹ Hepatocellular carcinoma (HCC) accounts for 75–85% of primary liver cancer.¹ The main risk factors for HCC are chronic infection with hepatitis B virus (HBV) or hepatitis C virus (HCV).² Currently, surgical resection and transplantation are the most effective methods for the treatment of liver cancer, but the recurrence rate is still high and the prognosis is poor.³ Surgical resection is the main intervention for HCC patients who are diagnosed in the early stages.⁴ However, early diagnosis and effective therapy of HCC remain a challenge, and most of HCC cases are diagnosed at an advanced stage, with less opportunity to take curable surgery, thus contributing to the poor survival rate.⁵ Consequently, it is urgent to reveal the underlying molecular mechanisms of HCC development from early stage to late stage.

In the current study, we downloaded mRNA, miRNA and lncRNA expression profiles of early (stage I–II) and late (stage III–IV) stage HCC patients from The Cancer Genome Atlas (TCGA) database. Then, the differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs), lncRNAs (DElncRNAs), key molecules associated with the prognosis, and important immune cell types in the late stage HCC were identified. Finally, the nomogram based on age, gender, stage, and all important factors was constructed to predict the survival of HCC.

Methods

Data Acquisition

Clinical data of 377 HCC patients were downloaded from TCGA database (collected before 25 January 2021). The mRNA and lncRNA expression profiles of 305 HCC patients, including 227 early (stage I–II) and 78 late (stage III–IV) stage HCC tissues, and miRNA expression profiles of 302 HCC patients, including 225 early and 77 late stage HCC tissues, were obtained from TCGA as well.

Differential Expression Analysis

First, the difficult-to-detect mRNAs/miRNAs/lncRNAs with a read count of 0 in more than 90% of the samples were filtered and deleted. Based on the read count of each sample, the DEmRNAs, DEmiRNAs and DElncRNAs in the late stage HCC compared to the early stage HCC were determined by edgeR with p-value <0.05. Hierarchical clustering analysis was performed using the R package “pheatmap”. To further analyze the biological functions of DEmRNAs in the late stage HCC, Gene Set Enrichment Analysis (GSEA) was performed by using GSEA (version 4.1.0). The h.all.v7.2.symbols.gmt selected as the reference gene set in the present study. P-value <0.05 was considered statistically significant.

DElncRNA-DEmiRNA-DEmRNA Regulatory Network

Given miRNAs tend to decrease the expression of their target mRNAs, we selected miRNA (down)-mRNA (up) and miRNA (up)-mRNA (down) relationships for further research. DEmiRNA-DEmRNA interaction pairs were obtained with miRwalk 2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2). Similarly, the DElncRNA-DEmiRNA interaction pairs were predicted with NPInter, Starbase and DIANA-LncBase v2. Then, the DElncRNA-DEmiRNA-DEmRNA regulatory network was constructed. Finally, Cytoscape was used to visualize the regulatory network.

Survival Analysis and Nomogram of Key Molecules in DElncRNA-DEmiRNA-DEmRNA Regulatory Network

Univariate Cox regression analysis was performed to identify the prognostic value of all biomarkers in DElncRNA-DEmiRNA-DEmRNA regulatory network. To avoid overfitting of the Cox model, lasso regression was performed on biomarkers significantly associated with prognosis. The receiver operating characteristic (ROC) curves were applied to assess the accuracy of the Cox model. Based on age, gender, stage and biomarkers significantly associated with prognosis, a nomogram was constructed using “rms” package in R. Following that, draw the calibration curve of the nomogram and observe the predictive value of the nomogram.

Immune Correlation Analysis

In order to further explore the molecular mechanism of HCC, we performed the deconvolution analysis with xCell to infer the fraction of immune cell types in HCC. xCell could estimate the abundance scores of 64 cell types based on the single sample Gene Set Enrichment Analysis (ssGSEA) method.⁶ Then, the immune cells with significant differences in the early and late stage HCC were identified. The patients were stratified into high-xCell score and low-xCell score groups based on the median value of the xCell. Univariate Cox analysis was performed in combination with xCell score of significantly different immune cells and clinical information of patients, followed by multivariate Cox analysis. Then, based on age, gender, stage and survival related immune cells, a nomogram was constructed using “rms” package in R, as well as calibration plot.

Correlation Analysis Between DEmRNAs and Survival Related Immune Cells

Pearson correlation analysis was done for DEmRNAs and survival related immune cells to determine DEmRNAs that were significantly associated with immune cells. Similar to the above process, univariate Cox analysis, Lasso regression analysis and ROC were performed. Then, age, gender, stage, immune cells and immune related DEmRNAs were included to construct the nomogram and calibration plot.

Survival Analysis and Nomograms of All Important Factors

All important factors, including key molecules in DElncRNA-DEmiRNA-DEmRNA regulatory network, survival related immune cells and DEmRNAs significantly associated with immune cells were incorporated into a nomogram.

Real-Time qPCR (RT-qPCR) Validation

Eight tissue samples of early HCC patients (n = 4) and late HCC patients (n = 4) were included in this study. We obtained the written informed consent and the approval from the ethics committee of the hospital. The Trizol reagent was used to separate total RNA from tumor tissues. FastKing cDNA first strand synthesis kit and miRNA first strand cDNA synthesis kit (tailing method) were reverse transcribed to synthesize mRNA and miRNA, respectively. SuperReal PreMix Plus (SYBR Green) and miRNA quantitative PCR kit (dye method) were used to perform RT-qPCR validation of mRNA and miRNA, respectively. The relative quantification of mRNAs and miRNAs was normalized to GAPDH and ACTB or U6 with the 2^−ΔΔCT method.

Validation and Survival Analysis of DEmRNAs in GEO

The gene expression matrix file of GSE14520 was downloaded from the GEO database and annotated according to the annotation file of the GPL3921 platform. Among them, 220 samples (early: late = 170: 50) with stage information were selected for further analysis. The expression of DEmRNAs was validated with the GSE14520 dataset. In addition, survival analysis was performed by using the survival package in R.

Statistical Analysis

R package was used to perform all the statistical analyses. R package edgeR was used for differential expression analysis of mRNA, miRNA and lncRNA. In RT-qPCR, one-way ANOVA was used to evaluate the statistical significance among different groups.

Results

Identification of DEmRNAs, DEmiRNAs and DElncRNAs

A total of 1516 DEmRNAs, 97 DEmiRNAs and 87 DElncRNAs were identified in the late stage HCC, of which 1082 mRNAs (71.4%), 92 miRNAs (94.8%) and 60 lncRNAs (69.0%) were up-regulated while others were down-regulated (Figure 1). The top 10 up- and down-regulated DEmRNAs, DEmiRNAs and DElncRNAs are displayed in Table 1. The GSEA results revealed two gene sets that were significantly up-regulated in the late stage HCC: G2M_checkpoint (p = 0.031) and mitotic_spindle (p = 0.027) (Figure 2).

Table 1 Top 10 Up- and Down-Regulated DEmRNAs, DEmiRNAs and DElncRNAs

Figure 1 The hierarchical clustering results of top 100 up- and down-regulated DEmRNAs (A), all DEmiRNAs (B) and all DElncRNAs (C) between the early and late stage HCC. Row and column represented DEmRNAs/DEmiRNAs/DElncRNAs and samples, respectively. The color scale represented the expression levels. Orange indicates above the reference channel (high expression genes). Blue indicates below the reference channel (low expression genes).

Figure 2 GSEA analysis illustrated up-regulated gene sets in the late stage HCC. (A) G2M_checkpoint; (B) mitotic_spindle.

DElncRNA-DEmiRNA-DEmRNA Regulatory Network and Survival Analysis

A DElncRNA-DEmiRNA-DEmRNA regulatory network including 78 mRNAs, 50 miRNAs and 1 lncRNA was established (Figure 3). Univariate Cox regression analysis performed on all 129 molecules in the above lncRNA-miRNA-mRNA network revealed that 50 molecules were significantly correlated with the prognosis of HCC (Figure S1A). Then, these 50 molecules were incorporated into a multivariate model. Lasso regression analysis indicated that 11 molecules were significantly correlated with the prognosis of HCC (Figure S1B and C). The ROC curve analysis indicated that the 1-, 3- and 5-year area under curve (AUC) values were 0.76, 0.72 and 0.71, respectively (Figure S1D). A nomogram was generated to predict the probability of 1-, 3- and 5-year overall survival (OS), by incorporating age, gender, stage and biomarkers significantly associated with prognosis. As shown in Figure S1E, each factor was assigned points in proportion to its risk contribution to survival. The calibration curve indicated that actual and predicted survival did not match well, and the C-index for the model was 0.559 (Figure S1F).

Figure 3 DElncRNA-DEmiRNA-DEmRNA regulatory network. The rectangle, inverted triangles and elliptical nodes indicate DElncRNAs, DEmiRNAs and DEmRNAs, respectively. Red and blue color represent up- and down-regulation, respectively.

Immunoinfiltration Analysis

With xCell, 19 significantly different immune cell types were identified between the early and late stage HCC (p-value <0.05) (Figure S2). Univariate Cox analysis performed in combination with xCell score of significantly different immune cells and clinical information of patients identified 11 survival related immune cell types. Then, multivariate Cox analysis of these 11 immune cells indicated that CD4+ Tem, Preadipocytes and Platelets were significantly correlated with survival. Kaplan–Meier curves of these 3 immune cells were presented in Figure S3A–C. The nomogram and calibration plot (C-index = 0.74) are shown in Figure S3D and E.

Correlation Analysis Between DEmRNAs and Survival Related Immune Cells

The correlation analysis between 3 immune cells and 1516 DEmRNAs was performed, showing that 82 DEmRNAs were significantly associated with Preadipocytes (|pcor| > 0.5 and p-value <0.05). Univariate Cox analysis further indicated that 70 DEmRNAs were significantly associated with survival. Then, Lasso regression analysis identified that 3 DEmRNAs (CFHR4, CYP2C9 and SAPCD2) were crucial biomarkers significantly associated with survival (Figure S4A and B). The ROC curve analysis indicated that the 1-, 3- and 5-year AUC values were 0.75, 0.74 and 0.71, respectively (Figure S4C). A nomogram was generated to predict the probability of 1-, 3- and 5-year OS, by incorporating age, gender, stage, xCell score of Preadipocytes and 3 DEmRNAs expression levels (CFHR4, CYP2C9 and SAPCD2) (Figure S4D). The calibration plot indicated that the C-index for the model was 0.73 (Figure S4E).

Survival Analysis and Nomograms of All Important Factors

Combined with the above analysis, all important factors, including xCell score of Preadipocytes and expression levels of SLC1A2, TECTB, CTNNA3, TRIM55, TRPM8, RANBP3L, RAMP1, PON1, hsa-miR-5581-3p, hsa-miR-7-5p, hsa-miR-767-5p, CFHR4, CYP2C9, and SAPCD2, were integrated into a Cox model. The nomogram is shown in Figure 4A. The calibration plot indicated that the C-index for the model was 0.758 (Figure 4B).

Figure 4 Nomogram (A) and calibration curve (B) in HCC.

RT-qPCR Validation

Two DEmRNAs (TRPM8 and TECTB) and four DEmiRNAs (hsa-miR-5581-3p, hsa-miR-7-5p, hsa-miR-767-3p and hsa-miR-490-3p) were selected to use for RT-qPCR analysis. Hsa-miR-5581-3p, hsa-miR-7-5p and hsa-miR-767-3p were up-regulated trend while TRPM8, TECTB and hsa-miR-490-3p were down-regulated trend in the late stage HCC. It is noted that the RT-qPCR results were in line with our integrated analysis (Figure 5). However, we found that RT-qPCR results lacked significance, which may be caused by the small sample size.

Figure 5 RT-qPCR results of DEmRNAs and DEmiRNAs.

Validation and Survival Analysis of DEmRNAs in GSE14520

The expression patterns of seven DEmRNAs (CFHR4, CTNNA3, CYP2C9, PON1, RAMP1, SLC1A2, and TRPM8) were validated in GSE14520 (Figure 6). The results revealed that except for CTNNA3 and RAMP1, the other five mRNAs were significantly down-regulated in the late stage HCC. Even so, CTNNA3 and RAMP1 showed the down-regulated trend in the late stage HCC. These findings were consistent with the results of our integrated analysis. In addition, survival analysis displayed that the expression of all these seven mRNAs was significantly correlated with the overall survival time of HCC patients (Figure S5).

Figure 6 Validation of DEmRNAs in GSE14520. (A) CFHR4, (B) CTNNA3, (C) CYP2C9, (D) PON1, (E) RAMP1, (F) SLC1A2, (G) TRPM8.

Discussion

At present, four different nomograms were constructed by different methods. Firstly, we constructed the DElncRNA-DEmiRNA-DEmRNA regulatory network, which included 78 mRNAs, 50 miRNAs and 1 lncRNA. Univariate Cox regression analysis and Lasso regression analysis further identified 11 key molecules associated with the prognosis, and a nomogram by incorporating age, gender, stage and 11 molecules was generated, which performed worse (C-index = 0.559). Secondly, 3 important immune cell types were identified by Univariate Cox analysis and multivariate Cox analysis to construct nomogram (C-index = 0.74). Thirdly, correlation analysis between DEmRNAs and survival related immune cells was performed, and Preadipocytes and 3 DEmRNAs expression levels (CFHR4, CYP2C9 and SAPCD2) were selected to construct nomogram (C-index = 0.73). Lastly, by incorporating all important factors identified above, a nomogram, including 1 immune cell, 11 mRNAs and 3 miRNAs, was generated to predict the survival of HCC, which achieved a good performance (C-index = 0.758) than the model constructed using only single factor.

Tripartite motif containing 55 (TRIM55) was suggested to act as a potential tumor suppressor to inhibit the malignant biological behavior of lung adenocarcinoma cells.⁷ A model of immune-related-gene-signature including TRIM55 was constructed to predict survival outcome and reflect the immune status of lung squamous cell carcinoma.⁸ It has been demonstrated that TRIM55 was decreased in HCC tissues, and its overexpression inhibited migration and invasion of HCC cells.⁹ Ding et al suggested that serum paraoxonase 1 (PON1) could potentially be used to diagnose microvascular invasion and to guide more personalized treatment strategy in HCC.¹⁰ It has been reported that low expression of PON1 was associated with poor survival in HCC patients, which may be helpful for the diagnosis of HCC.¹¹ Both TRIM55 and PON1 were detected to be decreased in the late stage HCC, suggesting their important roles in the development of HCC.

Cytochrome P2C (CYP2C) subfamily members are known to participate in clinical drug metabolism.¹² Cytochrome P450 family 2 subfamily C member 9 (CYP2C9) is a member of the CYP2C subfamily that is localized in a single gene locus on chromosome 10.¹³ Previous studies have revealed that dysregulation of CYP2C9 may participate in the development of HCC. An integrated bioinformatics analysis indicated that CYP2C9 was strongly associated with the risk of liver cancer development and progression in obese individuals.¹⁴ It has been reported that CYP2C9 is involved in the regulation of iron homeostasis and associated with the prognosis of HCC.¹⁵ Mu et al found that CYP2C9 was associated with poor OS in HCC patients.¹⁶ Suppression of CYP2C9 by hsa-miR-128-3p results in tumor cell invasion in HCC.¹² In this analysis, CYP2C9 was markedly reduced in the late stage HCC, indicating that CYP2C9 may exert momentous roles in the progression of HCC.

Aberrant expression of suppressor APC domain containing 2 (SAPCD2) has been indicated in a variety of cancers, including HCC. SAPCD2 level was correlated to the tumor size, Tumor Node Metastasis (TNM) stage and lymph node metastasis, and knocking down SAPCD2 inhibited migration and invasiveness ability of breast cancer cells.¹⁷ The expression of SAPCD2 was significantly increased in fibrosarcoma tissues, and silencing SAPCD2 inhibited the proliferation and lung metastasis of fibrosarcoma cells.¹⁸ It has been suggested that overexpression of SAPCD2 correlates with proliferation and invasion of colorectal carcinoma cells.¹⁹ Silencing SAPCD2 suppressed the proliferation, migration, and invasiveness of prostate cancer cells, suggesting that SAPCD2 may serve as an oncogene in prostate cancer.²⁰ A prognostic gene model consisting of four genes, including SAPCD2, was constructed to predict the OS of HCC patients.²¹ SAPCD2 exhibited increased expression in late stage HCC, which added to the evidence that SAPCD2 performs a critical function in HCC.

Hsa-miR-7-5p has been reported to suppress proliferation, migration and promote apoptosis in HCC cell lines by inhibiting SPC24 expression.²² It has been reported that hsa-miR-7-5p had independent prognostic significance for the OS of HCC.²³ Hsa-miR-767-3p was reported to be inhibited by hsa_circ_0000673 in the malignant progression of HCC, which evidenced the expected tumor-suppressive effect of miR-767-3p.²⁴ In addition, the up-regulation of miR-5581-3p can significantly enhance the viability and proliferation of cancer cells, and enhance the migration and invasion of cancer cells, suggesting that miR-5581-3p may play a promoting role in the progression of HCC by regulating cancer cells.²⁵ In this study, these three miRNAs were identified as significant signatures and included in the model to predict survival in HCC patients.

Studies have found that lncRNA X inactive specific transcript (XIST) participates in HCC development by targeting multiple miRNA/mRNA axes. It has been reported that lncRNA XIST accelerates HCC cell growth by inhibiting miR-488.²⁶ Wang et al suggested that lncRNA XIST up-regulates TRIM25 by targeting and binding to miR-192 to accelerate the occurrence and development of HCC.²⁷ Liu et al indicated that XIST could enhance the cell growth ability of HCC by targeting miR-200b-3p, which suggested that XIST may be a potential therapeutic target in HCC.²⁸ Similarly, hsa-miR-490-3p has been linked to HCC by different lnRNA/hsa-miR-490-3p/mRNA axes. For instance, BCYRN1/miR-490-3p/POU3F2 was suggested to promote tumor cell growth and metastasis of HCC.²⁹ SNHG15/miR-490-3p/HDAC2 axis was demonstrated to promote HCC progression.³⁰ In this analysis, XIST/hsa-miR-490-3p/VHLL axis was detected in the DElncRNA-DEmiRNA-DEmRNA regulatory network, suggesting its potential role in the progression of HCC.

However, there are some limitations in this study. Firstly, the lack of significance of RT-qPCR results may be caused by the small sample size, so it is necessary to collect samples for further study. However, the validation in GEO database increased the reliability of our analysis. Secondly, the detailed molecular mechanisms of the key DEmRNAs, DEmiRNAs and DElncRNAs identified remain unclear and require further study. Thirdly, additional prospective validation is required before the prognosis model can be considered for clinical application. The reliability of this model needs to be iteratively improved in long-term clinical applications with larger sample sizes. To further validate these findings and evaluate prognostic value and clinical applicability of the model, a multicenter prospective study with more patients and a longer follow-up period in the real world is required.

Conclusion

In summary, this study identified DEmRNAs/DEmiRNAs/DElncRNAs in the early and late stage HCC that may serve as factors related to the progression of HCC. Meanwhile, immune cell types related to the progression of HCC were obtained by immune infiltration analysis with ssGSEA. The Cox survival model was constructed to screen out the key factors significantly related to the prognosis of HCC. Finally, a Cox model with a differentiation degree of 0.758 including 14 molecules and 1 type of immune cells was obtained. This study is the first attempt to use immune cells as a variable to predict survival in HCC patients, providing new ideas for future research. In addition, the identification of new biomarkers provides potential targets for the diagnosis and management of early stage HCC.