Differential expression analysis of irlncRNAs (DEirlncRNAs)
We downloaded transcriptome data of HCC from the TCGA database, including 50 non-tumor tissues and 374 tumor tissues. The gene transfer format (GTF) files from Ensemble were used to annotate the accessed data, and co-expression of known ir-genes and lncRNAs was analyzed. Totally, 740 irlncRNAs were identified, of which 490 were classified as DEirlncRNAs. Among the DEirlncRNAs, 16 were down-regulated, and 474 were up-regulated (Figure 1A, 1B, Table S1, S2).
DEirlncRNA pairs screening and risk model construction
A total of 10344 valid DEirlncRNA pairs from 490 DEirlncRNAs were identified by iteration loop and 0 or 1 matrix. 1009 DEirlncRNA pairs were screened by a single factor test and modified LASSO analysis, of which 30 pairs were involved in the Cox model by the stepwise method. Results were shown in Figure 1C, 1D, and 1E. After that, the areas under the curve (AUC) for each receiver operating characteristic (ROC) curve of 1009 were calculated, and the curve was plotted. At the same time, we found that the maximum AUC value was obtained when the highest point equal to 0.941, then the the optimal DEirlncRNA pair was determined (Figure 2A). Our study also used the Akaike information criterion (AIC) values to determine the maximum inflection point as the cut-off point of the 5-year ROC curve (Figure 2B). ROC curves at 1, 3, and 5 years were drawing to verify the optimality, which indicated all AUC values exceeded 0.91 (Figure 2C). Moreover, AUC values between the 5-year ROC curve and some common clinical parameters were also compared. (Figure 2D). Furthermore, 343 HCC patients’ data were collected from the TCGA database, and the risk scores of all patients were calculated. Then the cut-off point was utilized to re-differentiate the high-risk group and low-risk group for verification.
Application of risk assessment model in the clinical evaluation
Fifty-nine cases and 284 cases were divided into high-risk and low-risk groups according to the cut-off point. Figure 3A and Figure 3B display the riskScores and survival rate of each case. These data indicate that patients’ clinical outcome in the high-risk group was inferior to that in the low-risk group. The survival of the high-risk group was poor than that of the low-risk group by Kaplan Meier analysis (p < 0.001) (Figure 3C). Next, a Chi-square test was conducted to explore the relationship between the risk of HCC and clinicopathological features. The stripping diagram (Figure 4A) and the scatter plots examined through Wilcoxon signed-rank test displayed T classification, tumor stage, tumor grade, Child-Pugh grade, ECOG, vascular invasion, and survival status (Figure 4B-4H) were significantly associated with the risk. Subsequently, univariate Cox regression analysis indicated that there were statistical differences in tumor stage (p < 0.001, HR = 1.627, 95% CI [1.218 – 2.173]), vascular invasion (p = 0.007, HR = 1.737, 95% CI [1.161 – 2.599]), ECOG (p < 0.001, HR = 1.945, 95% CI [1.346 – 2.811]), and riskScore (p < 0.001, HR = 1.028, 95% CI [1.017 – 1.039]), while tumor stage (p = 0.018, HR = 1.488, 95% CI [1.072 – 2.067]), ECOG (p < 0.001, HR = 1.905, 95% CI [1.299 – 2.792]), and riskScore (p < 0.001, HR = 1.024, 95% CI [1.015–1.034]) can be regarded as independent prognostic predictors by multivariate Cox regression analysis (Figure 4I).
Relationship between tumor-infiltrating immune cells, immune molecules and the risk model
Since lncRNAs were initially associated with ir-gene, we explored whether this model correlates with the tumor immune microenvironment. Our study found that there is a positive correlation between the high-risk group and tumor-infiltrating immune cells, such as B cells, Neutrophils, and macrophages, while negatively correlated with CD8 + T cells, CD4 + T cells, and monocytes. Spearman correlation analysis was carried out in detail, and the results were shown as a lollipop diagram (Figure 5A). Since ICIs currently play a significant role in the treatment of HCC, we explored the correlation between the risk model and ICI-related biomarkers. The results showed that high risk scores were associated with the expression of CD276 (p < 0.001), GSDME (p < 0.001), HAVCR2 (p < 0.01), and TNFRSF18 (p < 0.05) (Figure 5B-5E). There was no statistical difference between the high-risk scores and ir-genes, such as CTLA4, PDCD1, and LAG3 (All p > 0.05, Figure 5F-5H).
Correlation analysis of risk model and chemotherapy drugs
In addition to checkpoint blocking therapy, we also tried to explore the association between the risk and efficacy of commonly used chemotherapy drugs for HCC. We found that the high-risk score was along with a higher half inhibitory concentration (IC50) of chemotherapy drugs for erlotinib (p < 0.001), methotrexate (p = 0.0024), and rapamycin (p = 0.0043), while a lower IC50 for bleomycin (p < 0.001), doxorubicin (p = 0.0023), gemcitabine (p < 0.001), mitomycin (p < 0.001), and paclitaxel (p = 0.011), suggesting that this model can be served as a potential predictor for the sensitivity of chemotherapies (Figure 6A-6J).