Identification and validation of oxeiptosis-associated lncRNAs and prognosis-related signature genes to predict the immune status in uterine corpus endometrial carcinoma

Screening of prognostic lncRNA candidates

Through the correlation analysis with |R²| > 0.2 and p < 0.05, 723 lncRNAs significantly associated with oxeiptosis genes were identified (Supplementary Table 1) and used for subsequent explorations. Meanwhile, we identified 158 differentially expressed lncRNAs (Supplementary Table 2) and 22 prognosis-associated lncRNAs (Supplementary Table 3) using differential expression and univariate Cox regression analyses. Finally, 8 overlapping lncRNAs were selected as candidate lncRNAs for further prognostic analysis (Figure 1A).

Figure 1. Identification of prognostic oxeiptosis-associated lncRNAs. (A) Venn diagram of candidate oxeiptosis-associated lncRNAs determined by differential expression and univariate Cox analyses. (B) Correlation network of prognostic lncRNAs and their associated mRNAs. (C) Correlation network of hub lncRNAs.

Risk signature construction

Using the above-identified candidate lncRNAs, we conducted Lasso penalized Cox regression analysis and constructed a risk signature with 5 hub lncRNAs: HOXB-AS3, AC009097.2, AL359220.1, AC100861.1, AC245884.9 (Supplementary Table 4). The connections of these identified hub lncRNAs to the genes related to oxeiptosis are presented in Figure 1B, 1C.

Prognostic value analysis of oxeiptosis-associated lncRNAs in UCEC

When evaluating the gene expression levels of hub lncRNAs in UCEC, the results show that AC009097.2 (Figure 2A), AC100861.1 (Figure 2B), and HOXB-AS3 (Figure 2E) were significantly higher expressed in UCEC tissues compared to normal samples (p < 0.05). In contrast, the others, including AC245884.9 (Figure 2C) and AL359220.1 (Figure 2D), were expressed at significantly lower levels in UCEC tissues (p < 0.05). The KM survival analysis was further applied to estimate the correlation between lncRNA expression and UCEC prognosis. The results indicated that the higher AC009097.2 (Figure 2F), AL359220.1 (Figure 2I), and HOXB-AS3 (Figure 2J) expression subgroups had a significantly higher OS of UCEC patients (p < 0.05). On the other hand, the higher expression of AC100861.1 (Figure 2G) and AC245884.9 (Figure 2H) were dramatically connected to the poor prognosis in UCEC patients (p < 0.05). As a single diagnostic biomarker, the ROC analysis revealed that the area under the ROC curve (AUC) of HOXB-AS3 was 0.751, AC009097.2 was 0.743, AL359220.1 was 0.824, AC100861.1 was 0.807, and AC245884.9 was 652, indicating that all identified lncRNAs had a high predictive accuracy in UCEC patients (Figure 2K). Meanwhile, when all hub lncRNAs were combined into a prediction model, the ROC analysis showed that the predictive accuracy of UCEC increased to 0.949 (Figure 2L).

Figure 2. Clinical value of hub lncRNAs in UCEC. Gene expression levels of AC009097.2 (A), AC100861.1 (B), AC245884.9 (C), AL359220.1 (D), HOXB-AS3 (E) in risk subgroups. Survival curve of AC009097.2 (F), AC100861.1 (G), AC245884.9 (H), AL359220.1 (I), HOXB-AS3 (J) in UCEC. ROC curves of single diagnostic biomarkers (K) and prediction model (L).

Risk signature value in clinics

