Enrichment Analysis of Oxidative Stress-Related Genes (OS-DEGs)
We obtained 127 oxidative stress-related differentially expressed genes (OS-DEGs), including 15 down-regulated genes and 112 up-regulated genes (Fig. 1A and Supplementary Table 2). The gene enrichment results of Kyoto Encyclopedia of Genes and Genomes (KEGG)show that OS-DEGs are mainly enriched in endocrine resistance, lipid and atherosclerosis, fluid shear stress and atherosclerosis, chemical carcinogenesis-reactive oxygen species, Kaposi sarcoma-associated herpesvirus infection, IL-17 signaling pathway, Chagas disease, Hepatitis B, Prolactin signaling pathway, Pancreatic cancer (Fig. 1B). In the biological process category, Gene Ontology (GO)analysis showed that DEGs were mainly enriched in the response to oxidative stress, the response to xenobiotic stimulus, the cellular response to chemical stress, the cellular response to xenobiotic stimulus and the cellular response to oxidative stress. In the cellular components category, GO analysis showed that DEGs were mainly enriched in neuronal cell body, distal axon, cation channel complex, membrane raft and membrane microdomain. In the molecular function category, GO analysis showed that DEGs were mainly enriched in monooxygenase activity, heme binding, tetrapyrrole binding, oxidoreductase activity and steroid hydroxylase activity (Fig. 1C).
Construction of the Oxidative Stress-Related lncRNAs Predictive Signature
We found 583 oxidative stress-related lncRNAs by correlation analysis. After Univariate Cox regression analysis, we identified 68 lncRNAs associating with prognosis of HCC patients (Supplementary Table 3). Multivariate Cox regression analysis revealed that 8 oxidative stress-related lncRNAs (OS-lncRNAs) (TMEM22-AS1, AC012360.2, SNHG3, GASAL1, PCAT6, AC009006.1, AL031985.3, AC009403.1) were used to construct predictive signature(Supplementary Table 4). The interaction analysis between lncRNA and mRNA is visualized by Cytoscape as shown in Fig. 2A. The co-expression network contained lncRNA-mRNA were shown in Fig. 2B. LncRNA AC009005.1 is the target lncRNA which we were interested in. AC009005.1 had co-expressive relationship with 12 OS-DEGs (DDIT3, HSPB1, STIP1, NUDT1, CDKN2A, HRAS, BAK1, FANCD2, TRAF2, Fig. 3.
Evaluation of Predictive Signature
According to the calculations based on the risk score formula, all patients were divided into high- or low-risk groups based on the median risk score. The Kaplan-Meier survival curve was used to analyze the OS time of high and low risk groups. The results showed that compared to the high-risk group, the low-risk group had a shorter survival time (Fig. 4A, p < 0.001). It can be seen that as the risk score increases, the mortality rate of patients also increases (Fig. 4B). The risk scores of the high- and low-risk groups are shown in Fig. 4C. The AUCs of the 1-, 3-, and 5-years survival were 0.764, 0.715, 0.759, which meant our predictive signature has good predictive ability (Fig. 4D). Figure 4E showed the expression levels of OS-lncRNAs in different risk groups, we can see that the expression difference of lncRNA AC009005.1 is more significant between the high and low risk groups compared to other OS-lncRNAs.
Correlation between the Predictive Signature and the Prognosis of HCC Patients
To determine whether the prognostic model can serve as an independent prognostic factor for HCC patients, Cox regression analysis was conducted. Univariate Cox regression analysis showed that stage, T stage, and risk score were significantly correlated with OS in HCC patients (Fig. 5A). Multivariate Cox regression analysis showed that risk score were independent predictors of OS in HCC patients (Fig. 5B). The AUCs of the risk score, stage, T stage were 0.822, 0.721, 0.711, which was better than other clinicopathological variables in predictive signature of HCC patients (Fig. 5C). We visualized the clinical pathological variables of the high and low risk groups, and the results are shown in Fig. 6.
Nomogram of the Predictive Signature
To further predict the prognosis of HCC patients, we developed a nomogram which includes risk score and clinicopathological variables, that can predict the prognosis of HCC patients at 1-, 3-, and 5-years (Fig. 7A). According to the 1-, 3-, and 5-years calibration curve, we found that there is good consistency between the predicted survival rates and the actual OS rates (Fig. 7B-D).
Survival analysis of clinicopathological variables in predictive signature
To verify the impact and correlation of each clinicopathological variables on patients’ survival time, we plotted Kaplan-Meier survival curves for patients in high- and low-risk groups under different clinical clinicopathological variables conditions (Fig. 8). For each different clinical case variable, the OS of high-risk group patients was significantly shorter than that of low-risk group patients. Among them, male patients have a poor prognosis, which may be related to their preference for alcohol consumption (Fig. 8C-D). The poor significance of Stage Ⅲ-Ⅳ and T Stage 3–4 may be due to a small number of clinical samples (Fig. 8J and L). The above results indicate that predictive signatures can effectively predict the prognosis of HCC patients.
Internal Validation of the Predictive Signature
To verify that the prognostic signature we constructed is applicable to the entire TCGA-LIHC dataset (n = 377), after deleting the patients with OS time less than 30days, we randomly divided HCC patients (n = 343) into two groups, with 171 or 172 patients in each group. The demographic characteristics of patients in the two cohorts are shown in Supplementary Table 5. As expected, in all internal validation sets, the OS rate in the high-risk group was lower than that in the low-risk group (Fig. 9A-B). The ROC results show that the AUC values of all internal validation sets are above 0.6, demonstrating good predictive ability (Fig. 9C-D)
Gene Enrichment Analysis
Because of the different prognoses of patients in the high- and low-risk groups, we conducted GSEA to study the possible differences between the high- and low-risk groups. We found that in the KEGG pathway, genes in high-risk groups are mainly enriched in cell cycle related signaling pathways such as Base Exception Repair, Cell Cycle, p53 Signaling Pathway, and VEGF Signaling Pathway (Fig. 10A). In the GO pathway, genes in high-risk groups are mainly related to molecules that affect cell cycle, such as Nuclear Ubiquitin Ligase Complex, Positive Regulation of Cell Cycle Phase Transition, DNA Polymerase Complex, Regulation of Cell Cycle Phase Transition, etc. (Fig. 10B).
Analysis of immune cell infiltration, immune cell function, and immune checkpoint
To further guide clinical treatment, we analyzed the infiltration of immune cells, immune cell function, and expression of immune cell checkpoints between high- and low-risk groups. The results showed that there were significant differences in activated dendritic cells (aDCs), B cells, macrophages, and neutrals among the high-and low-risk groups, indicating that patients in the high- and low-risk groups may have different abilities in antigen uptake and presentation (Fig. 11A). In terms of immune cell function, we found that there were significant differences in cytolytic activity, MHC-Ⅰ molecule, and type II IFN response among patients in the high- and low-risk groups (Fig. 11B). Immune checkpoint analysis showed significant differences in the expression levels of all immune checkpoint molecules in the high and low risk groups, indicating that using different immune checkpoint therapy drugs for different risk groups may have better effects (Fig. 11C).
Sensitivity prediction of clinical therapeutic drugs
In order to provide further guidance for clinical treatment, and for better targeted treatment, achieve precise medical purposes, we selected clinical common anti-tumor drugs (Supplementary Table 6), and analyzed the sensitivity of high- and low-risk HCC patients to verify the possibility of constructed prognostic markers for clinical treatment. The prediction results show that in tyrosine kinase inhibitor drugs such as Nilotinib, Gefitinib, Erlotinib, Axitinib, etc., high-risk patients will be more sensitive. For patients in the low-risk group, drugs such as Tipifarnib, a farnesyl protein transferase inhibitor, and Mitomycin C, an antibiotic chemotherapy drug, are more effective (Fig. 12). The prediction results of the differences in drug sensitivity between the high-risk and low-risk groups may help to explore individualized clinical treatment plans for HCC.
The Expression of AC009005.1 and Relative mRNA in HCC Cell Lines and Tumor Tissues
The expression of AC009005.1 is one of the most upregulated OS-lncRNAs in prognosis. Therefore, in order to further evaluate the role of AC009005.1 in HCC, we used the TCGA database and RT-qPCR to evaluate the expression level of AC009005.1. Compared with adjacent normal liver tissues, AC009005.1 mRNA expression is higher in HCC tissues (Fig. 13A). Similar results were also found in different HCC cell lines (Fig. 13B). We also detected the expression of 12 related mRNA in the AC009005.1 co expression network, as shown in Fig. 13C. Interestingly, we found that most genes are associated with oxidative stress and cell cycle, so we believe that AC009005.1 may affect cell cycle by affecting intracellular ROS levels.
Downregulation of AC009005.1 increases ROS levels and G2/M phase of cells
As we mentioned earlier, excessive accumulation of ROS can have a toxic effect on cells, inducing them to die. So, we detected ROS levels in HCC cell lines. As the result showed, the ROS levels of HepG2 cells transfected with siRNA-3 were significantly higher than those of control transfected cells (Fig. 14A). We also tested the cell cycle and the results showed that the S phase of the cell cycle decreased while the G2/M phase increased, indicating that AC009005.1 may reduce S phase division and cause G2/M phase arrest, affect the cell cycle, inhibit cancer cell proliferation (Fig. 14B).
Downregulation of AC009005.1 Inhibits HCC Cell Migration, Invasion, and Proliferation
To further confirm the role of AC009005.1 in migration and invasion, transwell assays and wound healing assays were performed. Our results showed that the invasion and migration rates of HepG2 cells transfected with siRNA were significantly lower than that of the control-transfected cells (Fig. 15A). Wound healing assays revealed that downregulate AC009005.1 significantly repressed wound healing in HepG2 cells (Fig. 15B). Cell proliferation and colony formation assay also confirmed the above results (Fig. 15C-D).