The transcriptomic and clinical data (detailed information about the demographic characteristics of the patients in Supplementary Table S3) on CC (normal = 3,tumor = 306) was obtained from the TCGA database (https://tcga-data.nci.nih.gov/), and the immune-related genes were accessed from the Import database (https://www.immport.org/). Count values of raw data were converted to transcripts per kilobase million (TPM) values for subsequent analyses; count values were used only to identify the differential genes. To identify differentially-expressed lncRNAs, we compared different gene expressions between normal and tumor samples with a threshold of |log2 FC (Log2 Fold Change)| > 2 and FDR (false discovery rate) < 0.01. IRLs were obtained based on the relationship between the expression of lncRNAs and immune genes using the Person correlation test (correlation coefficients >0.6). By taking the intersections of DElncRNAs and IRLs, the relevant IRLs were obtained.

Machine learning is widely used in applications such as nearest neighbour search in large-scale data (Yan et al., 2021), dimensionality reduction of features, etc., We filtered the significant prognostic lncRNAs with p < 0.05 through univariate Cox analysis and identified the final lncRNAs using lasso regression analysis. In order to create a prognostic risk model, we used the coefficients obtained from the lasso to calculate the riskscore of each patient with CC. The calculation was as follows: r i s k S c o r e = ∑ i = 1 n C o e f i ∗ x i

The coefficient is Coefi, while xi is the count value of each DElncRNA. Based on the median riskscore (The advantages of the median are that it makes full use of all data information to reflect the centralized trend of a group of data, is not affected by extreme data, and is easy to find. It can clarify the middle level and is less affected by extreme data. Disadvantages: It is easily affected by extreme values), patients were divided into two groups, namely, the high-risk and the low-risk group. Furthermore, based on the survival duration of patients, a Kaplan-Meier analysis was performed to evaluate their prognostic value. Model-based receiver operating characteristic (ROC) curves were plotted for the first, third and fifth year; based on the survival time of patients, Kaplan-Meier analyses were performed, and survival curves were used to display the risk model results.

We collected several common oxidative-stress-related genes from published studies, including SOD1, SOD2 (Zelko et al., 2002), PON1(Teranishi et al., 2012), NOS3(Katkam et al., 2018), UCP2(Hu et al., 2019), GSR (Couto et al., 2016), GPX1 (Teranishi et al., 2012) and GSTM1(Cadoni et al., 2006). 8-hydroxy-2 deoxyguanosine (8-OHDG) is known as a key marker of oxidative stress (Reuter et al., 2010) and we collected 8-OHDG-related genes from GeneCards (https://www.genecards.org/) and obtained gene enrichment pathways using ClueGO. A protein-protein interaction network was then obtained through the String website (https://cn.string-db.org/). Four methods based on the naximal clique centrality (MCC), the density of maximum neighborhood component (DMNC), the maximum neighborhood component (MNC) and degree of cytoHubba were used to screen key genes, the top 10 of which were crossed. Correlations between the risk model and oxidative-stress-related genes were analyzed using the Spearman test.

The correlations between the model and the age, grade, clinical stage as well as TNM stage of patients with CC were assessed using the Chi-square test.

Microsatellite instability (MSI), tumor mutational burden (TMB) and homologous recombination deficiency (HRD) are common molecular characteristics of genomic instability, which are validated biomarkers for targeted therapies (26). We performed a Kaplan-Meier analysis for TMB, MSI and HRD to assess their prognostic values. A Chi-square test was also used to evaluate their associations with the risk model.

We downloaded the data on SNPs and CNVs of patients with CC from TCGA and UCSC databases. The SNPs and CNVs were visualized using circos (http://circos.ca) and R. The focus was on the demonstration of their locations on the chromosomes where the genes were present in the model. Significantly-mutated genes (p < 0.05) and gene mutation interactions between the high- and low-risk group were analyzed using the MAFTOOLS software. In both analyses, only the genes mutating more than 10 times in at least one group were considered, whose expression was probed using GEPIA. A statistical test for significant mutation rates was performed using a one-sided z-test. Copy number alterations among patients with CC were analyzed with GISTIC 2.0 (Mermel et al., 2011). The copy number gistic score, together with the percentage of patients in both risk groups, was also analyzed.

We obtained most of our immune cell data from XCELL (Aran et al., 2017), EPIC(Racle et al., 2017) and CIBERSORT (Newman et al., 2015). Next, the immune cell infiltration was quantified using ssGSEA for subsequent analyses (Barbie et al., 2009; Bindea et al., 2013). Further analyses were conducted on the correlations between immune cell types and immune cell content of the risk groups. Based on a Pearson correlation analysis, we analyzed how immune cells and IRLs interacted.

We used unsupervised subclass mapping methods (https://cloud.genepattern.org/gp) to predict the responses of different risk groups to immunotherapies (Lu et al., 2019).

PIK3CA mutations are more common in CC (Xiang et al., 2015). As a result, mutations in the PIK3CA gene can be used as a target biomarker for patients with CC. We segregated the patients with CC carrying PIK3CA mutations. After downloading data on drug sensitivity AUC value from the Cancer Therapeutics Response Portal (CTRP2.0), Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) repurposing dataset and the Cancer Cell Line Encyclopedia (CCLE) expression profile, we performed a differential drug response analysis on patients in both risk groups. Spearman correlation analysis (r < −0.30 for CTRP; r < −0.35 for PRISM) was used to screen compounds with negative correlation coefficients (Yang et al., 2021).

A total of 22 pairs of cancerous and non-tumorous tissues were collected from patients with CC at the First Affiliated Hospital of Zhengzhou University. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (Ethics Number: 2022-KY-0093-002), and an informed consent was obtained from patients. qRT-PCR was used to determine the expression of lncRNAs through the risk model. Primers for TFAP2A-AS1, AP000911.1, AL133215.2 and LINC02078 were designed using primer 5.0 (Supplementary Table S1). Total RNA was extracted with a trizol (CWBIO, China), and the first strand of cDNA was synthesized using a reverse transcription kit (Takara, Kyoto, Japan). Finally, cDNA was quantified through qRT-PCR using SYBR green master mix (Vazyme, China). GAPDH was used as an internal reference for calibration. The 2^−ΔΔCT method was chosen to calculate the relative expression of lncRNAs. Cellular functional assays were performed for AL133215.2, which was knocked down in the HeLa and SiHa cell line via transfection. Cells transfected with siRNA and controls were stained with an Annexin V-FITC apoptosis detection kit (Beyotime, Shanghai, China). The stained cells were then analyzed through flow cytometry. Relative cell viability was monitored 24, 48, 72, and 96 h after transfection using cell counting kit-8 (CCK-8, Beyotime, Shanghai, China).

We used GSVA to analyze the 50 hallmark pathways described in the molecular signature database (Subramanian et al., 2005; Hänzelmann et al., 2013). Next, we used a limma package to obtain pathways that differed significantly between patients in the high- and low-risk group. A GSEA (Subramanian et al., 2005) was conducted for both risk groups, and we selected significantly-enriched pathways based on p-values and FDR q-values that were below 0.05 and 0.25 respectively. We obtained the previously-reported gene sets related to immunotherapy from Hu et al. (2021). In addition to the immune-related gene sets we collected, gene ontology (GO) pathways associated with oxidative stress were also enriched using GSVA. Finally, we examined the associations between genes in the model and the enrichment scores.

The nomogram-integrated factors including the risk score, T, N, MSI. calibration curves and ROC were used to evaluate the accuracy and predictive ability of the nomogram; decision curve analysis (DCA) was used to evaluate the clinical effectiveness of the nomogram.

R version 4.4.1 was used to perform all statistical tests in this manuscript. The χ2 test was used for appropriate categorical data, and the two-sample Wilcoxon test was used for continuous data. Survival analyses were performed using the R package “survival”. Correlation analysis was performed using the Pearson correlation test. Statistical significance was defined as a p-value of less than 0.05 for all statistical analyses.

We obtained a total of 493 differentially-expressed lncRNAs; among them, 96 were immune-related lncRNAs (Supplementary Table S2). A total of four IRLs were identified and selected for risk modeling through univariate cox analysis (Figure 1A) and lasso analysis (Figures 1B, C); among them, TFAP2A-AS1 and AL133215.2 showed a significantly high expression, while LINC02078 and AP000911.1 had a significantly low expression in cancer tissues (Figure 1D). The risk score of each patient was computed as follows: Riskscore = TFAP 2 A − AS 1 * − 0.207 + AP 000911.1 * − 0.214 + AL 133215.2 * − 0.229 + LINC 02078 * − 0.308

The risk model had a good clinical predictive power, with a ROC value of 0.763, 0.645 and 0.678 for 1-, 3- and 5-year survival respectively (Figure 1E). The C-index and IBS of the risk model were 0.918 and 0.035, respectively, which also show that the risk model had a good predictive performance (details of the calculation process were in Supplementary Material S1). Patients in the low-risk group had a higher overall survival rate (Figure 1F). Furthermore, we determined the distribution of risk scores, the survival statistics of patients in different risk categories and the expression characteristics of the four IRLs (Figures 1G–I). As is shown in the graph, the low-risk patients showed an overexpression of these four protective lncRNAs.

Among the eight common oxidative stress-related genes, the risk model was positively correlated with SOD2 while negatively correlated with UCP2 (Figure 2A). One of the four lncRNAs, AL133215.2, correlated significantly positively with SOD2 while negatively with UCP2; another was TFAP2A-AS1, which correlated significantly positively with UCP2. (Figure 2B). There was a significant positive correlation between SOD1 and GPX1 among the eight genes related to oxidative stress (Figure 2C). We obtained 62 8-OHDG-related genes from GeneCards, and the protein interaction network showed more interactions among TP53, OGG1, SOD2, CAT and other proteins (Figure 2D). The pathway enrichment results showed that 8-OHDG-related genes were significantly enriched in the negative regulation of oxidative stress-induced intrinsic apoptotic signaling pathways, responding to oxidative stress and other pathways (Figure 2E). Five key 8-OHDG-related genes were obtained through MCC, DMNC, MNC and Degree, namely, SOD2, OGG1, TP53, NFE2L2 and CAT (Figure 2F). The risk model and OGG1 were significantly negatively correlated (Figure 2G), and AL133215.2, TFAP2A-AS1 as well as OGG1 were significantly positively correlated (Figure 2H).

Among all clinical characteristics, N and T stage were significantly associated with the risk model (Figure 3). More patients in the low-risk group were in Stage N0, more of whom in the high-risk group were in Stage T3 and T4. Perhaps this is one of the reasons for the shorter overall survival time of patients in the high-risk group.

Of the three markers, MSI and the risk groups were significantly correlated (Figure 4A). The high-risk group had a higher MSI value (Figure 4B). However, there was no significant difference between the two groups for TMB (Figure 4C) or HRD (Figure 4D). In a combined analysis on MSI and riskscore, highMSI-highRisk patients had a significantly shorter survival duration than lowMSI-lowRisk ones (Figure 4E). These results indicated that high-risk patients tended to have a higher MSI value, which contributed to a poorer prognosis.

From the analysis, TTN, PIK3CA and KMT2C had a high mutation frequency in both high- and low-risk group (Figure 5A, B). There were two genes with a high frequency of mutations in high-risk patients, namely, DNAH2 and AHNAK. There were 12 genes with a high mutation frequency in low-risk patients, namely, DNAH2 and AHNAK (Figure 5C). These genes exhibited a significant co-occurrence (Figure 5D). AHNKA, DMD, CACAN1H, KINA1109 and BRCA1 were differentially expressed between the CC and normal group (Figures 5E–I). CNV results from patients with CC showed a greater significant increase in gene copy number on chromosomes 1 and 3 (Figure 5J). The CNV of chromosomes 6, 10, 11, and 17, where the four lncRNAs TFAP2A.AS1, AL133215.2, AP000911.1, and LINC02078 were located respectively, are shown in Figure 5K. As is shown in the figure, there was a higher gene copy number loss in these chromosomes. A significant similarity in the chromosomal aberrations of patients in the high- and low-risk group was also observed (Figure 5L).

Among the immune cell types, assessed by several methods, a significant positive correlation was shown between neutrophils and riskscore (Figure 6A). When ssGSEA was used to quantify immune cells, most of them showed a strong positive relationship with each other (Figure 6B). Neutrophils, NK cells and pDCs were significantly different between the high- and low-risk group (Figure 6C). Among the four lncRNAs in the risk model, the expression of TFAP2A-AS1 and AL133215.2 was negatively correlated with that of the majority of immune cells, while LINC02078 and AP000911.1 showed opposite trends (Figure 6D). We used the submap analysis to compare the gene expression profiles of the defined high- and low-risk group with another dataset containing 47 melanoma patients who showed immunotherapeutic responses (Roh et al., 2017). Anti-PD-1 therapies were more likely to be effective for low-risk patients (Bonferroni corrected p = 0.015) (Figure 6E). We screened several drugs that showed sensitivity among low-risk patients, including bleomycin A2, dasatinib and afatinib in CTRP2.0 (Figure 6F) as well as NVP-AUY922, dinaciclib, pelitinib, obatoclax, echinomycin, dasatinib and dacomitinib in PRISM (Figure 6G). Based on these studies, we could develop individualized treatment plans for patients in different risk groups.

The qRT-PCR assay indicated that TFAP2A-AS1 and AL133215.2 showed a significantly high expression in cancer tissues, while LINC02078 and AP000911.1 had a significantly low expression, which were consistent with our previous predictions (Figures 7A–D). Given the above analyses, we found a close association between AL133215.2 and oxidative-stress-related genes, so AL133215.2 was selected for experimental validation. CCK-8 assay showed that the proliferation of CC cells could be significantly inhibited by knocking down AL133215.2 (Figures 7E, F). We analyzed CC cells through flow cytometry and showed that apoptosis was significantly accelerated in both Siha and Hela cells after AL133215.2 knockdown (Figures 7G–I). These experimental results demonstrated the ability of AL133215.2 to promote the progression of CC.

GSVA was performed on each patient to compare pathways between high- and low-risk group. Among high-risk patients, inflammatory response and oxidative phosphorylation were significantly enriched (Figure 8A). According to GSEA results of the HALLMARK and KEGG gene set, inflammatory response and oxidative phosphorylation were significantly enriched in high-risk individuals (Figures 8B–D). A positive correlation was found between AL133215.2 and most immunotherapeutic predictive pathways, while a negative one was seen in all immune precursor pathways (Figure 8E). The AL133215.2 low-expression group had a higher enrichment score in terms of the oxidative stress pathways (Figure 8F). This was consistent with our previous study, where the high-risk group corresponded to a lower expression of AL133215.2 and a more pronounced oxidative stress.

A nomogram was constructed to predict survivals after 3 and 5 years (Figure 9A). The calibration curve validated that the nomogram had a good accuracy in predicting patient survivals (Figure 9B). ROC and DCA showed that the risk score had a better clinical efficacy compared to several other characteristics (Figure 9C, D).