Two lncRNA signatures with cuproptosis as a novel prognostic model and clinicopathological value for endometrioid endometrial adenocarcinoma

Publicly available datasets and preprocessing

We downloaded the mRNA expression level by RNA sequencing (RNA-Seq), somatic mutations and clinical phenotype of 516 endometrial cancer patients from TCGA database(https://portal.gdc.cancer.gov/) [11]. The RNA sequencing data was then screened against the clinical data collected for patients with endometrial cancer based on their histological subtype. Endometrioid endometrial adenocarcinoma (EEA) is a type of endometrial cancer that originates from the endometrial glands and is characterized by estrogen-dependent growth and favorable prognosis. EEA variant data were retained. RNA-Seq data were normalized using the R software package limma [12]. According to the annotation file for the platform, probes were converted into gene symbols. If several probes corresponded to the same gene symbol, the probe with the greatest value was used. The identification and removal of outliers was conducted using boxplots and principal components analyses (PCAs) [13]. Finally, the expression and clinical data of 318 patients (TCGA-UCEC) were used for the study analysis (Table 1).

MAF was originally developed by The Cancer Genome Atlas (TCGA) as a standardized format for storing somatic variants detected in cancer samples, we analyze and visualize MAF Files by the MAFtools R package [14]. Tumor mutational burden (TMB) is defined as the total number of somatic, coding, base substitutions, and indels (insertions and deletions) per megabase of the genome examined, which were extracted and estimated by Perl scripts. Furthermore, the data about stemness scores, MSI status, and infiltration measure were summarized from a number of research papers, so as to provide a more comprehensive analysis. Stemness scores can be accessed at https://bioinformaticsfmrp.github.io/PanCanStem_Web [15]. The data on MSI status were acquired using the R package “TCGAbiolinks”. Infiltration measure data (ESTIMATE, CIBERSORT etc.) were obtained by an R package ‘IOBR’ [16].

Coexpression analysis between cuproptosis-related genes and lncRNAs

By reviewing the literature, we obtained 19 key genes including NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, GCSH, DLST associated with copper death signaling [17]. Then, we performed coexpression analysis between cuproptosis-related genes and lncRNAs using the “correplot” package in R [12]. Co-expressed genes can be identified by using the Pearson correlation coefficient with a threshold of at least 0.4 and a p-value lower than 0.001. R package “ggalluvial” was used to visualize the results. Finally, 941 cuproptosis-related lncRNAs were found.

Identification of cuproptosis-associated clusters and survival analysis

After performing univariate Cox regression analysis of CRLs, we retained 11 lncRNAs associated with cuproptosis. ConsensusClusterPlus R package was used to analyze consistency [18]. The maximum number of clusters was set to 6, and 100 samples were drawn for clusterAlg = “hc”, and innerLinkage = “average”. Then, we performed a survival analysis for the clusters.

Prognostic model construction and Nomogram prediction model

A risk signature for predicting the prognosis of EEA patients was constructed based on the lncRNAs related to cuproptosis. Least absolute shrinkage and selection operator (LASSO) regression was used to select predictors to avoid overfitting [19]. Multivariate Cox regression analysis was used to determine the candidate genes for final inclusion in the risk model. A risk signature for predicting the prognosis of EEA patients was constructed based on the lncRNAs related to cuproptosis. The mathematical formula for calculating the risk score was as follows:

$$\text{riskscore}={\displaystyle \sum _{i=1}^{N}Expression\_i\times coefficient\_i}$$

Patients were divided into high and low risk groups using the median risk score as the cutoff value.

Nomograms are widely used for cancer prognosis mainly because of their ability to reduce statistical predictive models to a single numerical estimate of event probabilities [20]. Based on the results of LASSO and multivariate COX regressions, we developed a nomogram by using the R packages, “regplot”, “RMS” and “survival”. Calibration curves were used to assess the accuracy of the nomogram.

Pathway enrichment analysis and gene set enrichment analysis

In this study, differentially expressed genes (DEGs) between high-risk and low-risk groups were identified using the R package “limma.” The threshold value was set to log2 |fold change|>1 and p-value < 0.05. We applied Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to elucidate molecular functions and key signaling pathways. At the same time, gene set enrichment analysis of EEA patients were performed to identify an association with disease phenotypes and pathways through R package clusterProfiler [21].

Tumor microenvironment and immune cell infiltration profile

According to the proportion of immune cells and stromal cells in the tumor, an estimation algorithm was used to calculate the stromal score, immune score, estimated score and tumor purity. We applied TIMER, MCP-counter, CIBERSORT, QUANTISEQ, and xCell to calculate the proportion and abundance of tumor-infiltrating immune cells (TICs) in EEA [22]. Using the single-sample gene set enrichment analysis (ssGSEA) algorithm, we evaluated tumor-infiltrating immune cell infiltration and function [23]. An algorithm called Tumor Immune Dysfunction and Exclusion (TIDE) was used to predict the response to immunotherapy [24].

Chemotherapeutic response and construction of the ceRNA network

The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity [25, 26]. The efficacy of chemotherapeutic drugs in EEA patients was predicted through the Genomics of Drug Sensitivity in Cancer (GDSC) database [27]. Drug sensitivity and gene expression profiling data of cancer cell lines in GDSC are integrated for investigation. The expression of each risk group in the gene set was performed by Spearman correlation analysis with the small molecule/drug sensitivity (IC50). The half maximal inhibitory concentration (IC50) was calculated by Prophetic software in R package.

According to the competitive endogenous RNA (ceRNA) theory, we constructed a network of lncRNAs, miRNAs, and mRNAs [28]. StarBase V2.0 was used to analyze the interaction between lncRNAs and miRNAs [29].

Cell culture and quantitative real-time PCR analysis

HEC-1A is a human endometrial cancer cell line obtained from the Cell Bank of the Chinese Academy of Sciences. The HHUA cells were generously donated by the Department of Obstetrics and Gynecology, the First Affiliated Hospital of Shantou University Medical College. The cells were cultured in the presence of 10 μM Cu2+ and 10 nM Elesclomol for 48 hours, followed by microscopic observation and imaging. Cell proliferation was measured using the CCK8 kit (Beyotime, Shanghai, China). Subsequently, cellular RNA was extracted for further analysis. A control group was included, in which cells were cultured without the addition of Cu²⁺ and Elesclomol. TRIzol reagent was used to extract total RNA. Total RNA was reverse transcribed using PrimeScript reverse transcription reagent (Takara, Otsu, Shiga, Japan). In accordance with the manufacturer’s protocol, miRNA reverse transcription was performed using miRNA stem loop reverse transcription kit (Shanghai Sangong Biological Company, China). Quantitative PCR analysis was carried out using the 7500 Fast Real-Time PCR System instrument with TB Green Premix Ex Taq II (Takara, Japan) according to the manufacturer’s protocol. In each sample, the experiment was repeated three times and 2-Ct was used to calculate the RNA expression. Table 2 shows the primer sequences used in this study.

Lentiviral transfections and Western blotting

To establish human endometrial cancer cell lines with overexpression of ACOXL-AS1, we used ACOXL-AS1-expressing lentiviral vector (OE-ACOXL-AS1) and negative control lentiviral vector provided by Akey Biotechnology Co. LTD (Guangzhou, China). Human endometrial cancer cells were seeded in each well of a 6-well plate at a density of 1 × 105 cells per well. When the cell confluence reached 60% in each well, we added lentivirus (MOI = 10) to the cell culture medium. Then, we selected stable ACOXL-AS1-overexpressing cells with puromycin (3 μg/ml). Finally, we measured the expression level of ACOXL-AS1 in stable cells by qRT-PCR.

Cell lysis was performed with RIPA cell lysate (absin) and the lysates were centrifuged at 14000 g/min for 5 min. The supernatants were collected, and the protein concentration was determined by a BCA protein assay kit (Thermo Scientific). Equal amounts of protein were separated by SDS-PAGE and transferred to PVDF membranes. The membranes were incubated with the following primary antibodies: rabbit anti-MTF1 (1:1000; ab236401; Abcam), mouse anti-β-actin (1:1000; TA-09; ZSGB-BIO). Fluorescent antibody Anti-mouse/rabbit 800 was used as the secondary antibody. The signals were detected by Bio-Rad ChemiDoc™ MP Chemiluminescence/fluorescence imaging system (Bio-Rad) and quantified by ImageJ software (1.4.3.67, NIH).

Statistical analysis

R (version 4.2.1) software was used for statistical analysis and visualization of the data. The following packages were used: “Maftools,” “ggplot2,” “limma,” “survminer,” “dplyr,” “plyr,” “survivalROC,” “clusterProfiler,” “ggplotify,” “cowplot,” “Hmisc,” “gridExtra,” “GSVA,” “corrplot,” “VennDiagram,” “pheatmap” and “tidyverse,”. When comparing two groups, the Wilcoxon test was used, and when comparing three groups, the ANOVA test was used. Prognosis was assessed using the Kaplan–Meier method and the Pearson method was used for correlation. We considered a p-value of 0.05 as statistically significant (^*p < 0.05; ^**p < 0.01; ^***p < 0.001).

Data and code availability statements

According to the TCGA project (https://portal.gdc.cancer.gov/) policies, all public access to these databases will be open. All code used in this manuscript are available at https://github.com/w28924461701/Cuproptosis-in-UCEC/.

Identification of CRLs and consensus clustering analysis

In the EEA cohort within TCGA, we obtained nineteen cuproptosis-related genes and 16,876 lncRNAs that were annotated by NCBI GenBank, Ensembl, and GENCODE. By co-expression analysis, 941 cuproptosis-related lncRNAs were found, and associations are visualized in the Sankey plot (Figure 1A). Eleven CRLs associated with prognosis were mined by univariate Cox regression analysis (Figure 1B, P < 0.05). The least absolute shrinkage and selection operator (LASSO) is an efficient gene selection method [30]. A total of 11 lncRNAs were selected to fit a LASSO regression model, the next step was to find the most appropriate values for λ (=0.00659) using 10-fold cross-validation (Figure 1C). Finally, two cuproptosis-related lncRNAs (AL512353.1, ACOXL-AS1) were identified by LASSO and Multivariate Cox regression analysis (Figure 1C, 1D). Cytoscape software was employed to visualize associations between lncRNAs and cuproptosis-related gene (Figure 1E).

Figure 1. Consensus clustering analysis of cuproptosis-related lncRNAs and LASSO regression model. (A) Coexpression analysis of IncRNA and cuproptosis. (B) Univariate Cox regression analysis of cuproptosis-related lncRNAs. (C) LASSO regression model for the prognostic analysis of the lncRNA signature. (D) Forest plot visualing the results of multivariate Cox regression. (E) Correlation between two lncRNAs and 19 cuproptosis-related gene. Network visualization and analysis of lncRNAs and cuproptosis-related gene with Cytoscape. (F) Consensus matrix heatmap defining three clusters (k = 4) and their correlation area. (G, H) Cumulative distribution function (CDF) when k = 2–6. (I) Kaplan–Meier curve of the 4 clusters of patients with EEA in the TCGA cohort.

Cluster analyses were conducted using the eleven prognosis-related lncRNAs. When EEA patients were divided into four subgroups, the clusters could be seen, and the subgroups showed good internal consistency and stability (Figure 1F–1H). Survival analysis revealed the prognostic significance of the clustering (Figure 1I, P = 0.025).

Construction and validation of cuproptosis-related risk score and consensus clustering analysis

To further understand the prognostic role of the two cuproptosis-related lncRNAs, we randomly divided the dataset into training and validation groups. The training dataset was used to train a risk model, while the testing dataset was used to evaluate the performance of the model. A multivariate Cox regression analysis was then applied to the remaining genes to establish a prognostic signature in the training dataset. Regression coefficient in multivariate Cox regression analysis were derived from the training data, and were the parameters further refined using the combined dataset after passing the testing dataset. The risk score (RS) equation was defined as RS = (AL512353.1 × 1.264) + (ACOXL-AS1 × 1.658). According to the mathematical formula for calculating the risk score, patients were classified into a high risk or low risk group based on the median value of RS.

Kaplan-Meier survival analysis in the training, validation and all samples set showed that the overall survival (OS) of patients in the high-risk group predicted poor survival (Figure 2A–2C, P < 0.05). The model yielded better accuracy and calibrated survival estimates in predicting 1-, 3-, and 5-year survival (AUC: 0.703 to 0.798, Figure 2D). The layout of risk scores and survival status of EEA patients is shown in Figure 2E, 2F. In addition, the expression levels of the two cuproptosis-related lncRNAs between high- and low-risk groups are also displayed.

Figure 2. Evaluation of the established cuproptosis-related lncRNAs signature in the training and test sets. (A–C) Kaplan-Meier survival analysis in the training (A), validation (B) and all samples (C) sets. (D) ROC curves for risk scores at 1, 2, and 3 years. (E, F) Expression heat map, risk score distribution, and relapse status.

Relationships between cuproptosis-related risk and clinicopathological features

The relationship between risk and clinicopathologic characteristics in EEA patients was further investigated. When comparing high-risk and low-risk groups with current clinical variables, we can use the chi-squared test to measure the association between categorical variables and the Wilcoxon rank-sum test to compare continuous variables. Risk classification was correlated with molecular subtypes, histologic grade, and myometrial invasion (p < 0.05), but other features did not obvious relevance (Figure 3A–3E, Table 3). The risk score was associated with molecular subtypes (p < 0.01), cancer status (p < 0.01) and age (p = 0.03), in EEA patients, POLE hypermutation and MSI subtypes have low risk values, while high CNV is associated with high-risk values (Figure 3A). For lymph node metastasis, there were no obvious differences in risk score between patients with lymph node metastasis and that without lymph node metastasis (Figure 3B–3D, p = 0.07). At the same time, in EAA patients in early-stage cancer (stage I, II), we also found that the patients in the predicted high-risk groups had shorter survival time than those in the low-risk groups (Figure 3E, p = 0.032).

Figure 3. Clinicopathological analysis and prognostic nomogram based on cuproptosis-related lncRNAs signature. (A–D) Correlation between risk score and molecular subtypes (A), cancer status (B), age (C), and lymph node metastasis (D). (E) Correlation between risk score and survival probability in patients at stage I–II. (F, G) Univariate and multivariate Cox regression analysis. (H) Nomogram for predicting probabilities of EEA patient 1-, 3-, 5-year OS in the TCGA cohort. (I) ROC curves showing the predictive efficiency of the nomogram. (J) Calibration plots of 1-, 3-, and 5-year OS predicted by the nomograms.

Development of a nomogram following univariate and multivariate Cox analysis

To investigate the influence of clinicopathological features on prognosis of EEA patients, univariate and multivariate Cox proportional hazards analyses were employed, and results visualized as a forest plot. Univariate Cox analysis showed that risk Score, distant metastasis and necrosis were related to OS (Figure 3F P < 0.05). Then, multivariate Cox regression analyses demonstrated that risk Score and necrosis were also independently associated with OS (Figure 3G P < 0.05).

A nomogram model including all seven clinicopathologic factors was constructed to predict survival probability by adding all the points associated with the seven factors (Figure 3H). The higher the score, the greater the likelihood of death. The result showed that risk, age, lymph node metastasis, necrosis and cancer status were correlated with OS (p < 0.05). To validate reliability of the nomogram, the AUC values were calculated. Internal validation showed an AUC value of 0.816 (Figure 3I), a total C-index value for predicting OS of 0.907, and a C-index value of risk scores of 0.782. Calibration curves of the 1-, 3- and 5-year overall survival showed that the predicted OS was in good agreement with the observed OS (Figure 3J).

Gene Ontology, KEGG pathway and gene Set enrichment analysis

Identification of differentially expressed genes (DEGs) can be used to gain mechanistic insights from diseases [31]. When R package limma was used to analyze differentially-expressed genes between the high-risk and low-risk groups, 354 genes (228 up-regulated and 126 down-regulated genes) showed significant changes (Figure 4A). The results of Gene Ontology (GO) analysis indicated that the differentially expressed genes (DEGs) were significantly enriched in various biological processes (BP), molecular functions (MF) and cellular components (CC). Specifically, the enriched BP terms included the production of molecular mediators of the immune response, while the enriched MF terms involved immunoglobulin complex formation, antigen binding, immunoglobulin receptor binding, and receptor-ligand activity (Figure 4B–4D). Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, the differentially expressed genes (DEGs) were associated with multiple pathways, such as neuroactive ligand-receptor interactions, ECM receptor interactions, and IL-17 signaling pathways (Figure 4E–4G). Furthermore, GSEA identified that biological pathways including positive regulation of triglyceride lipase activity and neuropeptide signaling pathway were differentially enriched in the high and low risk groups (Figure 4H).

Figure 4. Biological function, and mechanism analysis. (A) Volcano plot of differentially-expressed genes (DEGs) between the high-risk group and low-risk group. (B) GO Circle plot (red, upregulated genes; blue, downregulated genes) show that Top 15 enriched GO terms in the BP category. (C) and (D) Plot of the enriched GO terms Go enrichment analysis for associated mRNAs with risk grouping. Y-axis represents the enriched GO terms; X-axis (C) represents the amount of the related mRNAs enriched in GO terms; X-axis (D) represents the ratio of the related mRNAs enriched in GO terms. (E) represents the amount of the related mRNAs enriched in KEGG pathways, X-axis. (F) represents the ratio of the related mRNAs enriched in KEGG pathways. (G) DEGs associated with the significant KEGG pathway. (H) GSEA showing the top six most significantly enriched signaling pathways.

Immune infiltration analysis and Immunotherapy

The ESTIMATE algorithm was applied to infer tumor purity, immune score, and stromal score, which represent the level of immune cell infiltration in the tumor [32]. There was a correlation between the risk scores and both immune score, and stromal score (Figure 5A–5E, p < 0.01). The heat map showed infiltration of immune cells in each tumor sample, as determined by TIMER, CIBERSORT, quanTIseq, MCPcounter, XCell and Epic software (Figure 5F, 5G). The CIBERSORT algorithm was used to calculate the association of 22 immune cell proportions (Figure 5H). The CIBERSORT algorithm was also employed to identify the proportion of immune cells between high-risk and low-risk groups and showed that the signature in high-risk patients was associated with decreased CD8+ T lymphocyte infiltration, mast cells resting (Figure 5I). Meanwhile, ssGSEA was used to evaluate all 13 immune-related functions in the high- and low-risk groups, and the result showed obviously differences in immune check point, T cell co-stimulation, cytolytic activity, promoting of inflammation, and APC co stimulation, indicating a strong relationship between cuproptosis-related lncRNAs and immunity (Figure 5J).

Figure 5. Immune profile of cuproptosis lncRNAs based on EEA groups. (A–E), Level of immune cell infiltration in the tumor, determined by ESTIMATE. (F, G) Infiltration situation of immune cells in each tumor sample. (H) Calculate of the association of 22 immune cell proportions. (I) CIBERSORT algorithm identifying the proportion of immune cells. (J) Scores of 13 immune-related functions by ssGSEA, ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Immune checkpoint inhibition is an immunotherapy method that blocks the binding of immune checkpoint proteins to chaperone proteins [33]. Gene expression analysis of immune checkpoint between the risk groups showed that PDCD1 and TIGIT were associated with risk group and tended to be expressed at low levels in the high-risk group (Supplementary Figure 1A). The scatter patterns of the top 15 most frequently mutated genes can be seen by looking at somatic variants, and it was found that the high-risk group had a lower rate of genetic alterations (94% vs. 97.35%, Supplementary Figure 1B, 1C) than the low-risk group. There were lower TMB values in the high-risk group, but higher tumor stemness and TIDE scores in the high-risk group. (Supplementary Figure 1D–1I).

Prediction of chemotherapeutic response

By identifying molecular signatures of cancer, the GDSC database helps predict tumor’s response to antitumor therapy [34]. A large difference in IC50 was found between high-risk and low-risk groups for six chemotherapy agents, such as AKT inhibitors (Figure 6A–6F). Simultaneously, there was a correlation between six drugs and risk score (Figure 6G–6L, p < 0.001).

Figure 6. Chemotherapy drug response prediction for EEA in high- and low-risk EEA patients. (A) BMS345541. (B) NPK76-II-72-1. (C) AP-24534. (D) AKT inhibitor VIII. (E) BX-912. (F) XL-184. (G–L) Correlation analysis between risk group and 6 chemotherapy drugs.

Construction of LncRNA-miRNA-mRNA networks

Using starBase query results, miR-421, miR-483-3p, and miR-3200-5p can bind to their lncRNA ACOXL-AS1 targets (Figure 7A). Then, we conducted a differential comparison analysis of genes related to copper death, and the heatmap showed that 13 genes were differentially expressed between normal and tumor group in TCGA-EEA (Figure 7B). The online database miRDB provides functional annotation and predictions for miRNA targets [35]. Through the miRDB database, we found 3 miRNAs targeting 784 unique genes, and intersecting with 13 copper death-related genes in TCGA. Two of these genes, MTF1 and GLS were related to miR-421 and miR-483-3p, respectively (Figure 7C). At the same time, we also found that miR-421 was highly expressed in tumors (Figure 7D, p < 0.01). Thus, ACOXL-AS1 can serve as an endogenous “sponge” to regulate the expression of MTF1 by miR-421 (Figure 7E).

Figure 7. Construction of LncRNA-miRNA-mRNA networks. (A) StarBase V2.0 was used to analyze the interaction between lncRNAs and miRNAs. (B) Heatmap indicating cuproptosis-related genes show differential expression changes in normal and EEA tumor groups. (C) Venn diagram showing the common genes between TCGA and miRDB target genes. (D) miR-421 expression in normal or tumor group. (E) Schematic diagram for ACOXL-AS1-miR-421-MTF1 regulation.

Cuproptosis inhibited the cell proliferation and influenced on ACOXL-AS1, hsa-mir-421, and MTF1 expression in endometrial cancer cells

To investigate the impact of cuproptosis on the growth of endometrial cancer cells, we established a copper-induced cell death model. When HHUA and HEC-1A, two endometrial cancer cell lines, were treated with Cu2+ plus Elesclomol, their growth was significantly inhibited in comparison to the control group (p < 0.001, Figure 8A, 8B).

Figure 8. Cuprotosis cell model and expression of ACOXL-AS1, mir-421 and MTF1 in endometrial cancer cells. (A) Colony formation ability under differential conditions as determined by a colony formation assay. (B) Proliferation of endometrial cancer cells is detected by the Cell Counting Kit-8 assay. (C, D) Relative expression of ACOXL-AS1, mir-421 and MTF1 in endometrial cancer cells under different conditions. (E) overexpress the ACOXL-AS1 gene in HEC-1A and HHUA cells, while the “negative control lentiviral vector” is used as a reference (F) The impact of ACOXL-AS1 overexpression on the expression levels of MTF1 and miR-421, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001.

In addition, the Cu2+ plus Elesclomol-treated group showed significantly higher gene expression levels for ACOXL-AS1 and hsa-mir-421 (p < 0.01, Figure 8C), while low expression levels were observed for MTF1 gene and protein in two endometrial cancer cell lines (p < 0.01, Figure 8C, 8D).

Effect of ACOXL-AS1 overexpression on the expression of hsa-mir-421 and MTF1

To study the impact of ACOXL-AS1 on the expression of hsa-mir-421 and MTF1 in human endometrial cancer cell lines, two human endometrial cancer cell lines, HEC-1A and HHUA, were used in the study. Lentiviruses encoding ACOXL-AS1 cDNA (referred to as “ACOXL-AS1-OE”) and an empty vector were used to establish stable cells. The selection of stable cells was achieved using puromycin. ACOXL-AS1 expression in the stable cells was measured using quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR). The qRT-PCR results showed that the expression of ACOXL-AS1 increased significantly in the ACOXL-AS1-overexpressed (“ACOXL-AS1-OE”) cells (Figure 8E). ACOXL-AS1-overexpressing cells showed increased expression of hsa-mir-421 and reduced MTF1 expression compared with controls (Figure 8F).