Abstract
Background. Despite improvements in cancer screening and diagnosis, hepatocellular carcinoma (HCC) is still diagnosed at an advanced stage, and the prognosis is worse than that of early HCC patients. Therefore, better molecular markers and therapeutic targets in HCC are required.
Objectives. We investigated the predictive value of pyroptosis-related long noncoding RNAs (lncRNAs) in HCC and the effects of these lncRNAs on the immune microenvironment of HCC.
Materials and methods. RNA sequencing data of HCC patients were extracted from The Cancer Genome Atlas (TCGA) database to identify differentially expressed pyroptosis-related lncRNAs related to overall survival (OS). A model was established to verify the character of pyroptosis-associated lncRNAs in the tumor microenvironment, and their prognostic value was evaluated.
Results. A total of 721 PR lncRNAs were identified based on the analysis of the TCGA database. Univariate Cox analysis revealed 37 survival-related PRlncRNAs with prognostic values. As a result of least absolute shrinkage and selection operator (LASSO) regression analysis, ‘ELFN-AS1’, AC099850.3, AC073389.3, ‘HPN-AS1’, AC009283.1, and AL139289.1 showed prognostic value. Kaplan–Meier analysis indicated that the OS of the high-risk set was worse than those of the low-risk set in both the training and testing cohorts. Univariate and multivariate analyses revealed that the risk score was a better independent prognostic factor than the stage. The precision of the lncRNA signature was confirmed using receiver operating characteristic curve (ROC) analysis. Immune- and metabolism-related pathways were enriched in both the low- and high-risk groups. Gene set enrichment analysis suggested that the identified lncRNAs regulate HCC tumorigenesis and prognosis by modulating metabolism. Various algorithms were used to confirm the significant differences in immune cells between these 2 groups.
Conclusions. These findings could contribute to the development and validation of favorable biomarkers, improve the prognosis and survival of HCC, and help in developed individualized treatment plans for HCC.
Key words: prognosis, hepatocellular carcinoma, lncRNAs, pyroptosis
Background
Hepatocellular carcinoma (HCC) is the most commonly diagnosed primary liver cancer and the 3rd leading cause of cancer-related mortality worldwide.1 Surgical resection is still the preferred clinical treatment for early-stage HCC.2 Although recent technological advances have improved cancer screening and diagnosis, HCC, a destructive malignant tumor, is still generally diagnosed at an advanced stage with a 5-year overall survival (OS) of less than 18%.3 Immunotherapy with anti-PD-1, anti-PD-L1 and anti-CTLA-4 antibodies is the most popular systematic treatment for advanced HCC.4 KEYNOTE-224, a phase II trial, revealed a satisfactory curative effect of pembrolizumab against HCC after sorafenib failure. The overall response rate reached 17% with 1 complete response (CR), while the disease control rate was 61%.5 Moreover, the anti-CTLA-4 monoclonal antibody tremelimumab demonstrated a partial response and stable disease rates of 17.6% and 58.8%, respectively, in a small phase II pilot trial (NCT01008358) in advanced HCC patients.6
The major node in the gene regulatory networks was believed to be protein-coding genes. However, noncoding RNAs, with a length of >200 nts, have been reported to play an essential role in tumorigenesis and cancer development.7 Long noncoding RNAs (lncRNAs) mediate and regulate signaling transduction pathways in immune and tumor cells.8 Aberrant expression of lncRNAs can indicate the spectrum of disease progression, and some of these lncRNAs have diagnostic and prognostic value.9, 10 Expression levels of lncRNAs highly correlate with the tumor microenvironment and immune response. The sensitivity of low-risk groups to immunotherapy is high owing to the differences in lncRNA expression levels that lead to differential expression patterns of immune checkpoint genes.11 A detailed understanding of the pathogenesis and molecular mechanisms of HCC is essential to identify novel targets and explore immunotherapy efficacy.
Pyroptosis, a newly discovered process of adjusted cell death, is mediated by various inflammasomes. Pyroptosis initiates the inflammatory process and is involved in the pathogenesis of multiple inflammatory and immune-mediated diseases, such as systemic lupus erythematosus, myositis, acute myeloid leukemia, cervical cancer, and HCC.12, 13, 14, 15 Additionally, pyroptosis is mediated by activation of inflammatory caspases (caspase-1, 4, 5, and 11).16, 17 Several mechanisms regulate the pyroptosis of tumor cells,18, 19, 20 but the principal mechanism is related to several signal transduction pathways and pathologic molecular changes. However, currently, studies on the effects of pyroptosis-associated lncRNAs (PRlncRNAs) in the tumor microenvironment, especially the microenvironment of HCC cells, are limited. In this study, the functions of lncRNAs in pyroptosis were investigated, and their potential mechanisms of action in HCC were elucidated.
Objectives
Hepatocellular carcinoma has a high malignancy rate and poor prognosis, and the knowledge related to etiopathological mechanisms is limited. Our objective was to build a lncRNA model to verify the role of pyroptosis-related RNA in HCC, aiming to clarify immunotherapy responses, guide the treatment of HCC patients and predict the curative effect.
Materials and methods
Identification of PRlncRNAs using RNA sequencing data
The clinical and RNA sequencing (RNA-seq; fragments per kilobase million) data (comprising 374 cancerous and 50 noncancerous samples) of HCC cases were obtained from The Cancer Genome Atlas (TCGA) website (https://tcga-data.nci.nih.gov/tcga/). To identify PRlncRNAs, the correlation between lncRNAs and pyroptosis-related genes (PRGs) was examined employing Pearson’s correlation analysis. Pyroptosis-associated lncRNAs were identified based on the following criteria: p < 0.001 and correlation coefficient |R2| > 0.4.
Unsupervised clustering of PRlncRNAs
To identify different PRlncRNAs, patients were classified for further analysis, including unsupervised cluster analysis. The consistency of the clustering algorithm determines the number and stability of clustering. To this end, we used the ConsensusClusterPlus R package. To ensure the stability of the classification, 1,000 repeats were performed.
Construction and validation of the prognostic PRlncRNA signature
Data of the HCC cases from the TCGA database were categorized into training and testing cohorts using the caret R package (6.0–93). PRlncRNA prognostic significance was evaluated using univariate Cox regression analyses. Significant lncRNAs (as per the least absolute shrinkage and selection operator (LASSO) regression analyses) were selected as candidate characteristic lncRNAs. Furthermore, hub lncRNAs were identified using a multivariate-penalized Cox regression analysis to generate a lncRNA risk signature. The risk score of the prognostic signature was calculated as follows (Equation 1):
, (1)
where ai and xi are the expression and regression coefficient values of each PRlncRNA, respectively.
The risk score was calculated for each patient. Based on the median risk score, patients with HCC were divided into low-risk and high-risk groups. The accuracy of the prognostic signatures was predicted using receiving operating characteristic (ROC) curves. Thereafter, the effect of the PRlncRNA signature on the prognosis of patients was evaluated using Kaplan–Meier analysis. The validity of the model was subsequently verified using a calibration curve and consistency index (C-index). A hybrid nomogram comprised of the PRlncRNA signature and clinicopathological attributes was constructed using the rms R package (6.5–0) to predict the OS of HCC patients.
Gene set enrichment analysis
Gene Set Enrichment Analysis (GSEA; https://www.gsea-msigdb.org/gsea/index.jsp) aims to explore the possible biological functions of lncRNA. The contribution of target genes to a phenotype was evaluated using GSEA by investigating their distribution tendency between the 2 groups according to phenotype relevance. To this end, we utilized the ggplot2 R package (3.3.6) and pRRophetic R package; statistical significances were defined as a nominal p < 0.05 and a false discovery rate <0.25.
Immune response in the 2 risk groups
The heatmap of immune responses was generated using the XCELL, TIMER, QUANTISEG, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT algorithms. The correlation between the risk score values and tumor-infiltrating immune cells (TIICs) was determined using Spearman’s correlation analysis. A heatmap was drawn to depict the distribution differences in immune cells between the low- and high-risk groups and the risk score values.
Statistical analyses
Data analyses were conducted using R v. 4.2.0 (R Foundation for Statistical Computing, Vienna, Austria). The differences between normal and HCC tissues were analyzed using a Student’s t-test or Wilcoxon’s rank-sum test. The LASSO Cox and Cox regression analyses were used to construct a PRlncRNA prognostic model. We tested the proportional hazards assumption for the Cox model and generated plots of the proportional Schoenfeld residuals relative to the transformed time of the residuals for each covariate. Kaplan–Meier and ROC curve analyses aided in verifying the prognostic and predictive power of the risk model. The test comparing Kaplan–Meier curves in this study was the log-rank test. Spearman’s correlation analysis was used to further analyze the association between the TIIC ratio and risk scores. Using the “pRRophetic” package and the expression matrix of liver cancer patients, the minimum inhibitory concentration (IC50) of drugs in high- and low-risk groups of liver cancer patients was predicted, and the drugs with statistically different IC50 values were obtained, which may become candidates for the treatment of liver cancer. The Wilcoxon signed-rank test was applied to compare IC50 between the high- and low-risk patients. All statistical analyses involved a 2-tailed test. A p < 0.05 determined a statistically significant difference.
Results
Unsupervised cluster analysis
The workflow of the study is illustrated in Figure 1. The expression data of 52 PRGs from 50 noncancerous and 374 HCC tissues comprising the HCC cohort were retrieved from the TCGA. In total, 721 PRlncRNAs were identified based on the analysis of these 52 PRGs. Supplementary Table 1 shows the corresponding clinicopathological characteristics of 377 HCC patients. These patients were classified using the ConsensusClusterPlus R package based on the expression of 721 PRlncRNAs. The following 3 clusters were identified using the unsupervised clustering method: clusters 1, 2 and 3 with 145, 105 and 120 cases, respectively. Kaplan–Meier analysis revealed that the cases in cluster 2 were associated with adverse survival outcomes (Supplementary Fig. 1A). The heatmap of the top 231 lncRNAs relevant to favorable and poor prognosis and clinicopathological features is shown in Supplementary Fig. 1B, and the expression of lncRNA in cluster 2 was heterogeneous.
Construction and validation of the prognostic PRlncRNA signature
PRlncRNAs associated with the prognosis of HCC were screened using univariate Cox analysis, which revealed 37 survival-related PRlncRNAs with prognostic values (Figure 2A). The number of lncRNAs was narrowed down to 13 using the LASSO Cox regression analysis (Figure 2B,C). Multivariate Cox analysis identified the following 6 lncRNAs and their correlation coefficients, which were used for subsequent model construction: ‘ELFN-AS1’, AC099850.3, AC073389.3, ‘HPN-AS1’, AC009283.1, and AL139289.1 (Figure 2D). The results of our proportional hazards hypothesis test using the Cox proportional hazards model and the Schoenfeld residuals plot both indicated that the proportional hazards assumption was valid. Additionally, we conducted a linearity test for continuous variables within the Cox proportional hazards model using the easystats package, and the results confirm that the linearity assumption was valid. Consequently, we employed LASSO Cox and Cox regression analysis to construct the PRlncRNA prognostic model.
Prognostic value of PRlncRNAs in HCC
The study cohort was separated into training (256 cases) and testing (109 cases) cohorts in a ratio of 7:3. The risk score was calculated as follows:
Risk score = [‘ELFN-AS1’× 0.05564] + [AC099850.3 × 0.1184] + [AC073389.3 × 1.041] + [‘HPN-AS1’ × (–0.4541)] + [AC009283.1 × (–0.1988)] + [AL139289.1 × 0.6448]
Based on the risk score model, the prognostic risk of HCC patients could be measured, and the median risk score obtained through statistical analysis could be used as the cutoff value between high- and low-risk groups to guide subsequent analyses. Kaplan–Meier survival curve analysis of the 2 cohorts revealed that the PRlncRNA signature was associated with the prognosis of HCC patients, which could also be effectively predicted. In the training and testing cohorts, the prognosis of the high-risk group was significantly poorer compared to the low-risk group (Figure 3A,B). Additionally, the potential expression patterns of survival-related lncRNAs served as early predictors of HCC and were determined using ROC curves. The area under the ROC curve (AUC) at years 1, 2 and 3 was 0.807 (95% confidence interval (95% CI): 0.711–0.904), 0.799 (95% CI: 0.710–0.888) and 0.789 (95% CI: 0.693–0.884) in the training cohort, respectively, which had greater accuracy than that in the testing cohort (AUC at years 1, 2 and 3 was 0.708 (95% CI: 0.535–0.882), 0.724 (95% CI: 0.559–0.889) and 0.719 (95% CI: 0.562–0.876), respectively; Figure 3C,D). The C-index of the risk score model was 0.724 (95% CI: 0.661–0.779).
Supplementary Fig. 2A,C show the survival outcome and risk status of patients in the training cohort. The probability of mortality was higher among patients in the high-risk group than among patients in the low-risk group. Comparable results were achieved in the testing group (Supplementary Fig. 2B,D). The lncRNA expression levels in each group are shown in Supplementary Fig. 2E,F.
The correlation between clinicopathological factors, including clinical stage, risk score and prognosis of HCC patients, was examined using Cox regression analyses. The result revealed that the clinical stage, T stage, M stage, and risk score were independent prognostic factors during univariate Cox regression analysis (Figure 4A), whereas only T stage and risk score remained significant associations during multivariate Cox regression analysis (Figure 4B). The sensitivity and specificity of the lncRNA signature were evaluated using ROC curve analysis. The AUC values of the PRlncRNA signature in the training and testing cohorts were 0.789 (95% CI: 0.693–0.884) and 0.719 (95% CI: 0.562–0.876), respectively, which were more than those of other clinical indices, suggesting the robustness and validity of the lncRNA signature in determining HCC prognosis (Figure 4C,D).
A nomogram was constructed to determine the prognosis of HCC patients based on clinical characteristics and risk scores (Figure 5A). Calibration curves and a C-index were applied to evaluate the predictive capacity of the signature (Figure 5B). The heatmap of clinical features and risk subgroups revealed that the signature was markedly associated with the T stage and risk score (Figure 5C). Thus, the nomogram constructed based on clinical features and risk scores accurately indicated the prognosis of the patients.
GSEA
To identify potential enrichment pathways for PRlncRNA in the risk groups, we used GSEA for biological function analysis. The results showed that phenotypic enrichment of PRlncRNAs in both the high- and low-risk groups was mainly reflected in tumor and metabolism-related pathways, such as apoptosis and the mTOR/NOTCH signaling pathway. By contrast, phenotypic enrichment, although low, was mainly manifested in metabolism-related pathways, such as linoleic acid, drug metabolism and retinol metabolism (Supplementary Fig. 3A). This finding suggests that the prognostic features are related to tumor immunity and metabolism, which helps analyze the principles and mechanisms of the lncRNA signature affecting the metabolism of HCC cells. In addition, testing risk scores can improve the ability to predict chemotherapy efficacy. Among them, the score of the low-risk group was related to the half-maximal IC50 value of aminoimidazole carboxamide ribonucleotide and BI (p < 0.001). In addition, the high-risk group scores were correlated with the IC50 values of bortezomib, LFM-A13 and AKT inhibitors (p < 0.001). These results established that the lncRNA signal can predict the sensitivity to drugs in HCC (Supplementary Fig. 3B–I).
The differential immune cell composition and risk score values between the high- and low-risk groups were examined. Spearman’s correlation analysis was performed using different algorithms. As shown in Figure 6A, a lollipop-shape curve was obtained. Thus, we concluded that most immune cells were significantly and positively correlated with risk scores. Unanimously, compared with that in the low-risk group, the lncRNA signature in the high-risk group was mostly enhanced in immune-associated pathways. Thus, the prognostic PRlncRNA signature plays a critical role in predicting the tumor immune microenvironment, tumorigenesis and progression of HCC.
Additionally, we used distinctive algorithms in constructing the heatmap of immune responses, which is shown in Figure 6B. The proportion of most immune cells (CD4+ T cells, type 2 helper cells, monocytes, resting mast cells, myeloid dendritic cells, and macrophages) in the high-risk group was higher than that in the low-risk group. Compared with those in the low-risk group, the expression levels of immune checkpoint molecules, including CD200, TNFSF9, CD200R1, and TNFRSF4, were upregulated in the high-risk group (Supplementary Fig. 4), indicating that the signature could predict the immunotherapy response in HCC patients. Thus, the PRlncRNA signature may indicate the pyroptosis status of HCC patients and reveal potential biomarkers for clinical therapeutic intervention.
Discussion
The etiology of various tumors, including HCC, is complex. However, hepatitis is a major etiological factor for HCC. Chronic inflammation can lead to dysplasia of liver cells, resulting in HCC. Pyroptosis-mediated inflammatory necrosis may be the key factor regulating tumor cell necrosis.21 To explore the role of pyroptosis in the development of HCC, this study examined the characteristics of PRlncRNAs and PRGs and the clinical prognosis of HCC and elucidated the factors related to immune cell dysfunction in HCC patients.
In this study, 6 lncRNAs (‘ELFN-AS1’, AC099850.3, AC073389.3, ‘HPN-AS1’, AC009283.1, and AL139289.1) were selected and used to generate a risk signature. Additionally, a nomogram comprising the clinical stage and signature was constructed to predict the OS. Collectively, based on these findings, we concluded that this risk signature could be used to determine the prognosis of HCC patients and can be used to develop a staging system. The GSEA was performed to analyze the coordinated expression changes at the pathway levels. The PRlncRNA prognostic signature was correlated with tumor immunity and metabolism. The heatmap of immune cells generated using different algorithms and the expression status of ubiquitous immune checkpoint molecules predicted the immunotherapy response of HCC. Thus, this study demonstrated a strong correlation between gene expression levels and OS. Six differentially expressed lncRNAs were identified as independent prognostic factors for HCC. Among these 6 lncRNAs, some were involved in tumor pathogenesis.
ELFN-AS1, an autophagy-related and immune-related lncRNA, can determine the HCC prognosis.22, 23 AC099850.3 knockdown can promote cell apoptosis by inhibiting the proliferation and metastasis potential of HCC cells. AC099850.3, which mediates its effects through the PI3K/AKT pathway, was positively correlated with key immune checkpoint molecules (PD-1, PD-L1, PD-L2, and CTLA4). Thus, AC099850.3 is a potential immunotherapeutic target for HCC.24 HPN-AS1 can serve as a biomarker for hormone-related cancers, such as prostate, breast and renal clear cell cancers, as its expression level is significantly correlated with OS.25 Additionally, HPN-AS1 functions as a competitive endogenous RNA in HCC patients, and its expression is highly correlated with the tumor mutational burden. The role of AC009283.1 in HCC has not yet been reported. However, AC009283.1 is upregulated in the HER2-enriched cancer subtype. AC009283.1 knockdown modulated the expression of 158 genes, which were enriched in various pathways, including those associated with proliferation, the cell cycle and apoptosis, such as NOTCH3, TNFA, and FOSB, which are reported to regulate proliferation and the cell cycle in cancer.26 AC009283.1 is also considered a ferroptosis-related lncRNA that demonstrates a regulatory role in the tumor microenvironment and immune cell infiltration in colon cancer. Thus, AC009283.1 may aid in determining the prognosis and devising a personalized treatment for colon cancer patients.27 However, the exact mechanisms of AC009283.1 have not been elucidated. Additionally, the oncogenic roles of AC073389.3 and AL139289.1 have not been reported.
From the existing data, it appears that the underlying mechanism by which lncRNAs mediate the pathogenesis of multiple tumors is through the regulation of microRNAs that modulate proteins related to the pyroptosis signaling pathway. Furthermore, the indirect role of lncRNAs in the regulation of pyroptosis can also contribute to the development of tumors. Knockdown of RP1-85F18.6, a newly identified lncRNA with upregulated expression in colorectal cancer tissues, can promote the pyroptosis of colorectal cancer cells by cleaving GSDMD.28 By contrast, XLOC_000647 (an intergenic lncRNA located on chromosome 1) overexpression suppressed the proliferation, invasion and endothelial–mesenchymal transition in pancreatic cancer cells by downregulating NRLP3.29 These findings indicate that lncRNAs have contrasting roles in pyroptosis in different cancers. Nonetheless, within the scope of HCC, Li et al.30 developed a pyroptosis score using principal component analysis to uncover the association between pyroptosis and tumor immunity in individuals with hepatitis B virus-associated HCC. Their research demonstrated that patients with higher pyroptosis scores had a worse prognosis but were more responsive to anti-PD-L1 treatment.
Our study has revealed significant correlations between the signature and numerous immune-associated signaling pathways. Patients with higher risk scores exhibit elevated expression levels of various immune checkpoints, such as CD200, TNFSF9, CD200R1, and TNFRSF4. CD200 and its receptor, CD200R, play a crucial role in the immune system,31 where the CD200R pathway can lead to immune tolerance and inhibit immune function in patients with HCC.32 Huang et al. conducted a study to investigate the impact of the CD200/CD200R pathway on CD4+ T cell subsets in patients with HCC and evaluated the effect of thermal ablation on HCC in rats.33 Furthermore, TNFSF9 represents a potential novel tumor suppressor and therapeutic target for HCC. Activating TNFSF9 signaling could not only serve as a promising adjuvant strategy for immunomodulation therapy but also have a direct inhibitory effect on HCC.34 And based on their findings, TNFRSF4 indeed has a simultaneous impact on both immune cell infiltration and the frequency of gene mutations in HCC.35 Thus, further studies are needed to examine the roles of pyroptosis and lncRNAs in HCC.
Pyroptosis, a newly discovered method of cell death, has dual consequences on tumors. On one hand, key inflammatory bodies in pyroptosis can promote tumor cell death and inhibit the proliferation and metastasis of tumor cells. On the other hand, the accumulation of inflammasomes can create a microenvironment that is conducive to tumor cell growth, proliferation and metastasis. This study aimed to integrate several lncRNA biomarkers to determine their impact on patient outcomes, which could aid in the identification of novel biomarkers and precise therapeutic targets in HCC. Furthermore, the study has the potential to assist in the prognostic prediction, diagnosis and management strategies for HCC patients. However, it is crucial to conduct further independent studies to confirm the results of lncRNAs associated with predictive apoptosis. Although the study explored the potential of developing a prognostic model, further improvement is necessary.
Limitations
This study has some limitations. Currently, lncRNAs have not been utilized as therapeutic targets.36 Our study is established principally on bioinformatics, based on a small number of HCC cases from the TCGA database and limited clinical data; therefore, the results may be biased to a certain extent. The role of pyroptosis RNA in HCC needs to be further validated in clinical practice.
Conclusions
The findings of this study further elucidate the role of PRGs in HCC progression. The establishment of a prognostic model and elucidation of the role of the PRlncRNA signature in immunotherapy response will aid future studies on HCC molecular pathology. Detailed studies will help elucidate the regulatory mechanisms of PRlncRNAs in the pyroptosis pathway of HCC cells. In vivo, modification of lncRNAs using novel therapies is a potential strategy for determining the diagnosis, treatment and prognosis in HCC patients.
Supplementary data
The Supplementary materials are available at https://doi.org/10.5281/zenodo.11638461. The package includes the following files:
Supplementary Table 1. Clinicopathological characteristics of 377 HCC patients.
Supplementary Fig. 1. Identification of the clusters based on expression levels of 721 differentially expressed pyroptosis-related lncRNAs. A. Three clusters were identified using the unsupervised clustering method. B. The heatmap of the top 231 lncRNAs.
Supplementary Fig. 2. Risk survival status chart of TCGA-HCC cohort. The risk survival status plot of patients in the training (A) and testing cohorts (B). Survival outcomes of the patients in the training (C) and testing cohorts (D). Heatmap of the expression of 6 pyroptosis-related lncRNAs in the training (E) and testing (F) cohorts.
Supplementary Fig. 3. GSEA and IC50 values of chemotherapeutic drugs between the high- and low-risk groups based on the pyroptosis-related lncRNA signature in TCGA-HCC cohort. A. GSEA results suggested that the differentially expressed lncRNAs in 2 risk groups were mainly enriched in the tumor-related, immune-related and metabolism-related pathways. B–I Comparative analysis of the IC50 values of chemotherapeutics between the high- and low-risk groups.
Supplementary Fig. 4. Expression levels of immune checkpoints between high- and low-risk groups in patients with HCC.
Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.
Consent for publication
Not applicable.
