Regular Contribution
The Prognosis and Immunotherapy Prediction Model of Ovarian Serous Cystadenocarcinoma Patient was Constructed Based on Cuproptosis-Related LncRNA
2024 年 262 巻 2 号 p. 63-74
詳細
2024 年 262 巻 2 号 p. 63-74
Cuproptosis can serve as potential prognostic predictors in patients with cancer. However, the role of this relationship in ovarian serous cystadenocarcinoma (OV) remains unclear. 376 OV tumor samples were obtained from the Cancer Genome Atlas (TCGA) database, and long non-coding RNAs (lncRNAs) related to cuproptosis were obtained through correlation analysis. The risk assessment model was further constructed by univariate Cox regression analysis and LASSO Cox regression. Bioinformatics was used to analyze the regulatory effect of relevant risk assessment models on tumor mutational burden (TMB) and immune microenvironment. We obtained 5 lncRNAs (AC025287.2, AC092718.4, AC112721.2, LINC00996, and LINC01639) and incorporated them into the Cox proportional hazards model. Kaplan-Meier (KM) curve analysis of the prognosis found that the high-risk group was associated with a poorer prognosis. The receiver operating characteristic (ROC) curve showed stronger predictive power compared to other clinicopathological features. Immune infiltration analysis showed that high-risk scores were inversely correlated with CD8+ T cells, CD4+ T cells, macrophages, NK cells, and B cells. Functional enrichment analysis found that they may act via the extracellular matrix (ECM)-interacting proteins and other pathways. We successfully constructed a reliable cuproptosis-related lncRNA model for the prognosis of OV.
Ovarian cancer is one of the common malignant tumors of female reproductive system, and its mortality rate ranks first in gynecological tumors (Barani et al. 2021; Guo et al. 2021). The 5-year survival rate of advanced stage patients is only 30%-40%, while that in early stage patients can reach 70% (Carey et al. 2021). Therefore, finding biomarkers or efficient predictive models that are conducive to early screening will help improve the prognosis and survival of patients. Ovarian serous cystadenocarcinoma (OV) is the most common histological subtype of ovarian cancer, accounting for approximately 75% to 80% of all ovarian cancers (Kuroki and Guntupalli 2020; Penny 2020; Yang et al. 2020). At present, the discussion on OV targeted therapy has gradually become a research hotspot (Nowak and Klink 2020). The discovered therapeutic targets of OV mainly include PD-1, P53, HMGA2, CA125, SPEC-2 and PARP-1, etc. (Morand et al. 2021). Although the clinical role of these therapeutic targets has been widely explored, the effect is very limited, so the treatment of OV remains a difficult problem for clinicians. Therefore, deepening the understanding of molecular mechanisms and exploring new therapeutic targets are crucial for OV therapy.
Long non-coding RNA (lncRNA) is an important molecule in the non-coding RNA family with a length of more than 200 nt (Cao et al. 2019). Although it is not involved in encoding proteins, it has complex biological functions and can participate in the regulation of many pathophysiological processes in the body (McCabe and Rasmussen 2021). The study of lncRNA expression profile provides evidence for its function and role in ovarian cancer. A study by lncRNA expression profiling found that among 4,956 lncRNAs, 583 lncRNAs were up-regulated and 578 lncRNAs were down-regulated in ovarian cancer high metastatic cells compared with ovarian cancer parental cells, indicating that lncRNAs are involved in epithelial ovarian cancer metastasis (Braga et al. 2020). Furthermore, a recent study showed that 115 lncRNAs with significant changes, most of which are predicted to promote malignancy (Wu et al. 2020; Xing et al. 2021). These results suggest that lncRNA may be an efficient biomarker for tumor screening. However, there is often a partial bias in the prediction process of a single lncRNA, so finding an efficient composite model may be an important way to screen for OV markers.
A new study published in science was the first to propose a new copper-dependent method of cell death called “Cuproptosis” (Tsvetkov et al. 2022). Studies have shown that it differs from known methods of cell death such as apoptosis, pyroptosis or necroptosis, but rather a metal ion-induced regulation of cell death similar to zinc death and ferroptosis (Chen 2022; Li et al. 2022b). Recent studies have shown that overload or deficiency of copper ions can lead to a variety of diseases, such as age-related neurodegenerative diseases such as Alzheimer’s disease, inherited disorders of copper metabolism such as Wilson’s disease, and chronic liver and kidney diseases such as renal fibrosis (Cobine and Brady 2022; Tang et al. 2022). Studies have shown that the content of copper ions in serum and solid tumor tissues of tumor patients is significantly higher than that of healthy people, and the increase of copper ions can promote the proliferation and metastasis of tumor cells and promote angiogenesis (Li et al. 2022a). Copper complexes are coordination compounds formed by the combination of copper ions and some ligands with coordination bonds, which can interact with DNA and proteins to induce the generation of oxygen free radicals and destroy phosphodiester bonds, thus promoting the hydrolysis of DNA (Bao et al. 2022).
Studies have been conducted to explore the integrative analysis between lncRNA signatures associated with cuproptosis and cancer. For example, gene signatures associated with cuproptosis can serve as potential prognostic predictors in patients with clear cell carcinoma and may provide information for cancer therapy (Xu et al. 2022; Yang et al. 2022a; Yang et al. 2022b). However, the role of this relationship in OV and its regulatory role in the tumor microenvironment remain unclear. Therefore, this study attempted to construct a risk assessment model for the evaluation of patient prognosis and immunotherapy by analyzing the cuproptosis-related lncRNAs. These data may provide new horizons for the development of prognostic models and targeted therapy screening for OV patients.
The Cancer Gene Atlas (TCGA) database (https://portal.gdc.cancer.gov/) was used to obtain the clinical information of OV patients and the expression information of cuproptosis, lncRNA and other related genes. The clinical information for all groups is shown in Table 1. We identified and retrieved cuproptosis genes in TCGA expression matrix. In this study, the limma package was used to analyze the co-expression correlation between lncRNA and copper cytopenia-related gene expression profiles, Correlation coefficient |r| > 0.4, p < 0.001. The tumor mutational burden (TMB) was calculated by downloading the tumor mutation data using perl software, and the TMB in the samples was visualized on the RAWGraphs online website (https://www.rawgraphs.io).
Characteristic of ovarian serous cystadenocarcinoma (OV) patients.
The TCGA-OV samples were randomly divided into training cohort and test cohort in a 1:1 ratio. A lncRNA model associated with cuproptosis was constructed on the training set using Lasso-Cox regression analysis, and the model was validated using the entire set and test set. After conducting a single-factor Cox regression analysis, a Lasso-Cox regression was performed to eliminate the issue of multicollinearity between the factors. This step helped to identify lncRNAs that exhibited statistical significance. Next, a screening process was undertaken to identify lncRNAs that were associated with cuproptosis-related mortality. Using the optimal model parameters, features were constructed, and a risk score was calculated. High and low risk groups were distinguished by the median risk score. For Lasso-Cox regression, we use the R package “glmnet” to implement this step. Finally, a lncRNA risk model associated with cuproptosis was established. Risk Score = coef (lncRNA) × expr (lncRNA1) + coef (lncRNA 2) × expr (lncRNA 2) + ... + coef (lncRNA n) × expr (lncRNA n).
Construction of nomogramIn this study, a nomogram was established to predict the survival rate of OV patients by clinical features including information such as age, sex, grade, and stage. In addition, the agreement between predicted and actual survival was assessed by drawing a calibration curve.
Evaluation of the validity of the modelThe “survival” package was used to calculate the overall survival (OS) and progression-free survival (PFS) of patients with OV in different risk groups to assess the independent prognostic value of risk-predicting features. We used the survival receiver operating characteristic curve (ROC) package to compute the 1-, 3-, and 5-year areas of the features under the ROC curve for training, testing, and all groups.
Gene functional enrichment analysisGene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed, including cuproptosis-related lncRNAs, using the clustering “Profiler” package.
Tumor immune microenvironment analysis (TME analysis)Using the downloaded TCGA data, the single correlation map was implemented by the “ggstatsplot” package, and the multi-gene correlation map was displayed by “pheatmap” package.
As shown in Fig. 1, the 381 OV patients (376 patients with complete information) were included in this study. In addition, we extracted the expression data of 16,773 lncRNAs in the database and extracted 19 cuproptosis-related genes for subsequent analysis.
First, we used the Gene Expression Profiling Interactive Analysis (GEPIA) database to explore the expression levels of these 19 cuproptosis-related genes in OV tissues (Fig. 2A). In addition, correlation analysis showed that there was a positive correlation trend among them (Fig. 2B). We further screened out 113 lncRNAs that were significantly associated with cuproptosis using the r-values and p-values (Fig. 2C). Then, the significant lncRNAs in the univariate Cox analysis were screened and included in the LASSO-Cox regression analysis (Supplementary Fig. S1). Finally, we obtained 5 lncRNAs and incorporated them into the Cox proportional hazards model (Fig. 2D). The risk score formula: Risk score = Exp (LINC00996) × (−1.58) + Exp (AC025287.2) × (−1.02) + Exp (AC112721.2) × (0.54) + Exp (AC092718.4) × (−0.29) + Exp (LINC01639) × (0.91).
Schematic diagram of the analysis process.
TCGA, The Cancer Gene Atlas; KEGG, Kyoto Encyclopedia of Genes and Genomes.
Identification of cuproptosis-related lncRNAs.
A. Expression of 19 cuproptosis genes in ovarian serous cystadenocarcinoma (OV) tissues in the TCGA database. B. The correlation between each member of cuproptosis genes in OV. C. Sankey diagram showing the relationship between cuproptosis genes and lncRNAs. D. Corrplot displayed the relationship between the 5 cuproptosis-related lncRNAs with the 19 cuproptosis genes.
As mentioned above, we divided the training set and the test set in a 1:1 manner and there was no statistical difference in general indicators between the two sets (n = 188, Table 1). First, we verified the expression of these five lncRNAs in the two datasets, respectively (Fig. 3A, B). Further analysis of the distribution of risk scores and the survival status of patients showed that these 5 lncRNAs shortened the survival time and increased the number of deaths with increasing risk scores (Fig. 3C, E). These results show the same trend in the training and test sets. This suggests that our model is reliable in prognostic judgment (Fig. 3D, F).
Evaluation of lncRNAs in training and testing cohorts.
A, B. The heatmap was used to display the expression levels of lncRNAs in relevant samples in different risk groups. C, D. The risk score distribution map was used to show the distribution of patients’ survival in different risk groups. E, F. The scatter plot was used to depicted the relationship between the survival status and the risk score in ovarian serous cystadenocarcinoma (OV) patients.
Further, we used the Kaplan-Meier (KM) curve found that the high-risk group was associated with poorer prognosis. We combined the training and testing cohorts to draw survival curves and found that the prognostic survival rate of patients in the low-risk group was significantly better than that in the high-risk group (Fig. 4A). Among them, PFS refers to the survival period of patients with OV without further deterioration after treatment. Unlike the overall survival period, which uses death as the outcome indicator, the outcome indicator of the survival analysis of PFS is disease progression or progression. The corresponding improvement is no deterioration and no death, which can be considered clinically beneficial. In addition, we used ROC curves to verify the evaluation accuracy of the risk score and found that the areas under the curve (AUC) were 0.620 (1-year), 0.629 (3-year), and 0.654 (5-year). Further analysis found that the risk scoring model showed stronger predictive power compared to other clinicopathological features (AUC = 0.654, Fig. 4B). Finally, the calibration curve shows well agreement between the nomogram and the predicted results (Fig. 4C).
The Independent prognostic test was performed on 376 samples in the whole set, and it was found that the hazard ratio (HR) value of the risk model in the univariate survival analysis was 1.388 (p < 0.001), and the HR value in the multivariate survival analysis was 1.361 (p < 0.001) (Fig. 5A). Further, we explored the prognosis of OV patients with different age and tumor grades by risk score. The KM curve data showed that, except for G1-2 stage, the low-risk group showed better prognosis in different clinical categories. Furthermore, although the differences in the G1-2 data were not significant, similar trends were still seen (Fig. 5B). Finally, the principal component analysis (PCA) further verified the high degree of heterogeneity between high-risk and low-risk groups. This result is further confirmation of the accuracy of the model (Fig. 5C).
Prognostic impact of cuproptosis-related lncRNAs at different risks on ovarian serous cystadenocarcinoma (OV).
A. KM curves were used to predict the effect of high and low risk groups on overall survival (OS) and progression-free survival (PFS). B. The ROC curve was used to judge the prognosis and diagnosis of risk score at 1 year, 3 years and 5 years and the diagnosis of risk score, age and tumor grade. C. Calibration curve showed the consistency between the nomogram prediction and the actual survival.
The role of cuproptosis-related lncRNA in the prognostic evaluation of ovarian serous cystadenocarcinoma (OV) patients.
A. KM curves were used to predict the effect of high and low risk groups on overall survival (OS) and progression-free survival (PFS). B. The KM curve was used to verify the prognosis of high and low risk groups in terms of age and tumor grade. C. PCA analysis was used to detect the classification of different genes.
Current research shows that TMB was emerging as a potential biomarker. This study firstly shows the distribution of TMB scores for OV samples (Supplementary Fig. S2A). Further, this study looked at mutations in high-risk and low-risk groups and found that for most genes, the high-risk group had a higher mutation frequency than the low-risk group (Fig. 6A, B). Finally, we found that the higher the risk or the higher the incidence of TMB, the worse the OS level was in both TMB and the prognosis in the risk score (Fig. 6C, D).
Tumor mutational burden (TMB) and related prognosis in high and low risk groups.
A, B. SNP alterations were identified in high and low risk groups. C, D. The KM curve was used to verify the prognosis of different groups.
In this study, the immune status analysis in different groups found that the expression of most immune checkpoints was higher in the high-risk group. The heatmap shows that there were significant differences in risk scores for type I- or II-interferon (IFN) responses, MHC, inflammation promotion, checkpoints, T cell co-stimulation, etc. (Fig. 7A). In further immune infiltration analysis, the results showed that the high risk score was significantly negatively correlated with CD8+ T cells (r = −0.13, p = 0.035), CD4+ T cells (r = −0.16, p = 0.010), macrophages (r = −0.17, p = 0.004), NK cells (r = −0.18, p = 0.004) and B cells (r = −0.22, p = 0.000). These data also explain the poorer outcomes (Fig. 7B-H).
The immune microenvironment in different risk groups.
A. Heat map showing the enrichment of immune function in the gene set of ovarian serous cystadenocarcinoma (OV) samples with different risk scores. B. Correlation plot between the model and CD8+ T cell score. C. Correlation plot between the model and endothelial cell score. D. Correlation plot between the model and macrophage cell score. E. Correlation plot between the model and NK cell score. F. Correlation plot between the model and each uncharacterized cell score. G. Correlation plot between the model and B cell score. H. Correlation plot between the model and CD4+ T cell score.
Finally, this study first analyzed differentially expressed genes (DEGs) between high-risk and low-risk groups (Supplementary Fig. S2B). DEGs were then subjected to GO and KEGG enrichment to investigate their possible functional profiles. The top 3 GO terms for biological processes (BP) were ossification; sensory organ morphogenesis; ureteric bud development. The top 3 GO terms for cellular components (CC) were collagen-containing extracellular matrix; endoplasmic reticulum lumen; basement membrane. And the top 3 GO terms for molecular functions (MF) were extracellular matrix structural; receptor ligand activity; signaling receptor activator activity (Fig. 8A). The top 3 KEGG signaling pathways were ECM-receptor interaction; protein digestion and absorption; phagosome (Fig. 8B).
Functional enrichment analysis in different risk groups.
A. GO functional enrichment of differential genes. B. The KEGG pathway enrichment of differential genes.
The occurrence and development of OV is a complex biological process involving multiple factors. Exploring the role of a certain gene or signaling pathway in the occurrence and development of OV is only part of the overall molecular mechanism of OV (Armstrong et al. 2021; Coughlan and Testa 2021). Therefore, it is very important to explore the role of more molecules in the molecular mechanism of OV to clarify its molecular mechanism. Bioinformatics analysis combined with multi-gene prediction model construction is conducive to the development of customized personalized treatment plans and targeted drugs.
At present, changes in the expression levels of lncRNAs had been found in cancers such as breast cancer, prostate cancer, cervical cancer or ovarian cancer. Such changes in expression levels may serve as markers for cancer diagnosis and potential drug targets. The abnormal expression of lncRNA was related to the proliferation, invasion and metastasis of ovarian cancer cells. Among them, the up-regulated lncRNAs include PVT1, SOX2OT, BCYRN1, LSINCT5, etc., and the down-regulated lncRNAs include GAS5, SPRY4-IT1, linc00261, etc. (Wang et al. 2018; Wang et al. 2019; Tan et al. 2021). They enhance the ability of tumor cells to metastasize and become more aggressive by regulating processes such as EMT and immune escape (Li et al. 2016; Peng et al. 2017; Xing et al. 2021).
Current studies had confirmed that cuproptosis plays an important role in the regulation of tumor progression. However, few studies had paid attention to the relationship between them (Chen 2022; Wang et al. 2022). This study first used the GEPIA database to explore the expression of these 19 cuproptosis-related genes in OV tissues. A total of 113 lncRNAs were found to further identify cuproptosis-related lncRNAs using the above criteria. Finally, differentially expressed prognosis-related lncRNAs were selected by univariate Cox analysis and included in LASSO-Cox regression analysis. Finally, we obtained 5 lncRNAs (AC025287.2, AC092718.4, AC112721.2, LINC00996, and LINC01639) and constructed Cox proportional hazards model for them. Studies have shown that LINC01639 plays an important role in non-alcoholic steatohepatitis-induced liver cancer (Ilieva et al. 2022). In addition, other studies had shown that LINC00996 is involved in the regulation of multiple tumor processes including colorectal cancer, head and neck tumors, etc. (Ge et al. 2018; Zhou et al. 2020; Lina 2021; Yan et al. 2021; Cai et al. 2022). More importantly, AC092718.4 has been confirmed in an article on OV prognostic identification (Lin et al. 2021; Chen et al. 2022). However, our study also reported for the first time the role of AC025287.2 in the identification of OV prognosis. These data suggest that in-depth exploration of the potential regulatory mechanisms of these molecules may contribute to the targeted therapy and prognosis of OV.
Tumor microenvironment (TME) refers to the close relationship between the occurrence, growth and metastasis of tumor and the internal and external environment of tumor cells, which not only includes the structure, function and metabolism of tumor tissues, but also is related to the internal environment of tumor cells themselves (nuclear and cytoplasm) (Yang et al. 2020). As an important environment for regulating tumor development, invasion and metastasis, it is mainly composed of infiltrating immune cells, chemokines and cytokines (Yu et al. 2021). Under the combined stimulation of TME and tumor driver genes, tumor cells can reprogram normal cells of the host to achieve a more suitable phenotype and function for tumor growth. In TME, CD8+ T cells infiltrate the tumor, namely CD8+ tumor infiltrating lymphocytes (TILs), which are the symbol of immune recognition (Nowak and Klink 2020; Olalekan et al. 2021).
CD8+ TILs suggested that epithelial ovarian cancer patients had better prognosis. Regulation of CD8+ T cells is critical in ovarian TME by secreting granulozyme B, tumor necrosis factor, and interferon-γ (IFN-γ) to clear tumor cells. In the cancer environment where antigen exposure persists, T cells gradually lose their effector function, up-regulate the expression of inhibitory receptors, including PD-1, and exhibit poor memory (Bi et al. 2021; Guo et al. 2021). Such depletion of T cells can be reactivated by the inhibition receptor of PD-1 antibody. In addition, bone marrow-derived suppressor cells (MDSCs) have a significant ability to inhibit T cell responses, which are increased in ovarian cancer patients and play an important role in disease progression. Tumor-associated macrophages (TAMs) are considered to be the most abundant invasive immune cells in epithelial ovarian cancer (EOC) tissue and ascites (Nowak and Klink 2020). They possess an immunosuppressive M2 phenotype, characterized by expressions of CD163, CD204, CD206, and IL-10, and are associated with tumor progression. Natural killer cells (NK) exhibit excellent antitumor activity, including IFN-γ production and direct cytotoxicity to cancer cells. Fully mature human NK cell subsets were found to overexpress PD-1 in the ascites of ovarian cancer patients (Jiang et al. 2020; Maniati et al. 2020).
We further analyzed the relationship between the model and TMB, immunity and so on. TMB is an emerging cancer biomarker characterized by microsatellite instability, considered a predictor of response to tumor immunotherapy, and expressed differently in different cancers (Liu et al. 2019; Cheng et al. 2021; Ju et al. 2021). This study looked at mutations in both groups and found that for most genes, the high-risk group had a higher frequency of mutations than the low-risk group. TP53 was tumor suppressor genes, which is reported to inhibit the progress of cancer cells by regulating the cell cycle. Titin (TTN) was a basalin gene encoding sarcomere, which may be associated with dilated cardiomyopathy, and its role in tumorigenesis needs to be further studied (Yang et al. 2022a; Yang et al. 2022b). Current studies have shown that CD4+ T cells and CD8+ T cells are involved in killing tumor cells to inhibit tumor growth. Patients with higher infiltration levels of immune cells were more likely to have better immunotherapy effects and prognosis (Dumauthioz et al. 2018; Fujii and Shimizu 2019; Efimova et al. 2020; Deets and Vance 2021). Our results showed a trend of negative correlation between the scores of high-risk groups and the degree of infiltration of CD4+ and CD8+ T cells. It could explain the poorer prognosis of high risk scores. Finally, we try to deeply explore the potential functional mechanisms in different risk groups. We first screened the DEGs between these groups, and then performed enrichment analysis on the DEGs by GO and KEGG. Our results for GO enrichment found that the BP were ossification; sensory organ morphogenesis; ureteral bud development. The CC were collagen-containing extracellular matrix; basement membrane; endoplasmic reticulum lumen. The MF were extracellular matrix structure; receptor ligand activity; signaling receptor activator activity. At present, studies have shown that functions including extracellular matrix structure; receptor ligand activity and signaling receptor activator activity were significantly related to the OC process (Haslene-Hox et al. 2015; Zhao et al. 2018; Huang et al. 2022). This suggests that AC025287.2, AC112721.2, AC092718.4, LINC00996 and LINC01639 in the model of this study may function through this mechanism, which will be another evidence for the later mechanism exploration. Further KEGG results showed that DEGs significantly related pathways including extracellular matrix structure; receptor ligand activity and signaling receptor activator activity pathways. Studies had reported that the ECM-receptor interaction pathway was associated with a variety of tumor metastasis (Gao et al. 2018; Bao et al. 2019; Han et al. 2020; Nersisyan et al. 2021). Therefore, this pathway is still the main direction of our later discussion. Our data have confirmed the predictive efficiency and reliability of the model for OV. However, we still lack the necessary clinical cohorts to validate the potential of this model for clinical application. At present, our research group has started to construct a relevant cohort model, and we will further explore the mechanism of these 5 lncRNAs in future research to better understand their mechanism of action. In conclusion, this study successfully constructed a reliable cuproptosis-related lncRNA model for the prognosis of OV. These data provide a new insight into predicting the survival rate of OV patients and the effectiveness of clinical treatment.
This work was supported by Natural Science Foundation of Sichuan Province（Young Scholars Project, 2022NSFSC1532), National Natural Science Foundation of China for Young Scholars (No.82003865), Application research of the Science and Technology Department in Sichuan province, 2021YJ0461 (21YYJC2810) and the Fundamental Research Funds for the Central Universities（SCU2022F4080).
The authors declare no conflict of interest.
