Immune-Related lncRNA Pairs Clinical Prognosis Model Construction for Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the most common primary liver cancer in the world and the third most common cause of cancer-related death.¹ However, the traditional treatment methods for HCC, such as surgery, radiotherapy and chemotherapy, are not ideal. The high postoperative recurrence rate is also one of the main causes of death of HCC. Therefore, finding effective biomarkers for the prognosis of HCC is still a major challenge in the diagnosis and treatment of HCC.² As an inflammation-related tumor, HCC has an immunosuppressive tumor microenvironment (TME),^3,4 which can promote immune tolerance through various mechanisms.⁵ A series of progress has been made in immunotherapy.^6,7 Immunotherapy that activates tumor-specific immune response has brought new hope for the treatment of HCC.⁸

The liver is the largest immune organ in the body, rich in a variety of immune cells.⁹ HCC tumor immune microenvironment is extremely complex, which is composed of cross-linking between tumor components (such as cancer cells, stromal cells and immune cells). The imbalance of pro-inflammatory and anti-inflammatory cells in the TME promotes tumor progression.¹⁰

Long non-coding RNA (lncRNA) is an RNA molecule composed of more than 200 nucleotides,¹¹ which is closely related to the occurrence and development of tumors and can be used as a potential marker for tumor diagnosis and prognosis.¹² At the same time, a large amount of evidence shows that lncRNA is an active participant in all stages of tumor immunity.¹³ In recent years, with the continuous development of gene sequencing technology, bioinformatics databases have been widely used in the field of medical research.^14,15 Among them, the Cancer Genome Atlas (TCGA) can be used to quantitatively study and analyze the changes in expression profile during tumor occurrence and development through large-scale second-generation sequencing. The relevant research results provide a reliable basis for tumor diagnosis and prognosis evaluation. Zhou et al identified 6-lncRNAs and established a model to predict the survival rate of HCC.¹⁶ The purpose of most current studies was to build prognosis prediction models composed of several single lncRNAs. In tumor diagnosis and prognosis prediction, the combination of two biomarkers is better than one biomarker.¹⁷ Therefore, based on the TCGA database, this study constructed a new immune-related lncRNA (irlncRNA) model through a new modeling algorithm, pairing and iteration, which is not affected by the expression level. Therefore, it is of great significance to make a relatively accurate prediction of the prognosis of HCC and formulate personalized treatment strategies for HCC.

Materials and Methods

Data Download and Processing

