The development and validation of a novel senescence-related long-chain non-coding RNA (lncRNA) signature that predicts prognosis and the tumor microenvironment of patients with hepatocellular carcinoma

Abstract

Background: The epigenetic regulators of cellular senescence, especially long non-coding RNAs (lncRNAs), remain unclear. The expression levels of lncRNA were previously known to be prognostic indicators for tumors. We hypothesized that lncRNAs regulating cellular senescence could also predict prognosis in patients with hepatocellular carcinoma (HCC) and developed a novel lncRNA predictive signature.

Methods: Using RNA sequencing data from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) database, a co-expression network of senescence-related messenger RNAs (mRNAs) and lncRNAs was constructed. Using univariate Cox regression analysis and a stepwise multiple Cox regression analysis, we constructed a prognostic HCC senescence-related lncRNA signature (HCCSenLncSig). Kaplan-Meier analysis was used to compare the overall survival (OS) of high- and low-risk groups stratified by the HCCSenLncSig. Furthermore, the HCCSenLncSig risk score and other clinical characteristics were included to develop an HCC prognostic nomogram. The accuracy of the model was evaluated by the time dependent receiver operating characteristic (ROC) and calibration curves, respectively.

Results: We obtained a prognostic risk model consisting of 8 senescence-related lncRNAs: AL117336.3, AC103760.1, FOXD2-AS1, AC009283.1, AC026401.3, AC021491.4, AC124067.4, and RHPN1-AS1. The HCCSenLncSig high-risk group was associated with poor OS [hazard ratio (HR) =1.125, 95% confidence interval (CI): 1.082-1.169; P<0.001]. The accuracy of the model was further supported by ROC curves (the area under the curve is 0.783, sensitivity of 0.600, and specificity of 0.896 at the cut-off value of 1.447). The HCCSenLncSig was found to be an independent prognostic factor from other clinical factors in both univariate and multivariate Cox regression analyses. The prognostic nomogram shows HCCSenLncSig has a good prognostic effect for survival risk stratification. Finally, we found that a higher number of immunosuppressed Treg cells infiltrate in high-risk patients (P<0.001 compared to low-risk patients), possibly explaining why these patients have a poor prognosis. On the other hand, the expression of immunotherapy markers, such as CD276, PDCD1, and CTLA4, was also up-regulated in the high-risk patients, indicating potential immunotherapy response in these patients.

Conclusions: The development of HCCSenLncSig allows us to better predict HCC patients' survival outcomes and disease risk, as well as contribute to the development of novel HCC anti-cancer therapeutic strategies.

Keywords: Hepatocellular carcinoma (HCC); cellular senescence; immunotherapy; long-chain non-coding RNA (lncRNA); prognosis.

Figure 1

Flowchart of the present study. TCGA, The Cancer Genome Atlas; RNA-seq, RNA sequencing; HCC, hepatocellular carcinoma; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNAs, long non-coding RNAs; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; ROC, receiver operating characteristic; TMB, tumor mutation burden.

Figure 2

Identification of senescence-related DEGs. (A) Heatmap of 71 senescence-related DEGs in normal and tumor HCC tissues. (B) Volcano plot of 279 senescence-related genes in normal and HCC tumor tissues. Pink dots represent genes that are up-regulated, while blue dots represent genes that are down-regulated. (C) KEGG analysis of senescence-related DEGs. (D) GO analysis of senescence-related DEGs. N, normal tissues; T, tumor tissues; FDR, false discovery rate; FC, fold change; PD-L1, programmed cell death ligand-1; PD-1, programmed cell death-1; BP, biological process; CC, cellular components; MF, molecular function; DEGs, differentially expressed genes; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Figure 3

Identification of senescence-related lncRNAs related to HCC prognosis and lncRNA-mRNA co-expression network construction. (A) Forest plot showing 54 lncRNAs with hazard ratios, 95 percent confidence intervals, and P values for their associated HCC prognosis based on univariate Cox proportional hazards analysis. (B) Heatmap depicting the expression levels of 8 senescence-related lncRNAs identified by multivariate Cox regression analysis as being associated with HCC prognosis. (C) The lncRNA-mRNA co-expression network of the prognostic senescence-related lncRNA signature. Yellow squares represent prognostic lncRNAs, while green ellipses represent senescence-related mRNAs. Levels of expression of the 8 senescence-related lncRNAs were linked to the levels of 120 senescence mRNAs. (D) Sankey diagram showing the associations between prognostic senescence-related lncRNAs, mRNAs, and risk type. N, normal tissues; T, tumor tissues; lncRNAs, long-chain non-coding RNAs; HCC, hepatocellular carcinoma; mRNA, messenger RNA.

Figure 4

Prognostic value of the risk score determined by the HCCSenLncSig predictive model. (A) Kaplan-Meier curves for OS in the high- and low-risk groups stratified by the median of risk scores determined by the HCCSenLncSig. (B) Risk curve based on the risk score for each sample, where the yellow dot indicates a high-risk and blue dot indicates a low-risk. (C) Scatterplot based on the survival status of each sample. Yellow and blue dots indicate survival and death, respectively. (D) Forest plot for univariate Cox regression analysis. (E) Forest plot for multivariate Cox regression analysis. (F) ROC curve of the risk score and other clinicopathological variables. (G) Time-dependent ROC curves for 1-, 3-, and 5-year survival for the predictive signature. T, tumor; M, distant metastasis; N, lymph node metastasis; CI, confidence interval; AUC, area under the curve; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; OS, overall survival; ROC, receiver operating characteristic.

Figure 5