Based on the median value of calculated risk scores, we categorized UCEC patients into two subgroups with low- and high-risk scores (Figure 3A, 3B). While evaluating the expression level of lncRNAs in two risk subgroups, Figure 3C disclosed that AC100861.1 and AC245884.9 were expressed at significantly higher levels in high-risk subgroup (p < 0.05). In contrast, the others, including AC009097.2, AL359220.1, and HOXB-AS3, were all significantly lower expressed in high-risk subgroup compared to low-risk subgroup (p < 0.05). The UCEC patients with high-risk scores showed a significantly lower OS than those with low-risk scores (p < 0.05; Figure 3D). Univariate Cox regression analysis also confirmed that the risk signature was significantly connected with OS of UCEC patients (Figure 3E). Additionally, Figure 3F also demonstrated that this risk signature could be used as an independent factor for predicting UCEC patients. A high predictive accuracy of this risk signature at 1 (AUC = 0.849), 3 (AUC = 0.730), and 5 (AUC = 0.760) years were found using ROC curve analysis (Figure 3G). Compared with other traditional clinicopathological features (including age, TNM stage, and cancer grade), a significantly higher accuracy of this risk signature was observed by ROC curve analyses at 1 year (Figure 3H), demonstrating the sensitivity and specificity of this risk signature for OS prediction of UCEC.

Figure 3. Associations between risk signature and UCEC prognosis. Risk score distribution (A) and survival status (B) analysis of TCGA-UCEC cohort. (C) Expression level of hub lncRNAs in risk subgroups. (D) Survival curve of UCEC patients. Univariate (E) and multivariate Cox (F) regression of clinicopathological features. TimeROC (G) and ClinicalROC (H) curves to forecast overall survival of patients.

Additionally, compared to patients with TNM stage III-IV, patients with stage I-II exhibited significantly lower risk scores (p < 0.05; Figure 4A). Meanwhile, compared to patients with grade 1 or 2, obviously higher risk scores were found in patients with grade 3 (p < 0.05; Figure 4B). Furthermore, the risk signature’s value for prognosis in UCEC patients with diverse clinical features was investigated. As a result, it was revealed that there existed critical significant differences among low- and high-risk subgroups in patients younger than 60 years old (Figure 4C), over 60 years old (Figure 4D), patients with TNM stage I (Figure 4E), and patients with grade 3 (Figure 4F). In these two subgroups, all high-risk signatures displayed a significant OS disadvantage when compared with the low-risk signature. Finally, to predict the outcome of UCEC patients, we constructed a nomogram using this identified risk signature (Figure 5A), and the calibration curves at 1-, 3-, and 5-year follow-up showed that our nomogram had a substantial agreement (Figure 5B). Overall, due to the close association of the risk signature with UCEC development, our established risk signature might be a valuable tool for managing UCEC patients in clinics.

Figure 4. Associations between risk signature and clinicopathological factors. Correlations between risk scores and TNM stage (A) and grade (B). The prognosis of risk signature under the stratifications of (C, D) age ≤60 and age >60; (E) TNM stage I; and (F) grade 3.

Figure 5. Construction of nomogram. (A) Nomogram for predicting UCEC 1-, 3-, and 5-year overall survival. The red dashed line represented a sample of UCEC patient's death probability by year 1, 3, and 5. (B) Decision curve analysis of risk signature.

Functional enrichment analyses

Using GSEA and GSVA analyses, significant enrichments of the hub identified lncRNAs were identified in pathways such as cell adhesion molecules, cell cycle, and chemokine signaling (Figure 6A–6J). Meanwhile, several immune-related pathways were also enriched, such as primary immunodeficiency, natural killer cell-mediated cytotoxicity, antigen processing and presentation, and immune network for IgA production.

Figure 6. Functional enrichment analysis of hub lncRNAs. GSEA analysis of lncRNAs AL359220.1 (A), AC245884.9 (B), AC100861.1 (C), AC009097.2 (D), HOXB-AS3 (E). GSVA analysis of lncRNAs AL359220.1 (F), AC245884.9 (G), AC100861.1 (H), AC009097.2 (I), HOXB-AS3 (J).

Immune infiltration of Hub lncRNAs

