[go: up one dir, main page]
Skip to main content
Download PDF

Development and validation of the potential biomarkers based on m6A-related lncRNAs for the predictions of overall survival in the lung adenocarcinoma and differential analysis with cuproptosis

BMC Bioinformatics volume 23, Article number: 327 (2022) Cite this article

Abstract

Background

The treatment and prognosis of lung adenocarcinoma (LUAD) remains a challenge. The study aimed to conduct a systematic analysis of the predictive capacity of N6-methyladenosine (m6A)-related long non-coding RNAs (lncRNAs) in the prognosis of LUAD.

Methods

594 samples were totally selected from a dataset from The Cancer Genome Atlas. The identification of prognostic m6A-related lncRNAs were performed by Pearson correlation analysis and Cox regression analysis. Systematic analyses, including cluster analysis, survival analysis, and immuno-correlated analysis, were conducted. A prognosis model was built from the optimized subset of m6A-related lncRNAs. The assessment of model was performed by survival analysis, and receiver operating characteristic (ROC) curve. Finally, the risk score of patients with LUAD calculated by the prognosis model was implemented by the analysis of Cox regression. Differential analysis was for further evaluation of the cuproptosis-related genes in two risk sets.

Results

These patients were grouped into two clusters according to the expression levels of 22 prognostic m6A-related lncRNAs. The patients with LUAD in cluster 2 was significantly worse in the overall survival (OS) (P = 0.006). Three scores calculated by the ESTIMATE methods in cluster 2 were significantly lower. After the least absolute shrinkage and selection operator algorithm, 10 prognostic m6A-related lncRNAs were totally selected to construct the final model to obtain the risk score. Then the area under the ROC curve of the prognosis model for 1, 3, and 5-year OS was 0.767, 0.709, and 0.736 in the training set, and 0.707, 0.691, and 0.675 in the test set. The OS of the low-risk cohort was significantly higher than that of the high-risk cohort in both the training set (P < 0.001) and test set (P < 0.001). After the analysis of Cox regression, the risk score [Hazard ratio (HR) = 5.792; P < 0.001] and stage (HR = 1.576; P < 0.001) were both considered as independent indicators of prognosis for LUAD. The expression levels of five cuproptosis-related genes were significantly different in two risk sets.

Conclusions

The study constructed a predictive model for the OS of patients with LUAD and these OS-related m6A-lncRNAs might have potential roles in LUAD progression.

Peer Review reports

Background

Lung cancer, a malignancy which has been identified to have the highest mortality and the second highest morbidity worldwidely, has a poor 5-year survival since the rate reaches only 10–20% [1, 2]. The most common subtype is lung adenocarcinoma (LUAD) [3, 4]. Despite therapeutic options including targeted therapy, radiotherapy, chemotherapy and surgery progressed rapidly, the prognosis of LUAD patients remains unsatisfactory [3,4,5]. LUAD, the pathogenesis of which involves multiple molecular mechanisms through various pathways, requires more in-depth mechanistic research to develop more promising therapies.

N6-methyladenosine (m6A) modification is a crucial RNA post-transcriptional modification in most eukaryotic long non-coding RNAs (lncRNAs) and mRNAs [6, 7]. It can modulate the stabilization, splicing, degradation, translation, and processing of the target RNA via three type of regulators, which include writers (methyltransferases), erasers (demethylases), and readers (binding proteins) [8, 9]. Some studies have indicated that m6A modifications are correlated with the oncogenesis and development of malignancies, including LUAD [10,11,12]. Wang previously showed that the 13 m6A regulators, such as KIAA1429 and FTO, were aberrantly overexpressed in tumor samples [11]. Moreover, several studies also found that the modification of m6A has involvement in the regulation of the immune response and tumor microenvironment infiltrating cells, such as dendritic cells [13, 14]. Therefore, it is necessary to study these m6A modifications to gain a comprehensive understanding their functions.

LncRNAs, a group of RNAs which is non-coded and has over 200 nucleotides in length, functions in various biological processes including growth, apoptosis, and regulation of cell development [15]. Aberrant regulation of lncRNAs has been investigated and found to be related to different malignancies, including in LUAD [16, 17]. Ding et al. found that high expression of lncRNA OGFRP1 leads to significantly worse survival outcomes and constructed a prediction model whose area under the receiver operating characteristic (ROC) curve (AUC) reaches 0.766 [18]. Studies focusing on the mechanism by which the m6A modification acts on the occurrence and progression of lncRNA-dependent LUAD are limited and the whole roles of the three types of m6A regulators in the aberrant lncRNAs regulation is still undefined [19]. Separately, the cell death was also involved in the modification of m6A for lncRNAs in the development of LUAD and cuproptosis was new form of cell death, which was found to be associated with the progression of cancer [20,21,22,23]. However, the relationship between cuproptosis and m6A-related lncRNAs in the development of LUAD has not been reported.