Visualization of the expression levels of the 8 prognostic lncRNAs based on clinicopathological variable stratification, and PCA of the gene sets performed to classify patients into different risk groups. (A) Eight prognostic senescence-related lncRNAs and clinicopathological variables were distributed in a heatmap for high- and low-risk groups. PCA of low- and high-risk groups based on (B) whole-genome mRNA transcripts, (C) senescence-related mRNAs, (D) senescence-related lncRNAs, and (E) risk model including the 8 HCCSenLncSig senescence-related lncRNAs. Patients with high-risk scores are indicated in red, and those with low-risk scores are indicated in green. N, lymph node metastasis; M, distant metastasis; T, tumor; PC1, principal components 1; PC2, principal components 2; PC3, principal components 3; mRNA, messenger RNA; lncRNAs, long-chain non-coding RNAs; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; PCA, principal component analysis.

Figure 6

Kaplan-Meier survival curves for high- and low-risk groups of patients sorted by clinicopathological variables. (A,B) Age; (C,D) sex; (E-H) grade; (I-L) overall stage; (M-P) T stage. T, tumor.

Figure 7

Construction and verification of the nomogram. (A) Nomogram combining clinicopathological variables and risk score predicts HCC patient 1-, 3-, and 5-year survival probabilities. The calibration curves assess the consistency between the actual and predicted OS at (B) 1 year, (C) 3 years, and (D) 5 years. T, tumor; N, lymph node metastasis; M, distant metastasis; HCC, hepatocellular carcinoma; OS, overall survival.

Figure 8

HCCSenLncSig internal validation for overall survival using the entire TCGA dataset. (A) Kaplan-Meier survival curves for the first internal cohort. (B) Kaplan-Meier survival curves for the second internal cohort. (C) In the first internal cohort, the ROC curve and AUCs at 1-, 3-, and 5-year survival were calculated. (D) In the second internal cohort, the ROC curve and AUCs at 1-, 3-, and 5-year survival were calculated. AUC, area under the curve; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

Figure 9

HCCSenLncSig-based GSEA of different risk groups. (A) According to the GSEA results, KEGG genes were differentially enriched for senescence-related lncRNA expression. Five KEGG items, namely the cell cycle, pathogenic Escherichia coli infection, homologous recombination, Fc gamma R-mediated phagocytosis, and oocyte meiosis, were significantly differentially enriched in the high expression phenotype. Drug metabolism-cytochrome p450, tryptophan metabolism, fatty acid metabolism, tyrosine metabolism, and peroxisome were enriched in the low-risk group based on the NES, NOM P value, and FDR value. (B) Differential enrichment of genes in GO with senescence-related lncRNAs (5 GO items, namely positive regulation of spindle midzone, chloride channel complex, protein depolymerization, regulation of microtubule cytoskeleton organization, and microtubule) revealed significant differential enrichment in the high expression phenotype. The other 5 GO terms, namely acylglycerol homeostasis, monocarboxylic acid catabolic process, organic acid catabolic process, lipid oxidation, and fatty acid catabolic process, were found to be significantly enriched in the low expression phenotype based on the NES, NOM P value, and FDR value. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; BP, biological process; CC, cellular components; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; GSEA, gene set enrichment analysis; lncRNA, long non-coding RNA; Fc, fold change; NES, normalized enrichment score; NOM, nominal; FDR, false discovery rate.

Figure 10

The relationship between HCCSenLncSig risk scores, epigenetic mutation, and TMB. Waterfall plot showing the genetic mutations between high-risk (A) and low-risk (B) HCC patients. (C) Difference in TMB between patients from the low- and high-risk score subgroups. (D) Kaplan-Meier curves for the high and low TMB groups. (E) Kaplan-Meier curves for patients stratified by both TMB and risk scores. The P value represents the ANOVA test between the subgroups. TMB, tumor mutation burden; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; HCC, hepatocellular carcinoma; ANOVA, analysis of variance.

Figure 11

Immune cell infiltration and immune-related functions in different risk groups. (A) The ssGSEA algorithm was used to compute the levels of infiltration of 16 immune cells in the high- and low-risk groups. (B) The relationship between risk score and 13 immune-related functions. (C) The Wilcoxon rank-sum test was used to explore whether there were differences between the 22 kinds of immune cells in different groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, non-significant. aDCs, activated dendritic cells; iDCs, immature dendritic cells; NK, natural killer; pDCs, plasmacytoid dendritic cells; Tfh, T follicular helper; Th1, T helper type 1; Th2, T helper type 2; TIL, tumor-infiltrating lymphocyte; Treg, T regulatory cell; APC, antigen-presenting cell; CCR, C-C chemokine receptor; HLA, human leukocyte antigen; MHC, major histocompatibility complex; IFN, interferon; ssGSEA, single-sample gene set enrichment analysis.

Figure 12

Comparison of immune checkpoints, TIDE scores, and sensitivity of chemotherapy and targeted therapy drugs in high- and low-risk groups. (A) The expression of 27 immune checkpoint genes differed between the high- and low-risk groups. Red boxes represent high-risk patients, while blue boxes represent low-risk patients. (B) The online software “TIDE” predicted the TIDE score of the outcome of subgroups of HCC patients treated with anti-PD-1 or anti-CTLA4 therapies. The IC50 values for (C) 5-fluorouracil, (D) XL-184 (cabozantinib), (E) sunitinib, (F) gemcitabine, (G) paclitaxel, (H) imatinib, (I) bortezomib, and (J) erlotinib in the high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, non-significant; P<0.05 indicates statistical significance. TIDE, tumor immune dysfunction and exclusion module; HCC, hepatocellular carcinoma; PD-1, programmed cell death-1; IC50, half-maximal inhibitory concentration.