As shown in Figure 7A, several immune cells were significantly differential infiltrated between UCEC and control samples by the CIBERSORT algorithm analysis. The results indicated that memory B cells, activated CD4 memory T cells, helper follicular T cells, M0 and M1 macrophages, activated dendritic cells, eosinophils, and neutrophils were significantly upregulated in UCEC samples; however, the proportions of naïve B cells, resting CD4 memory T cells, activated NK cells, M2 macrophages, and resting mast cells, were significantly increased in control samples. The correlation between hub identified lncRNAs and immune infiltration in UCEC was determined. Figure 7B shows that HOXB-AS3, AC009097.2, AL359220.1, AC100861.1, AC245884.9 were all strongly associated with the content of immune cells, indicating that all these lncRNAs may be prognostic targets for UCEC immunotherapy.

Figure 7. Immune infiltration analysis. (A) The proportion of 22 types of immune cells between normal control and UCEC samples. (B) Correlation heatmap depicting correlations between infiltrated immune cells and hub lncRNAs in UCEC.

Relationship of the risk signature with immunity, tumor stemness, and m6A-related genes

Based on the TIMER, CIBERSORT, QUANTISEQ, MCP counter, XCELL, and EPIC analyses, a close connection of the risk signature to several immune cells was detected (Figure 8A). A significant reduction of various immune cell subpopulations and functions, including aDCs, DCs, macrophages, Th2 cells, Treg, APC co-inhibition, APC co-stimulation, CCR, MHC class I, and parainflammation, was found in the patients with low-risk scores compared to those with high-risk scores (Figure 8B, 8C) (p < 0.05). To further clarify the connections of the risk signature to immune components, we evaluated the immune infiltrates associated with the promotion and suppression of tumors [15], such as wound healing (C1), INF-g dominant (C2), inflammatory (C3), and lymphocyte-depleted (C4) subtypes. Among these subtypes, a significantly lower risk score was found for the C3 subtype (Figure 8D).

Figure 8. Potential role of risk signature in UCEC immune status. (A) Heatmap for immune responses based on EPIC, XCELL, MCP counter, QUANTISEQ, CIBERSORT, and TIMER among two risk subgroups. Boxplots of scores of immune cells (B) and immune-associated functions (C) in risk subgroups. Associations between risk signature and immune infiltration subtypes (D).

For immune checkpoints, various identified immune-associated genes were significantly differentially expressed in the two risk subgroups (p < 0.05) (Figure 9A). Meanwhile, the patients with high-risk scores exhibited a dramatically low expression of several genes, including TNFRSF14, CD200, TNFRSF25, VTCN1, HHLA2, CD40LG, TNFSF14, and BTLA, but except for CD276, TNFSF9, CD80, PDCD1LG2, CD40, TNFSF4, TNFRSF8, and CD274. Moreover, we comprehensively analyzed the connection of PD-L1 loci to the risk signature. We observed a dramatically increased PD-L1 (Figure 9B) expression in high-risk score patients compared to low-risk score patients. Meanwhile, an obvious positive correlation was also found between PD-L1 (Figure 9C) expression with the calculated risk score.

Figure 9. Potential role of risk signature in immune checkpoints. (A) Expression of immune checkpoints among two risk subgroups in UCEC patients. (B) Expression levels of genes PD-L1 in risk subgroups. (C) Correlation analysis between risk score and PD-L1.

Tumor stemness, such as DNA methylation pattern and RNA stemness score, and m6A-associated genes are also key regulators of tumor progression. Compared to the subgroup with low-risk scores, the subgroup with high-risk scores exhibited dramatically decreased YTHDC2 expression and increased RBM15 expression (p < 0.05) (Figure 10A). Considering tumor stemness, a significantly positive correlation of the risk signature was observed with RNA methylation patterns (RNAss; p < 0.05) but not with DNA methylation patterns (DNAss; p > 0.05) (Figure 10B, 10C).