Therefore, the objective of our research was to perform a bioinformatic analysis to distinguish those survival-related m6A-lncRNAs based on data of LUAD patients from The Cancer Genome Atlas (TCGA) and conduct a systematic analysis. Furthermore, a model was constructed to predict the overall survival (OS) of these patients by utilizing the optimized set of survival-related m6A-lncRNAs. Differential analysis was for further evaluation of the relationship between the performance of the model and cuproptosis.

Methods

Datasets and m6A-related genes

The transcriptome profiling normalization data by fragments per kilobase of transcript per million mapped reads and the relevant clinical features were acquired from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov). A total of 515 cases with 594 samples were involved in this study from a TCGA project, and 59 samples were in the normal group and 535 in the tumor group. In addition, 23 m6A regulators were totally determined on the basis of previous studies, as shown the Additional file 1, which could be divided in to writers, erasers, and readers. Annotation of the lncRNAs and mRNAs in the TCGA project were conducted using the annotation file of Genome Reference Consortium Human Build 38 (GRCh38) acquired from the GENCODE website (https://www.gencodegenes.org). A total of 14,086 lncRNAs were identified by the Ensemble IDs of the genes in the TCGA dataset.

Selection of m6A-related lncRNAs and cluster analysis

m6A-related lncRNAs were initially identified by Pearson correlation analysis with the P < 0.001 and |Pearson R|> 0.5. Then, the prognostic m6A-related lncRNAs were determined by the univariate analysis of Cox regression. Wilcoxon tests were used for analyzing the prognostic m6A-related lncRNA expression between the normal group and the tumor group. Cluster analysis was used for analyzing the expression of the prognostic m6A-related lncRNAs. The OS of different clusters were compared by Kaplan–Meier (KM) curves and the log-rank test. Then, the prognostic m6A-related lncRNAs expression and clinical features in different clusters were analyzed by differential analysis.

Differential analysis and correlation analysis of CD274 expression

The CD274 (programmed cell death 1 ligand 1, PD-L1) expression in the normal group and the tumor group, and in the two clusters, were compared using the Wilcoxon test. Moreover, CD274 expression and the expression of prognostic m6A-related lncRNAs were analyzed by correlation analysis.

Tumor-infiltrating immune cells (TIIC) evaluation and gene set enrichment analysis (GSEA)

According to the CIBERSORT algorithm, standard gene expression data and the Wilcoxon test, the relative proportions and the difference of 22 TIIC subpopulations were evaluated [24]. The immune score, estimate score, and stromal score were computed using the ESTIMATE algorithm to further predict tumor purity and to analyze the tumor microenvironment [25]. The Wilcoxon test was used to analyze these three indicators of the tumor microenvironment between the two clusters. Then we used GSEA software (version 4.1.0) to investigate gene set enrichment in the two clusters of lung adenocarcinoma.

Establishment of the prognosis model and survival analysis

The patients with LUAD were grouped into a training cohort and a test cohort at a ratio of 5:5 randomly. The optimized subset of prognostic m6A-related lncRNAs was chosen using the least absolute shrinkage and selection operator (LASSO) algorithm to build the final prognosis model. After the number of lncRNAs was determined, the corresponding coefficients of each lncRNA of the optimized subset were calculated. The assessment of prognosis model was then performed by ROC analysis and a risk plot. The natural logarithm of risk score was calculated by summing the chosen lncRNAs expression, weighted by their corresponding coefficients. Then patients with LUAD were divided into low-risk set and high-risk set by the median value of the risk scores. The OS of patients with LUAD in the two risk sets were compared by the KM method and log-rank test in the training cohort and test cohort. The evaluation of the risk score and clinical features for the OS of patients with LUAD were implemented by univariate and multivariate analysis of Cox regression. Then, subgroup analysis of each clinical feature, based on the two risk sets, was conducted using the KM method and log-rank test.

Differential analysis of risk sets and risk score

The differential analysis of immune score, clinical features and cluster status between the two risk sets were analyzed by chi-squared test. The differential analysis of PD-L1 expression and the cuproptosis-related genes expression were conducted using the Wilcoxon test between the two risk sets. Correlation analysis was also implemented between the risk score and each TIIC of samples using Spearman correlation analysis.

Statistical analysis

The R statistical software (version 4.0.5) was used for all statistical analyses. The “limma” package was used to assess the transcriptome profiling data. The survival analysis was conducted by the “survival” and, “survminer” packages. The “pheatmap”, “reshape2”, “ggpubr”, and “ggplot2” packages were used to generate heatmaps, boxmaps and risk plots. The “ConsensusClusterPlus” package was to conduct the cluster analysis. The “corrplot” package was to conduct the correlation analysis of PD-L1 expression. The “e1071”, “parallel”, and “preprocessCore” packages were used to evaluate the TIICs. The “estimate” package was to obtain the stromal score, immune score, and estimate score of sample. The “vioplot”, “ggpubr” “ggplot2”, and “ggExtra” packages were used to investigate the differences among the TIICs in the two clusters and to conduct correlation analysis between the risk score and the TIICs. The “caret”, “glmnet”, “timeROC”, “survival”, and “survminer” packages were to build the model and to get the ROC. A two-side P < 0.05 was considered statistically significant.

Results

Determination of m6A-related lncRNAs in patients with LUAD

1558 m6A-related lncRNAs were totally determined using Pearson correlation analysis (with a P < 0.001 and |Pearson R|> 0.5). After clinical survival data integrated, 504 patients were further involved in the prognostic analysis from a TCGA project. The information for 504 patients was shown in Table 1. Then, a total of 22 prognostic m6A-related lncRNAs were identified by univariate Cox regression analysis (with P < 0.01, Fig. 1a). These prognostic m6A-related lncRNAs had significantly different expression levels between normal and tumor samples by the Wilcoxon test (Fig. 1b). In tumor samples, the expression levels of certain prognostic m6A-related lncRNAs, such as AC099850.4, AL606489.1, AC010999.2, and AC034102.8, were higher than in normal samples; however, the expression levels of other prognostic m6A-related lncRNAs, such as PAN3-AS1, AF131215.5, AC024075.1, MIR99AHG, AC005884.1, and AC090617.5, were lower than in normal samples (Fig. 1b).

Table 1 Clinical features of 504 patients with LUAD
Full size table
Fig. 1
figure 1

The selection of m6A-related lncRNAs for prognosis and cluster analysis. a Forest plot of prognostic m6A-related lncRNAs of patients with LUAD using univariate Cox regression analysis (P < 0.01, CI: confidence interval); b Boxplot showing the expression levels of prognostic m6A-related lncRNAs in the normal and tumor groups. c Heatmap showing the relationship between clinical features, prognostic m6A-related lncRNA expression, and the two clusters. d Kaplan–Meier curves showing that cluster 2 group had worse overall survival than cluster 1 by log-rank test (P = 0.006). lncRNA, long non-coding RNA; LUAD, lung adenocarcinoma. “***”: P < 0.001, “**”: P < 0.01, “*”: P < 0.05

Full size image

Cluster analysis of m6A-related lncRNAs

According to the expression of the 22 prognostic m6A-related lncRNAs, all samples were divided into two clusters by cluster analysis. Then, the differences in clinical features and the 22 prognostic m6A-related lncRNAs between the two clusters were investigated. A heatmap showed that the stage classification was significantly different between the two clusters (Fig. 1c). Moreover, the OS between two clusters were compared. As shown in Fig. 1d, the OS of patients with LUAD was significantly different between the two clusters (P = 0.006).

Immuno-correlated analysis

The PD-L1 expression in the tumor group was significantly lower than that in the normal group (Fig. 2a). The PD-L1 expression of cluster 1 was significantly lower than that in cluster 2 (Fig. 2b). Figure 2c shows the correlation between PD-L1 expression and the expression levels of the 22 prognostic m6A-related lncRNAs. The infiltration of 22 immune cells in the LUAD microenvironment were analyzed based on the CIBERSORT algorithm. As shown in Fig. 2d, the 22 kinds of TIICs were compared between the two clusters. The violin plot showed that the proportions of 13 kinds of immune cells were significantly different between the two clusters. Then the stromal score, immune score, and estimate score assessed by the ESTIMATE methods were compared between two clusters and all three were significantly different (Fig. 2e–g).

Fig. 2
figure 2

Immuno-correlated analysis. a Boxplot showing that the CD274 (PD-L1) expression was significantly different between the normal set and the tumor set. ***: P < 0.001. b Boxplot showing that PD-L1 expression was significantly different between the two clusters. ***: P < 0.001. c Plot showing the correlation between CD274 (PD-L1) expression and the expression of 22 prognostic m6A-related lncRNA. The red circle indicates a positive correlation and the blue circle indicates a negative correlation. “”: P < 0.05. d Violin plot showing the differences in the percentage of 22 kinds of tumor-infiltrating immune cells between the two clusters. e, f, g: These boxplots show that the estimate score, immune score, and stromal score of cluster 2 were significantly lower than cluster 1 (P = 0.0058, 0.019, and 0.005, respectively). PD-L1, programmed cell death 1 ligand 1; lncRNA, long non-coding RNA

Full size image

GSEA analysis

GSEA was conducted and the results showed that in cluster 2, the genes set was mainly enriched for the oocyte meiosis, cell cycle, and ubiquitin mediated proteolysis (Fig. 3a–c). In cluster 1, the genes set was mainly enriched for arachidonic acid metabolism, linoleic acid metabolism, and asthma (Fig. 3d–f).

Fig. 3
figure 3

Gene set enrichment analysis of the two clusters. a–c: The top three gene sets enriched in cluster 2 included oocyte meiosis (NES = 2.65, Norm P < 0.001, FDR q < 0.001), cell cycle (NES = 2.60, Norm P < 0.001, FDR q < 0.001), and ubiquitin mediated proteolysis (NES = 2.52, Norm P < 0.001, FDR q < 0.001). d–f: The top three gene sets enriched in cluster 1 included arachidonic acid metabolism (NES =  − 2.11, Norm P < 0.001, FDR q = 0.014), linoleic acid metabolism (NES =  − 2.08, Norm P < 0.001, FDR q = 0.013) and Asthma (NES =  − 2.07, Norm P = 0.002, FDR q = 0.009). NES, normalized enrichment score; Norm, normalized; FDR, false discovery rate

Full size image

Construction of the prognosis model and survival analysis

The 504 patients were randomly grouped into a training cohort and a validation cohort. Then 10 prognostic m6A-related lncRNAs chosen by the LASSO algorithm were to obtain the risk score, including AC087501.4, L3MBTL2-AS1, AL606489.1, AC007613.1, AC090617.5, AC073316.3, AC010999.2, AC005884.1, TSPOAP1-AS1, and ADPGK-AS1(Fig. 4a). The formula is shown in the Additional file 2. Built by the LASSO algorithm, the AUC of the prognosis model for 1-year, 3-year, 5-year OS was 0.767, 0.709, and 0.736 in the training cohort and 0.707, 0.691, and 0.675 in the validation cohort (Fig. 4b–c). Then, according to the median value of the risk scores of the training cohort, the survival status of the training cohort and the validation cohort were analyzed based the two risk sets (Fig. 4d–g). Moreover, the OS of low-risk set were significantly higher than the high-risk set in both the training cohort (P < 0.001) and validation cohort (P < 0.001) (Fig. 4h–i).

Fig. 4
figure 4

Building the prognosis model and survival analysis. a The process of least absolute shrinkage and selection operator regression to establish the final prognosis model. b, c The receiver operating characteristic of the prognosis model for OS in the training cohort (b) and validation cohort (c). d, e Patients were grouped into two risk sets by mean value of the risk score of training cohort in the training cohort (d) and validation cohort (e). f, g The distribution of patients with survival status based on two risk groups in the training cohort (f) and validation cohort (g). h, i: Kaplan–Meier curves of two risk sets in the training cohort (h) and validation cohort (i)

Full size image

Cox regression analysis and subgroup analysis

After univariate and multivariate analysis of Cox regression, the risk score and stage were both considered as independent predictive factors for the OS of patients with LUAD and the risk score was the key one in all groups (Table 2). After the subgroup analysis, the low-risk set showed significantly higher OS than the high-risk set in these 11 subgroup patients, as shown in Fig. 5. Then, T stage, TNM stage, N stage, sex, immune score, and cluster of patients were significantly different in the two risk groups (Fig. 6). Figure 6 also showed the difference expression levels of the 10 prognostic m6A-related lncRNAs in the two risk groups. The PD-L1 expression levels of the low-risk set and high-risk set were not significantly different (P = 0.54), as shown in the Fig. 7. The expression levels of five cuproptosis-related genes were significantly different in two risk sets, as shown in the Fig. 8. Moreover, the 6 kinds of immune cells, such as plasma cells (R =  − 0.2), correlated significantly and negatively with the risk score, whereas the other 6 kinds of immune cells, such as M0 macrophages (R = 0.22), correlated significantly and positively with the risk score (Fig. 9).

Table 2 Univariate and multivariate Cox regression analysis
Full size table
Fig. 5
figure 5

Survival analysis of clinical subgroups of patients. a–k: Kaplan–Meier curves and the log-rank test showed that the overall survival of the high-risk set was significantly worse than that of the low-risk set in these 11 subgroups of patients (a–k)

Full size image
Fig. 6
figure 6

Heatmap analysis of the two risk sets in all patients.The figure shows that the differences of patient with N stage, T stage, Stage, Sex, Immune Score, and Cluster was respectively significant in two risk sets. The heatmap also shows the expression of the 10 prognostic m6A-related lncRNAs in the two risk sets.lncRNA, long non-coding RNA.“***”: P < 0.001, “**”: P < 0.01, “*”: P < 0.05

Full size image
Fig. 7
figure 7

The difference of CD274 (PD-L1) expression in the two risk sets. The boxplot shows that the PD-L1 expression of the high-risk set and low-risk set were no significant difference (P = 0.54). PD-L1, programmed cell death 1 ligand 1

Full size image
Fig. 8
figure 8

Differential analysis of cuproptosis-related genes in the two risk sets.The boxplot shows that the difference of the expression of five cuproptosis-related genes were significant in the two risk sets. Log-transformation was only performed to draw the boxplot for visualization.“ns”: not significant; “***”: P < 0.001; “**”: P < 0.01; “*”: P < 0.05

Full size image
Fig. 9
figure 9

The correlation between immune cells of samples and the risk score.The figure shows 12 types of immune cells from samples with significantly related with risk scores

Full size image

Discussion

The transcriptome information and clinical data of patients that have been diagnosed with LUAD were totally collected from the TCGA dataset, with the aim to investigate those m6A-related lncRNAs that associated with LUAD prognosis. After the analysis of univariate Cox regression, 22 prognostic m6A-related lncRNAs were confirmed and used for clusters analysis. The OS of the patients in cluster 1 was significantly higher than in cluster 2. Then 10 m6A-related lncRNAs were selected from the previously mentioned 22 prognostic biomarkers, based on which a prediction model aiming at the OS of patients was established. The low-risk set patients with LUAD had significantly higher OS than the other set (P < 0.001). The risk score was considered as an independent risk factor for OS after the univariate and multivariate analyses of Cox regression. In the later subgroup analysis, the low-risk set had a significantly higher OS revealed by groups with M0 stage and divided by sex, age, TNM stage, N stage, and T stage. Moreover, the difference of T stages, TNM stages, N stages, sex, immune scores, five cuproptosis-related genes and clusters of these patients were significant in two risk sets.

Researches have shown that m6A modification of lncRNAs could regulate the oncogenesis and progression of malignancies [26,27,28]. Lan et al. showed that VIRMA, which was also named KIAA1429, could regulate lncRNA GATA3 to drive hepatocarcinogenesis, progression, and metastasis based on an m6A modification [26]. Yang et al. showed that METTL14, which is one of m6A “writer” regulators, could increase the m6A-methylated rate of lncRNA XIST and decrease the levels of XIST to suppress the proliferation and metastasis of colorectal cancer cells [27]. Ni et al. showed that YTHDF3 could promote m6A-related lncRNA GA55 degradation and then accelerate the development of colorectal cancer [28]. These studies indicated that m6A-related lncRNAs could contribute to not only the oncogenesis but also progression of malignant tumors and lncRNAs might be involved in competing against endogenous RNAs, affecting the invasiveness of tumors. However, taking a thorough look at how m6A modifications act on the progression of lncRNA-dependent LUAD is challengeable. Our present study found the gene set of m6A-related lncRNAs in LUAD mainly enriched in the oocyte meiosis, cell cycle, and ubiquitin mediated proteolysis in cluster 2, whose OS was obviously poorer than in cluster 1. These results indicated that lncRNAs with the m6A modification could influence the progression of LUAD and then the OS of patients with LUAD. Therefore, we believe these lncRNAs will be confirmed as potential therapeutic targets of LUAD by studying the functions and interaction of these lncRNAs and their m6A modifications in a future study.

In the present study, 10 of 22 m6A-related lncRNAs were selected by the LASSO algorithm, including AC087501.4, L3MBTL2-AS1, AL606489.1, AC007613.1, AC090617.5, AC073316.3, AC010999.2, AC005884.1, TSPOAP1-AS1, and ADPGK-AS1, which could predict the prognosis of LUAD. Several of the 10 lncRNAs have been investigated in different tumors [29,30,31]. Giulietti et al. showed that the methylation level of the TSPOAP1-AS1 promoter was higher in pancreatic ductal adenocarcinoma than in normal tissues and the aberrant methylation level of the lncRNA might be considered as an indicator for the diagnosis of pancreatic ductal adenocarcinoma [29]. Yang et al. showed that the expression level ADPGK-AS1 in the adjacent non-tumor tissues of patients with breast cancer was lower than in the tumor tissues, which might be regarded as an indicator for the prognosis of breast cancer patients because ADPGK-AS1 could facilitate cell proliferation and migration, suppress cell apoptosis, and induce epithelial-mesenchymal transition [30]. The study of Xu et al. showed that the risk score comprised of 12 m6A-related lncRNAs could predict the 1-year OS of patients with LUAD, and the AUC of its risk score was 0.759 [19]. The group of 12 m6A-related lncRNAs was different with our study and there might be some reasons as below. Firstly, our study included 23 m6A-related genes but the study of Xu et al. only included 21 m6A-related genes [19]. Therefore, our study performed further analysis with two more m6A-related genes than the study of Xu et al. [19]. Secondly, the Pearson's correlation coefficient of their study and our study was more than 0.3 and 0.5, respectively [19]. Thirdly, the prognosis of m6A-related lncRNAs selected by univariate Cox regression analysis in their study and our study was different (P < 0.05 and P < 0.01, respectively) [19]. Fourthly, sample size was insufficient. Moreover, the version of the annotation files might be different. What’s more, the study of Zhao et al. used 13 m6A-related lncRNAs to construct the prediction model of patients with LUAD [32]. The two studies were similar to some extent, but there were some differences between the results of our study and the study of Zhao et al. [32]. The main reasons are as follows. Firstly, 3 of 23 m6A-related genes are different in our study and the study of Zhao et al. [32]. IGF2BP1, IGF2BP2, IGF2BP3 were also included in the study of Xu et al., which are consistent with our study but not with the study of Zhao et al. [19, 32]. Therefore, there are different in m6A-related genes between our study and the study of Zhao et al. [32]. Secondly, there are different in m6A-related lncRNAs between our study and the study of Zhao et al. because of the above mentioned the differences of m6A-related genes [32]. Thirdly, a total of 504 patients were included into our study to construct the model, but 468 patients were included into the study of Zhao et al. [32]. Although two studies were derived from the TCGA lung adenocarcinoma dataset, there were differences in the transcriptome data due to the inconsistency of the patients included. Fourthly, 22 m6A- related lncRNAs were obtained by the univariate Cox regression analysis (with P < 0.01) in our study, while 91 m6A- related lncRNAs were obtained by the univariate Cox regression analysis (with P < 0.05) in the study of Zhao et al. [32]. Although there have been few studies about the relationship of how m6A-related lncRNAs act on the occurrence and progression in LUAD, based on our results, we still believe these m6A-lncRNAs have roles in LUAD tumorigenesis and progression, which need to be confirmed in future studies.

Several studies have investigated the relationship between TIICs and the prognosis of lung cancer [33,34,35]. Li et al. found that in the low-risk set, the percentage of neutrophil infiltration was lower while patients in the high-risk set had poorer prognosis in nonsquamous non-small cell lung cancer [33]. Sun et al. showed that LUAD groups with lower immune score, stromal score, or estimate scores could have worse OS than those with higher scores [34]. The finding is in accordance with the results of this study. We observed no difference in PD-L1 expression between the two risk sets. Similarly, Zhang et al. showed that there was no significant difference in PD-L1 expression among patients with adenocarcinoma in situ, minimally invasive adenocarcinoma, and invasive adenocarcinoma and this finding is also in line with ours [35]. However, evidence suggest that immune checkpoint inhibitors targeting the PD-1/PD-L1 interaction could lead to a reversal of the lung cancer–induced immunosuppressive microenvironment, bringing about an effective host antitumor immune response and also significant improvement in survival [36,37,38]. Therefore, further research into the undying association between the CD274 (PD-L1) expression and prognosis of LUAD is required. Moreover, our study showed that the expression levels of five cuproptosis-related genes were significantly different in two risk sets. The result indicated cuproptosis might have the involvement in the modification of m6A for lncRNAs in the development of LUAD, which might help to identify tumor therapeutic targets, and the further investigation is also required.

Several limitations are worth mentioning. Firstly, the study only included a TCGA dataset and the prognostic m6A-related lncRNAs need to be validated using various datasets or independent LUAD cohorts in further research. Secondly, the detailed mechanism and function of the prognostic m6A-related lncRNAs in tumorigenesis and the progression of LUAD need to be confirmed by in vitro and in vivo experiments.

Conclusions

In conclusion, the present study investigated the association between m6A-related lncRNAs and the prognosis in patients with LUAD, establishing a prognosis model for the OS of these patients.

The difference in the expression levels of five cuproptosis-related genes were confirmed in two risk sets. These results suggested potential therapeutic targets of LUAD, which might be confirmed by studying the functions and mechanism of m6A-related lncRNAs in the future.

Availability of data and materials

All data, including the transcriptome profiling normalization data and the relevant clinical features, were acquired from the Genomic Data Commons Data Portal of TCGA (https://portal.gdc.cancer.gov/repository), which is the open access data. The data analysis was also included into the Additional file 3.

Abbreviations

AUC:

Area under the receiver operating characteristic curve

CI:

Confidence interval

FDR:

False discovery rate

GRCh38:

Genome reference consortium human build 38

GSEA:

Gene set enrichment analysis

HR:

Hazard ratio

LASSO:

Least absolute shrinkage and selection operator

lncRNAs:

Long non-coding RNAs

LUAD:

Lung adenocarcinoma

m6A:

N6-methyladenosine

NES:

Normalized enrichment score

Norm:

Normalized

OS:

Overall survival

PD-L1:

Programmed cell death 1 ligand 1

ROC:

Receiver operating characteristic

TCGA:

The Cancer Genome Atlas

TIIC:

Tumor-infiltrating immune cells

KM:

Kaplan–Meier

References

  1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021;71(1):7–33. https://doi.org/10.3322/caac.21654.

    Article  PubMed  Google Scholar 

  2. Allemani C, Matsuda T, Di Carlo V, et al. Global surveillance of trends in cancer survival 2000–14 (CONCORD-3): analysis of individual records for 37 513 025 patients diagnosed with one of 18 cancers from 322 population-based registries in 71 countries. Lancet. 2018;391(10125):1023–75. https://doi.org/10.1016/S0140-6736(17)33326-3.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Herbst RS, Morgensztern D, Boshoff C. The biology and management of non-small cell lung cancer. Nature. 2018;553(7689):446–54. https://doi.org/10.1038/nature25183.

    Article  CAS  PubMed  Google Scholar 

  4. Lee JJ, Park S, Park H, et al. Tracing oncogene rearrangements in the mutational history of lung adenocarcinoma. Cell. 2019;177(7):1842-57.e21. https://doi.org/10.1016/j.cell.2019.05.013.

    Article  CAS  PubMed  Google Scholar 

  5. Hirsch FR, Scagliotti GV, Mulshine JL, et al. Lung cancer: current therapies and new targeted treatments. Lancet. 2017;389(10066):299–311. https://doi.org/10.1016/S0140-6736(16)30958-8.

    Article  CAS  PubMed  Google Scholar 

  6. Yue Y, Liu J, He C. RNA N6-methyladenosine methylation in post-transcriptional gene expression regulation. Genes Dev. 2015;29(13):1343–55. https://doi.org/10.1101/gad.262766.115.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Wang X, Lu Z, Gomez A, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505(7481):117–20. https://doi.org/10.1038/nature12730.

    Article  CAS  PubMed  Google Scholar 

  8. He RZ, Jiang J, Luo DX. The functions of N6-methyladenosine modification in lncRNAs. Genes Dis. 2020;7(4):598–605. https://doi.org/10.1016/j.gendis.2020.03.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m6A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 2018;28(6):616–24. https://doi.org/10.1038/s41422-018-0040-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Cheng Y, Wang M, Zhou J, Dong H, Wang S, Xu H. The important role of N6-methyladenosine RNA modification in non-small cell lung cancer. Genes (Basel). 2021;12(3):440. https://doi.org/10.3390/genes12030440.

    Article  CAS  Google Scholar 

  11. Wang H, Zhao X, Lu Z. m6A RNA methylation regulators act as potential prognostic biomarkers in lung adenocarcinoma. Front Genet. 2021;12: 622233. https://doi.org/10.3389/fgene.2021.622233.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Li Y, Gu J, Xu F, et al. Molecular characterization, biological function, tumor microenvironment association and clinical significance of m6A regulators in lung adenocarcinoma. Brief Bioinform. 2020. https://doi.org/10.1093/bib/bbaa225.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Shulman Z, Stern-Ginossar N. The RNA modification N6-methyladenosine as a novel regulator of the immune system. Nat Immunol. 2020;21(5):501–12. https://doi.org/10.1038/s41590-020-0650-4.

    Article  CAS  PubMed  Google Scholar 

  14. Xu R, Pang G, Zhao Q, et al. The momentous role of N6-methyladenosine in lung cancer. J Cell Physiol. 2021;236(5):3244–56. https://doi.org/10.1002/jcp.30136.

    Article  CAS  PubMed  Google Scholar 

  15. Statello L, Guo CJ, Chen LL, Huarte M. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol. 2021;22(2):96–118. https://doi.org/10.1038/s41580-020-00315-9.

    Article  CAS  PubMed  Google Scholar 

  16. Li Y, Jiang T, Zhou W, et al. Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nat Commun. 2020;11(1):1000. https://doi.org/10.1038/s41467-020-14802-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Salavaty A, Rezvani Z, Najafi A. Survival analysis and functional annotation of long non-coding RNAs in lung adenocarcinoma. J Cell Mol Med. 2019;23(8):5600–17. https://doi.org/10.1111/jcmm.14458.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Ding Y, Liu JH. The signature lncRNAs associated with the lung adenocarcinoma patients prognosis. Math Biosci Eng. 2019;17(2):1593–603. https://doi.org/10.3934/mbe.2020083.

    Article  PubMed  Google Scholar 

  19. Xu F, Huang X, Li Y, Chen Y, Lin L. m6A-related lncRNAs are potential biomarkers for predicting prognoses and immune responses in patients with LUAD. Mol Ther Nucl Acids. 2021;24:780–91. https://doi.org/10.1016/j.omtn.2021.04.003.

    Article  CAS  Google Scholar 

  20. Liu L, Li H, Hu D, et al. Insights into N6-methyladenosine and programmed cell death in cancer. Mol Cancer. 2022;21(1):32. https://doi.org/10.1186/s12943-022-01508-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Tsvetkov P, Coy S, Petrova B, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science. 2022;375(6586):1254–61. https://doi.org/10.1126/science.abf0529.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Tang D, Chen X, Kroemer G. Cuproptosis: a copper-triggered modality of mitochondrial cell death. Cell Res. 2022. https://doi.org/10.1038/s41422-022-00653-7.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Ge EJ, Bush AI, Casini A, et al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat Rev Cancer. 2022;22(2):102–13. https://doi.org/10.1038/s41568-021-00417-2.

    Article  CAS  PubMed  Google Scholar 

  24. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7. https://doi.org/10.1038/nmeth.3337.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612. https://doi.org/10.1038/ncomms3612.

    Article  CAS  PubMed  Google Scholar 

  26. Lan T, Li H, Zhang D, et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer. 2019;18(1):186. https://doi.org/10.1186/s12943-019-1106-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yang X, Zhang S, He C, et al. METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST. Mol Cancer. 2020;19(1):46. https://doi.org/10.1186/s12943-020-1146-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Ni W, Yao S, Zhou Y, et al. Long noncoding RNA GAS5 inhibits progression of colorectal cancer by interacting with and triggering YAP phosphorylation and degradation and is negatively regulated by the m6A reader YTHDF3. Mol Cancer. 2019;18(1):143. https://doi.org/10.1186/s12943-019-1079-y.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Giulietti M, Righetti A, Principato G, Piva F. LncRNA co-expression network analysis reveals novel biomarkers for pancreatic cancer. Carcinogenesis. 2018;39(8):1016–25. https://doi.org/10.1093/carcin/bgy069.

    Article  CAS  PubMed  Google Scholar 

  30. Yang J, Wu W, Wu M, Ding J. Long noncoding RNA ADPGK-AS1 promotes cell proliferation, migration, and EMT process through regulating miR-3196/OTX1 axis in breast cancer. In Vitro Cell Dev Biol Anim. 2019;55(7):522–32. https://doi.org/10.1007/s11626-019-00372-1.

    Article  CAS  PubMed  Google Scholar 

  31. Tu Z, Wu L, Wang P, et al. N6-Methylandenosine-related lncRNAs are potential biomarkers for predicting the overall survival of lower-grade glioma patients. Front Cell Dev Biol. 2020;8:642. https://doi.org/10.3389/fcell.2020.00642.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Zhao J, Lin X, Zhuang J, He F. Relationships of N6-methyladenosine-related long non-coding RNAs with tumor immune microenvironment and clinical prognosis in lung adenocarcinoma. Front Genet. 2021;20(12): 714697. https://doi.org/10.3389/fgene.2021.714697 (PMID:34777460;PMCID:PMC8585518).

    Article  CAS  Google Scholar 

  33. Li B, Cui Y, Diehn M, Li R. Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol. 2017;3(11):1529–37. https://doi.org/10.1001/jamaoncol.2017.1609.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Sun S, Guo W, Wang Z, et al. Development and validation of an immune-related prognostic signature in lung adenocarcinoma. Cancer Med. 2020;9(16):5960–75. https://doi.org/10.1002/cam4.3240.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhang C, Zhang J, Xu FP, et al. Genomic landscape and immune microenvironment features of preinvasive and early invasive lung adenocarcinoma. J Thorac Oncol. 2019;14(11):1912–23. https://doi.org/10.1016/j.jtho.2019.07.031.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Borghaei H, Paz-Ares L, Horn L, et al. Nivolumab versus docetaxel in advanced nonsquamous non-small-cell lung cancer. N Engl J Med. 2015;373(17):1627–39. https://doi.org/10.1056/NEJMoa1507643.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Brahmer J, Reckamp KL, Baas P, et al. Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer. N Engl J Med. 2015;373(2):123–35. https://doi.org/10.1056/NEJMoa1504627.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Garon EB, Rizvi NA, Hui R, et al. Pembrolizumab for the treatment of non-small-cell lung cancer. N Engl J Med. 2015;372(21):2018–28. https://doi.org/10.1056/NEJMoa1501824.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank the contributions of the TCGA project.

Funding

This study was supported by the Medical and Health Science and Technology Project of Zhejiang Province (2022KY230, 2019KY117), National Natural Science Foundation of China (82102128), "Pioneer" and "Leading Goose" R&D Program of Zhejiang (2022C03046), Research Project of Zhejiang Chinese Medical University (2022JKJNTZ19, 2021JKGJYY007) and the Natural Science Foundation of Zhejiang Province (LSY19H180003).

Author information

Authors and Affiliations

  1. Department of Radiology, The First Affiliated Hospital of Zhejiang Chinese Medical University (Zhejiang Provincial Hospital of Traditional Chinese Medicine), 54 Youdian Road, Hangzhou, China

    Chen Gao, Ning Kong, Fan Zhang, Maosheng Xu & Linyu Wu

  2. The First School of Clinical Medicine, Zhejiang Chinese Medical University, Hangzhou, China

    Chen Gao, Ning Kong, Fan Zhang, Maosheng Xu & Linyu Wu

  3. Department of Surgery, The Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China

    Liuzhi Zhou

Authors
  1. Chen Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Ning Kong
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Fan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Liuzhi Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Maosheng Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Linyu Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

All authors contributed to this paper. CG: Investigation, Data Curation, Software, Writing—Original Draft. NK: Investigation, Visualization. FZ: Formal analysis, Visualization. LZ: Data Curation, Software, Validation. MX and LW: Conceptualization, Project Administration, Methodology, Supervision, Writing—Reviewing and Editing. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Maosheng Xu or Linyu Wu.

Ethics declarations

Ethics approval and consent to participate

TCGA belongs to public databases. The patients involved in the database have obtained ethical approval. Users can download relevant data for free for research and publish relevant articles. Our study is based on open access data, so there are no ethical issues and other conflicts of interest. All methods were carried out in accordance with relevant guidelines and regulations.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. The official full name of 23 m6A regulators.

Additional file 2

. The formula of risk score.

Additional file 3

. The data analysis.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gao, C., Kong, N., Zhang, F. et al. Development and validation of the potential biomarkers based on m6A-related lncRNAs for the predictions of overall survival in the lung adenocarcinoma and differential analysis with cuproptosis. BMC Bioinformatics 23, 327 (2022). https://doi.org/10.1186/s12859-022-04869-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12859-022-04869-7

Keywords

BMC Bioinformatics

ISSN: 1471-2105

Contact us