The transcriptome profiling of human HCC was downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The relevant data were extracted by removing the incomplete data from the survival data. The gene encoding protein and lncRNA of the gene expression matrix were annotated using the human gene annotation file of Ensembl website (http://asia.ensembl.org/). The transcriptome data were divided into mRNA matrix and lncRNA matrix according to biotypes by Perl script.

The Analysis of Immune Cell Infiltration Between HCC and Normal Tissues

Using TCGA sequencing data, the CIBERSORT deconvolution algorithm was used to predict the composition and content of immune cells. R was used to draw a histogram and heatmap to show the composition of immune cells in different samples. The “corrplot” package was used to draw related heatmaps to show related immune cells. The “vioplot” package was used to draw a violin map showing the comparison of immune cells between hepatocellular carcinoma and normal tissues. The immune cell infiltration of HCC and normal tissues were analyzed by the CIBERSORT method.

Extraction and Differential Expression Analysis of irlncRNA

A dataset of identified immune genes was acquired from the ImmPort database (http://www.immport.org) and was utilized to identify irlncRNAs through a coexpression strategy. The “limma” package was used to screen the differentially expressed irlncRNA (DEirlncRNA) between normal tissues and tumor tissues based on the standard of “FDR < 0.05 and log| (fold change)| > 2”.

The Construction of DEirlncRNA Pairs

The DEirlncRNAs were paired circularly to construct a matrix of 0 or 1. Let C equal the ratio of lncRNA A to lncRNA B. If the expression level of lncRNA A was higher than that of lncRNA B, C is defined as 1; otherwise, C is defined as 0. Because some unspecified levels cannot correctly predict the survival outcome of patients, it was considered to be an effective pair whenever the number of lncRNA pairs with the expression of 0 or 1 accounts for more than 30% but less than 70% of the total number of pairs.

The Construction of the HCC Prognosis Prediction Model

The inclusion criteria of this study included HCC patients with complete survival data, and they were randomly divided into a training group and a test group. Univariate Cox regression analysis was used to evaluate the correlation between each effective lncRNA pair and survival data. LncRNA pairs with P < 0.01 were identified as candidate lncRNA pairs. Then, candidate lncRNA pairs were further analyzed by LASSO regression and stepwise multiple Cox regression analysis to select a group of lncRNA pairs, so as to construct the prognosis prediction model of HCC of irlncRNA pairs (IRLP).

The Validation and Evaluation of the IRLP Prognosis Prediction Model

The prediction performance of this model was evaluated using the time-dependent receiver operating characteristic curve (ROC) and multi-index ROC curve in the training group and test group. The survival curve and the risk curve were used to compare the 10-year survival rate of HCC patients in the high-risk and low-risk groups. Univariate and multivariate Cox regression models were used to evaluate the relationship between clinical variables, risk score and prognosis. Multivariate Cox proportional hazards regression model was used to adjust age, sex, grade and TNM stage to evaluate the independent prognostic value of the model.

Clinical Correlation Analysis of IRLP Prognosis Prediction Model

The Chi-square test was used to evaluate the relationship between the model and clinical characteristics, and Wilcoxon signed-rank test was used to calculate the difference in risk scores between different groups with these clinical characteristics.

The Correlation Analysis of Immune Cell Infiltration

Immune cell data was downloaded from the immune cell database (http://timer.comp-genomics.org). XCELL, TIMER, QUANTISEQ, McCounter, EPIC, CIBERSORT ABS and CIBERSORT were used to calculate immune cell infiltration status in samples. The difference in immune infiltrating cell content between the high-risk and low-risk groups in the model was analyzed by Wilcoxon signed-rank test. Spearman correlation analysis was used to calculate the relationship between risk score and immune infiltrating cells. The significance threshold was set to P < 0.05.

The Gene Correlation Analysis and Drug Sensitivity Analysis

R was used to obtain differentially expressed genes between normal samples and tumor samples, as well as differentially expressed genes between high-risk groups and low-risk groups in tumor samples. The Venn diagram was drawn by “Venny” package to obtain and determine the intersection of the two groups of differential gene sets. The website was then analyzed interactively using a gene expression profile (http://gepia.cancer-pku.cn/), and single-gene survival analysis was performed. The drug sensitivity was analyzed by “pRRophetic” package.

Statistical Analysis

All statistical analyses were performed with R statistical software. Univariate Cox regression analysis, LASSO regression and stepwise multiple Cox regression analysis were used to construct the IRLP prognosis prediction model. Chi-square test and Wilcoxon signed-rank test were used in the clinical correlation analysis of IRLP. Wilcoxon test was used to compare the differences between high-risk and low-risk groups. A P-value of less than 0.05 was considered statistically significant.

Results

The Analysis of Tumor Immune Cell Infiltration

The entire process flowchart of this study is shown in Figure 1. Transcriptome and clinical data were downloaded from the TCGA database. Transcriptome data included 50 normal samples and 374 tumor samples, and clinical data included 377 HCC patients. A total of 19,604 mRNAs and 14,086 lncRNAs expression profiles were isolated from the transcriptome.

Figure 1 Flowchart of this article.

Barplot was used to show the composition of immune cells between samples (Figure 2A). It could be found that the contents of macrophages M0 and T cells CD8 in the HCC group were higher than those in normal sample tissues from the heatmap (Figure 2B). The correlation heatmap showed that there was a significant positive correlation between resting NK cells and gamma delta T cells, and a significant negative correlation between macrophages M0 and T cells CD8 (Figure 2C). The CIBERSORT algorithm was used to measure the relative proportion of immune cells between tumor tissues and normal tissues, and it was represented by a violin diagram. As shown in the figure, blue represented the normal sample, and red represented the HCC sample. Correlation test showed that T cells regulatory (Tregs) and macrophages M0 increased significantly in the HCC group; B cells naive and monocytes decreased significantly (Figure 2D).

Figure 2 Immune cell infiltration analysis of tumor and normal tissues. (A), Histogram of the composition of immune cells in different samples. (B), Heatmap of the content of immune cells in different samples. (C), A correlation heatmap was generated to show the correlation of immune cells. (D), A violin map was graphed to show the difference of immune cells between the HCC group and the normal group. The blue column represents normal tissues, and the red column represents HCC tissues.

The Extraction and Differential Expression Analysis of irlncRNA

Based on the immune-related genes obtained in the immune gene database, 1700 immune-related HCC gene data were extracted by R software. At the same time, the correlation between immune-related HCC genes and lncRNAs was tested, and a total of 821 irlncRNAs were screened. A total of 153 DEirlncRNAs were screened by “limma” package as differentially expressed lncRNAs in HCC, of which 151 DEirlncRNAs were up-regulated, and 2 DEirlncRNAs were down-regulated (Figure 3A and B).

Figure 3 Identification of DEirlncRNAs. (A), Differential expression heatmap of DElncRNAs.Blue represents normal tissues, and orange represents tumor tissues. (B), Volcano plot of DEirlncRNAs. The green and red dots represent the downregulated and upregulated DEirlncRNAs, respectively.

The Construction of the HCC Prognosis Prediction Model

5223 DEirlncRNA pairs were obtained by iterative cycle and 0 or 1 matrix screening. Except for patients with incomplete survival data, a total of 365 HCC patients were included in the prognostic analysis. They were randomly divided into training group (n = 183) and test group (n = 182). Through univariate Cox regression analysis, a total of 234 DEirlncRNA pairs that were closely related to patient survival were obtained. They were included in lasso regression analysis, and the cross-validation method was used for iterative analysis. On the basis of ensuring the acceptable range of the partial likelihood of deviance, the construction of the prognosis prediction model minimized DEirlncRNA pairs for the convenience of clinical application, 9 DEirlncRNA pairs were selected (Figure 4A and B). Then, multivariate Cox analysis was carried out to screen 7 DEirlncRNA pairs related to the prognosis of HCC, and an IRLP model for predicting the prognosis of HCC was constructed (Figure 4C). Specifically, the risk score model for predicting the prognosis of HCC based on 7 irlncRNA pairs markers was shown as follows: risk score = expression value of AC099850.4|LMNTD2-AS1 * 1.0672 + expression value of ZFPM2-AS1|AL356234.2 * 0.7442 + expression value of AC004816.1|AS010280.2 * 0.7298 + expression value of TBX2-AS1|AC007773.1 * (−0.5426) + expression value of AC253536.6|AC0159123* (−0.6527) +expression value of SERTAD4-AS1|LINC02041 * (−0.7458) +expression value of AL031186.1|AC009275.1 * (−0.8633).

Figure 4 Risk assessment model construction and validation. (A), A coefficient profile plot was constructed to obtain the optimal parameter in the LASSO model. (B), LASSO coefficient profiles of the 9 DEirlncRNA pairs. (C), Forest plot of 7 DEirlncRNA pairs identified by the Cox proportional hazards regression model for multivariate analyses using a forward stepwise method. (D), The AUC value of the training group was 0.821, and the cutoff value was 2.307. (E), The AUC values were more than 0.80 in one-, two- and three-year ROC curves in the model. (F), The survival curve of the training group based on Kaplan–Meier analysis.

The Validation and Evaluation of the IRLP Prognosis Prediction Model

The IRLP prognostic model was used to analyze the time-dependent ROC of 183 sample data in the training group. The results showed that the area under the curve (AUC) of the IRLP prognostic model was 0.821, suggesting that the model had a certain ability to distinguish and predict the prognosis of patients with HCC. The corresponding value of the cutoff point at the maximum Jordan index was taken as the cutoff value of 2.307 (Figure 4D). The AUC values were more than 0.80 in one-, two- and three-year ROC curves of the model (Figure 4E). It showed that the 7 DEirlncRNA pairs in the IRLP prognostic model were highly sensitive and specific for survival prediction.

Based on the IRLP prognosis prediction model, the risk score of each case in the training group was calculated. According to the cutoff value (risk score = 2.307), the patients were divided into 64 cases into the high-risk group and 119 cases in the low-risk group. The survival curve showed that the 10-year survival rate of the high-risk group was significantly lower than that of the low-risk group (P < 0.001) (Figure 4F), indicating the effectiveness of DEirlncRNA pairs.

At the same time, the time-dependent ROC analysis and survival analysis of the test group were carried out simultaneously. The results showed that the AUC value was 0.724 (Figure 5A). The survival curve showed that the 10-year survival rate of the high-risk group was lower than that of the low-risk group, and the difference was statistically significant (P < 0.05) (Figure 5B), which confirmed the reliability of the IRLP model.

Figure 5 Evaluation of the risk model. (A), The AUC of the test group was 0.724. (B), The survival curve of the test group displaying 10-year survival rate of the high-risk group was lower. (C), Risk curve and survival outcome of each case. (D and E), Forest map by Cox univariate regression analysis and multivariate regression analysis.

The risk curve showed that the mortality rate of HCC patients in the high-risk group was higher than that in the low-risk group (Figure 5C). Patients without age, sex, grade and tumor-stage information were excluded, and Cox univariate and multivariate analysis were performed on samples of the training group to evaluate independent risk factors. In Cox univariate regression analysis, the forest map showed that there were significant differences in clinical stage [HR1.789, 95% CI (1.342–2.387), P < 0.001] and IRLP risk score [HR1.190, 95% CI (1.135–1.248), P < 0.001] (Figure 5D). In multivariate regression analysis, the IRLP risk score [HR1.162, 95% confidence interval (1.101–1.227), P < 0.001] (Figure 5E) was an independent risk factor affecting the prognosis of patients with HCC.

Clinical Correlation Analysis of IRLP Prognosis Prediction Model for HCC

The clinical correlation heatmap showed that the risk score was correlated with grade and T stage (Figure 6A). It was found that in HCC, with the increase in grade, the risk value showed an upward trend. The difference between G1, G2 and G3 was statistically significant (P < 0.05) (Figure 6B), and there was a significant difference between clinical stage I and stage III between high- and low-risk groups (P < 0.001) (Figure 6C).

Figure 6 Clinical correlation analysis of the risk model. (A–C), The risk score was correlated with grade, T stage and clinical stage in the clinical correlation heatmap and box diagrams.

The Correlation Analysis and Difference Analysis of Immune Cells

Since lncRNAs and immune-related genes were interrelated, we further explored the relationship between the model and tumor immune microenvironment. The results showed that the high-risk group in the model was significantly negatively correlated with the immune cell infiltration of Endothelial cell, Hematopoietic stem cell, Macrophage and T cell CD4 + memory resting (Figure 7A–D).

Figure 7 The correlation analysis and difference analysis of immune cells. (A–D), The high-risk group was negatively correlated with the immune cell infiltration of endothelial cell, hematopoietic stem cell, macrophage and T cell CD4 + memory resting.

Gene Correlation Analysis and Drug Sensitivity Analysis

The results of gene difference analysis showed that a total of 6 differential genes were screened. Through the single-gene survival analysis of the six differential genes on the gene expression profile interaction analysis website (http://gepia.cancer-pku.cn/), it was found that the high-risk group in the model was positively correlated with various immune checkpoint molecules, such as BIRC5 (P < 0.001), CKLF (P < 0.05), HDAC1 (P < 0.01), NRAS (P < 0.01), PSDM14 (P < 0.01) and SPP1 (P < 0.05) (Figure 8A).

Figure 8 Gene correlation analysis and drug sensitivity analysis. (A), The high-risk group in the model was positively correlated with BIRC5, CKLF, HDAC1, NRAS, PSDM14 and SPP1. (B), The high-risk group was associated with low IC50 of gemcitabine, the low-risk group was associated with low IC50 of erlotinib and axitinib.

Drug sensitivity analysis showed that the high-risk group was associated with low IC50 of chemotherapeutic drug gemcitabine (P < 0.001), while the low-risk group was associated with low IC50 of molecularly targeted drugs epidermal growth factor receptor (EGFR) erlotinib and vascular endothelial growth factor receptor (VEGFR) antagonist axitinib (P < 0.001) (Figure 8B). It showed that the model could be used as a potential predictor of drug sensitivity for the treatment of HCC.

Discussion

As the sixth most common cancer in the world, hepatocellular carcinoma (HCC) is one of the leading causes of cancer death in the world,¹⁸ with more than 700 thousand cases confirmed every year. HCC is a primary inflammation-driven tumor.¹⁹ A large number of studies results have confirmed its immune correlation and further clarified the characteristics and key pathogenic role of immune cells in tumor development.²⁰ At the same time, recent studies have shown that lncRNA plays an important role in regulating pigment, gene expression, growth, differentiation and development.²¹ More importantly, it plays an indispensable role in regulating T cells, B cells, dendritic cells, bone marrow-derived stem cells, TME, epithelial–mesenchymal transformation and microbial population balance.²² Therefore, to explore the role of irlncRNA in the occurrence and development of HCC and predict the prognosis of patients, providing promising biomarkers and immune related treatment targets are expected to become an ideal treatment target for precision medicine.

At present, the research on establishing lncRNA model to predict and evaluate the prognosis of patients with malignant tumors has attracted more and more attention.²³ However, most prognostic models still focus on the combination of several single lncRNAs.^24,25 Usually, there are many inconveniences in application due to sequencing differences and the need for standardized processing of expression profile data. Therefore, the focus of the research is to combine the two lncRNAs and establish a prognosis prediction model for HCC. They are not affected by the expression level. At the same time, because only the relative size of the two lncRNAs is involved and the compatibility is wider, PCR and immunohistochemical data can be applied.

This study was based on the TCGA database (https://portal.gdc.cancer.gov/), and transcriptome data were obtained, and mRNA and lncRNA expression profile data were isolated. Firstly, the CIBERSORT method was used to analyze the immune cell infiltration of liver cancer and normal tissues, and it was confirmed that liver cancer was closely related to immune cell infiltration. Then, the differences between HCC-related irlncRNAs were analyzed, and DEirlncRNA pairs were extracted by iterative cycle and 0 or 1 matrix screening. Then, we screened these DEirlncRNA pairs using univariate and multivariate Cox proportional hazards regression analysis, LASSO regression analysis and cross iterative analysis, and then constructed the prognosis prediction model for HCC. Finally, we verified and evaluated the model from the aspects of survival rate, clinical correlation, invasive immune cells, immune gene correlation and drug sensitivity. In addition, the model was further verified from the perspective of ROC analysis and survival analysis of the experimental group. The results showed that the prognostic model had an excellent predictive effect.

Recent evidence has suggested that lncRNAs contribute to the malignant phenotypes of cancer not only through genomic or transcriptomic alterations but also by altering the immune microenvironment because lncRNAs directly expressed genes related to the activation of immune cells.^26,27 Direct the expression of genes associated with the activation of immune cells.^28,29 This study constructed an IRLP HCC prognosis prediction model including 7 DEirlncRNA pairs, in which ZFPM2-AS1 was a newly discovered lncRNA involved in tumor progression. Zhang et al found that STAT1 up-regulated ZFPM2-AS1 through the mir-653/GOLM1 axis promotes the growth and metastasis of liver cancer cells.³⁰ STAT1 mediated up-regulation of ZFPM2-AS1 promotes the growth and metastasis of HCC cells through mir-653/GOLM1 axis. In the study of Luo et al, it was observed that there were significant differences in the expression levels of lncRNA ZFPM2 AS1 and ZFPM2 genes between HCC tissues and normal liver tissues, and these expression levels were also related to the prognosis of HCC.³¹ The risk assessment model based on 11 lncRNAs was constructed by Li et al to predict the prognosis of HCC, of which AC010280.2 also participated in the construction of the model.³² These findings suggested that ZFPM2-AS1 and AC010280.2 might be potential therapeutic targets and prognostic biomarkers of HCC, and provided theoretical support for the prognostic prediction model constructed in this paper. At present, there are limited verifiable research results on the structure and function of the remaining lncRNA in the model and its role and mechanism in the occurrence and development of HCC,³³ and further exploration is needed. As an inflammation-related tumor, HCC has an immunosuppressive tumor TME,³⁴ which can promote immune tolerance through various mechanisms.³⁵ Since traditional treatment methods for HCC, such as surgical resection, radiofrequency ablation and transcatheter arterial chemoembolization, are effective for localized HCC,^36,37 a variety of systemic drugs targeting TME can be used for unresectable HCC.³⁸ Recently, Liu et al proposed the combination of molecular targeted therapy and immune checkpoint therapy for HCC, which is a new standard treatment for advanced HCC.³⁹ Therefore, in this study, we not only analyzed the correlation and difference of immune cells but also analyzed the drug sensitivity of HCC patients in high- and low-risk groups to provide a certain basis for the further diagnosis and treatment of HCC. The results suggested that there was a certain correlation between the IRLP model and tumor immune microenvironment, and it was found that BIRC5 and other genes were highly expressed in high-risk groups. Xu and others claimed that BIRC5 related to epithelial–mesenchymal transition (EMT) played an essential role in the pathogenesis of HCC.^40,41 This was also consistent with the conclusions of this study. The results of tumor drug sensitivity analysis showed that gemcitabine was more sensitive in the high-risk group, and vascular endothelial growth factor receptor (VEGFR) antagonist axitinib and epidermal growth factor receptor (EGFR) erlotinib were more sensitive in the low-risk group. Therefore, in this study, it was concluded that improving the immune microenvironment and immune response can guide the choice of immunotherapy or combined therapy, effectively enhance the efficacy of immunotherapy, and provide future research direction for accurate treatment and individualized medical treatment of HCC.

In conclusion, our study found a new model that did not need a specific lncRNA expression level. It could predict the prognosis of HCC patients and had significant clinical application value in anti-tumor immunotherapy.

Conclusions

The prognosis prediction model established by irlncRNA pairs has a certain guiding significance for the prognosis prediction of HCC. It may provide valuable clinical applications in antitumor immunotherapy.