Figure 10. Potential role of risk signature in m6A-related genes and tumor stemness. (A) Expression of m6A-related genes among two risk subgroups in UCEC patients. Associations between risk signature and DNAss (B) and RNAss (C).

The role of lncRNA HOXB-AS3 in UCEC cells

Since the sequence of identified lncRNAs AC009097.2, AL359220.1, AC100861.1, AC245884.9 has not been clarified in NCBI database, we finally validated the expression level of HOXB-AS3 in HEC1A cell line. The lncRNA HOXB-AS3 was significantly higher expressed in HEC1A cells than in hEM15A cells (p < 0.05; Figure 11A), which was completely consistent with the bioinformatic analysis.

Figure 11. Role of lncRNA HOXB-AS3 in UCEC cells. (A) QRT-PCR analysis of SNCG. (B) The expression level of HOXB-AS3 in HOXB-AS3 inhibition model. (C) MTT assay. (D) Wound-healing assay.

To investigate the biological function of HOXB-AS3 in UCEC cells, the model of HOXB-AS3 inhibition was achieved by transfection of siRNA-HOXB-AS3 into HEC1A. As shown in Figure 11B, the expression of HOXB-AS3 was dramatically reduced by siRNA-HOXB-AS3 (p < 0.05), and no significant difference was observed between control and siRNA-NC subgroups (p > 0.05). MTT assay revealed that the proliferative ability of HEC1A was distinctly hampered by HOXB-AS3 deficiency as compared with that in victor transfected cells (p < 0.05; Figure 11C). Cell migration are another important aspect of cancer progression. Wound-healing assay results showed that HEC1A with HOXB-AS3 depletion exhibited significantly lower scratch healing rate than those in siRNA-NC and control subgroups (p < 0.05; Figure 11D).

Raw data acquisition

The normal endometrial cases and UCEC RNA sequencing datasets TCGA-UCEC (UCEC samples = 554, normal samples = 35) with reliable sources were collected from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. Samples were obtained from Homo sapiens. According to previous studies [11], five oxeiptosis-associated genes (PGAM5, KEAP1, AIFM1, NRF2, and AIRE) were screened out and applied for further analysis. All public databases in this study were searched following relevant guidelines, and no ethical approval was required from the Ethics Committee of the First People’s Hospital of Linping District.

Construction of the lncRNA Signature for prognosis

After evaluation of the connections between UCEC and oxeiptosis-associated lncRNAs and Pearson correlation analysis (|R²| > 0.2, p < 0.05), the “limma” R package was applied to identify differentially expressed lncRNAs related to oxeiptosis. Candidate lncRNAs were defined as the false discovery rate (FDR) < 0.05 and |log₂ fold change (FC)| > 1 between tumor and normal tissues. Then, the “survival” R package was applied, and the univariate Cox regression analysis was performed to select the prognostic oxeiptosis-associated lncRNAs from all lncRNAs with a cutoff p < 0.01. We selected the overlapping lncRNAs between lncRNAs related to prognosis and differentially expressed as our candidate oxeiptosis-related lncRNAs. The “VennDiagram” package was applied to visualize the results in a Venn diagram.

After that, to identify the hub lncRNA and generate the lncRNA risk signature, we integrated these selected lncRNAs into a Lasso penalized Cox regression analysis. The risk score of the hub oxeiptosis-associated lncRNA risk signature was constructed using the following formula:

$$\text{risk score}=\text{Σexplnci}\times \text{βi}$$

where explnci represents the relative expression of hub oxeiptosis-associated lncRNA i, and β is the regression coefficient. Then, UCEC patients were separated into low- and high-risk subgroups according to the median value of the constructed risk score.

Predictive value of the lncRNA signature

For investigating the distribution of these two risk subgroups, the “ggplot2” and “Rtsne” packages were employed. Based on the levels of the risk scores, the prognostic ability was compared using Cox regression and survival analyses. Then, the “timeROC” package was applied to calculate the accuracy of this risk signature for prediction. For predicting UCEC patients’ outcomes, we constructed a nomogram using the “rms” package according to the risk scores, and the decision curve analysis was conducted to evaluate the accuracy and discrimination.

Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

GSEA was applied to the gene expression matrix using the Hallmark and C7 gene sets v7.4. Enriched gene sets were used to detect KEGG pathways. Gene sets with p.adjust < 0.05 were considered significantly enriched after 1000 substitutions. GSVA was performed for each gene set and scoring. According to the GSVA score matrix, the changes at the gene-level were converted into changes at the pathway-level by the R package “GSVA”, and the potential biological functions were ultimately evaluated.

Immune and stem cell-like features and m6A correlation analysis

The relationship between the expression of hub lncRNAs and immune cells infiltration was evaluated by CIBERSORT analysis. For exploring the immune functions and comparing the infiltration of immune cells between the two subgroups, a single-sample GSEA was performed. The association of the risk score with immune infiltration subtypes was tested by two-way analysis of variance (ANOVA). The relationship of the risk signature with the immune-associated genes was determined using the potential immune checkpoints retrieved from a previous study [25]. Next, we evaluated the associations of the risk signature with PD-L1. The relationships among m6A-related genes, tumor stemness, and risk score were assessed using Spearman correlation analyses.

Cell culture

The human endometriosis cell line, hEM15A, and UCEC cell line, HEC1A, were purchased from American Type Culture Collection (ATCC) and cultured in Dulbecco’s Modified Eagle’s medium (DMEM; Gibco, Thermo Fisher Scientific, Waltham, MA, USA). The 10% FBS- (Hyclone, USA) and 1% penicillin/streptomycin- (Solarbio, China) contained medium was used for cell culture at 37°C under 5% CO₂.

Cell transfection

The short interference RNA that targets HOXB-AS3 (si-HOXB-AS3; 5′-UGCUUGUCUGGAGAUGGAGCCACUA-3′) were synthesized by GenePharma Corporation (Shanghai, China), and HEC1A cell lines were divided into three subgroups: siRNA-HOXB-AS3 (transfected with siRNA-HOXB-AS3), siRNA-NC (transfected with nonspecific scrambled), and control (cells without transfection). Cells were performed with Metafectene transfection reagent (Invitrogen, USA) according to the manufacturer’s protocols. Briefly, siRNA was formulated with Metafectene transfection reagent and added directly to the cells after diluted into culture medium. The transfection efficacy was confirmed by qRT-PCR analysis.

Cell proliferation analysis

Cells were seeded in 96-well plates at a density of 3500 cells/well. After 1 to 5 days of culture, cell viability was then assessed through MTT (3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) analysis, and the absorbance was assessed at 490 nm.

Wound-healing assay

Cells were seeded in a 6-well plate and incubated for 24 h. A wound was created by the sterile 200 μL pipette tip at cell surface. The wound closure was quantified at 0 h and 24 h after the wound was created using Image J software (NIH, Bethesda, MA, USA).

qRT-PCR analysis

Trizol reagent (Invitrogen, USA) was used to isolate RNA from cell lines. The RNeasy mini kit (QIAGEN, USA) was used to isolate RNA. Then, the RNA was reverse transcribed using a Transcriptor First-strand cDNA synthesis kit (Roche, Switzerland). SYBR green master mix was used to perform qRT-PCR using an Applied Biosystems 7500 Real-Time Cycler (Applied Biosystems, USA). Gene expression was measured using the 2^−ΔΔCT method. GAPDH (Forward: CTGCCTCGATGGGTGGAGTC; Reverse: GAGTTAAAAGCAGCCCTGGTG) was used as a normalization control. Primer sequences for HOXB-AS3 as follows: Forward: TGCTTGTCTGGAGATGGAGC; Reverse: GATAAGAGCGATGAGGCGCT.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.