Exploration of methylation-driven genes for monitoring and prognosis of patients with lung adenocarcinoma.

BACKGROUND: As one of the most common malignant tumors in humans, lung cancer has experienced a gradual increase in morbidity and mortality. This study examined prognosis-related methylation-driven genes specific to lung adenocarcinoma (LUAD) to provide a basis for prognosis prediction and personalized targeted therapy for LUAD patients.
METHODS: The methylation and survival time data from LUAD patients in the TCGA database were downloaded. The MethylMix algorithm was used to identify the differential methylation status of LUAD and adjacent tissues based on the β-mixture model to obtain disease-related methylation-driven genes. A COX regression model was then used to screen for LUAD prognosis-related methylation-driven genes, and a linear risk model based on five methylation-driven gene expression profiles was constructed. A methylation and gene expression combined survival analysis was performed to further explore the prognostic value of 5 genes independently.
RESULTS: There were 118 differentially expressed methylation-driven genes in the LUAD tissues and adjacent tissues. Five of the genes, CCDC181, PLAU, S1PR1, ELF3, and KLHDC9, were used to construct a prognostic risk model. Overall, the survival time was significantly lower in the high-risk group compared with that in the low-risk group (P < 0.05). In addition, the methylation and gene expression combined survival analysis found that the combined expression levels of the genes CCDC181, PLAU, and S1PR1 as well as KLHDC9 alone can be used as independent prognostic markers or drug targets.
CONCLUSION: Our findings provide an important bioinformatic basis and relevant theoretical basis for guiding subsequent LUAD early diagnosis and prognosis assessments.

Background

Lung cancer is one of the most malignant tumors in the world, with high morbidity and mortality [1]. Epidemiological statistics show that the incidence of lung adenocarcinoma (LUAD) in lung cancer is increasing, and it is more common among women and non-smokers [2, 3]. Due to the existence of tumor heterogeneity factors and the different pathological and molecular types of patients, the personalized treatment of LUAD faces large challenges [4]. In recent years, with the advancement of molecular biology, several LUAD driving genes such as EGFR, KARS, and TP53 have been discovered one after the other [5], and some targeted therapeutic drugs have been developed, which has greatly improved the treatment of some patients with LUAD [6, 7]. However, not all patients have benefited from these drugs, such as those who cannot tolerate the therapeutic side effects of targeted drugs or those who are resistant to drugs in the short term. Therefore, identification of new LUAD-associated driver genes by bioinformatics analysis and construction of risk model are necessary for the prognosis evaluation and post-treatment of patients.

Molecular mechanism research based on bioinformatic analysis is an important method for cancer research. Not only can this research explore the molecular pathogenesis of tumors but also it can identify biomarkers for early diagnosis, treatment and prognosis of tumors [8]. As one of the core elements in epigenetic modification, DNA methylation is an important signaling tool that regulates genomic function and is one of the important features mediating carcinogenesis [9]. DNA methylation can participate in many cellular processes such as cell differentiation, genome stability, and gene imprinting [10, 11]. Biological processes, specifically, changes in DNA methylation, can provide an important basis for early diagnosis and prognosis of cancer and a new ideas for further clinical applications.

In recent years, sufficient research evidence has shown that the occurrence and development of lung cancer is a multi-factor, multi-stage and multi-gene change process in which the inactivation of tumor suppressor gene methylation is one of its important mechanisms [12-15]; consequently, in-depth research on the lung cancer-related methylation mechanism has become a focus of great concern. A large number of studies have shown that the methylation of certain genes can affect their expression, and this phenomenon is closely related to the diagnosis and prognosis of lung cancer [16]. For example, Feng et al. showed that RAR-β methylation levels were abnormally expressed in non-small cell lung cancer (NSCLC) patients, while those with positive APC methylation status in tumor tissues survived longer than those with negative APC methylation, indicating that the methylation of RAR-β or APC is a promising diagnostic or prognostic marker of NSCLC [17]. In addition, Zhang et al. found that the methylation level of PAX6 in non-small cell lung cancer tissues was higher than that in normal tissues and that the methylation status of PAX6 was associated with a poor overall survival rate in cancer tissues. It has been suggested that methylated PAX6 may be a useful biomarker for the prognosis assessment of NSCLC [18].

The Cancer Genome Atlas (TCGA) database [19] provides researchers around the world with open data on cancer genetic and epigenetic profiles, enabling researchers to efficiently capture and analyze relevant data and genomic changes. MethylMix, an algorithm implemented with R, identifies methylation status based on a β-mixed model to identify disease-specific hypomethylated and hypermethylated genes to obtain disease-related methylation-driven genes [20]. At present, related studies on methylation-driven genes have been reported [21]; in this study, the methylation and mRNA expression data of LUAD patients were extracted from the TCGA database, and the methylation-driven genes related to LUAD were obtained by use of the MethylMix algorithm. The Cox–Kaplan–Meier-survival method was used to construct a survival model, assess the methylation-driven genes related to LUAD prognosis, explore the correlation between DNA aberrant methylation and LUAD genome level, and provide a scientific basis for personalized medicine.

Materials and methods

Data processing and analysis

We downloaded methylation and mRNA expression data from LUAD patients from the TCGA database. Among them were methylation data from 507 samples, including 32 normal samples and 475 cancer samples, as well as mRNA expression data from 594 samples, including 59 normal samples and 535 LUAD samples. First, based on the LIMMA package, the downloaded data were normalized and analyzed for differences to obtain aberrant methylated genes and differentially expressed genes. Then, based on the MethylMix algorithm implemented by R, we calculated the correlation between gene methylation level and gene expression. Next, we determined genes that were significantly related and identified the disease-specific hypomethylation and hypermethylation genes by constructing the β-mixed model. Finally, screening for methylation-driven genes was done. In addition, we screened 244 samples with stage I LUAD,which both have expression and clinical information for further testing. The data provided by TCGA is public and did not require the approval of a local ethics committee.

Functional and pathway enrichment analysis of methylation-driven genes

The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (http://david.abcc.ncifcrf.gov/) serves as an open source platform for determining associations between target molecules [22]. ConsensusPathDB (http://cpdb.molgen.mpg.de/) integrates interaction networks in Homo sapiens, including binary and complex signaling, gene regulatory and drug-target interactions, as well as biochemical pathways [23, 24]. To gain insight into the biological functions of these methylation-driven genes, the genes were subjected to functional and pathway enrichment analyses based on DAVID and Consensus PathDB online software, and P < 0.05 was set as the cutoff criterion.

Construction of risk assessment model and risk score calculation

To further screen for prognosis-related methylation-driven genes, a linear risk assessment model for LUAD methylation-driven genes was constructed using a single-factor, multivariate Cox analysis [25]. The prognostic index was defined as follows:

In the formula, Exp is the expression level of each methylation-driven gene in the specimen, and β is the multi-factor COX regression analysis coefficient of each methylation-driven gene in the COX model. The prognostic risk value of each sample was calculated according to the formula, and then the median of the index value was cut off. The patients were divided into high and low risk groups [26, 27] and verified with a time-dependent ROC curve. The Kaplan–Meier survival curve method was used to evaluate the overall survival rate of patients in the high- and low-risk groups. Log-rank test was used to determine whether there was any difference in the overall survival rate between the high-risk and low-risk groups. P < 0.05 was considered statistically significant. In addition, we applied the risk prediction model to patients with stage I LUAD who had both the expression and clinical information provided by the TCGA database to further test the validity and practicability of the model to predict prognosis.

Mapping of Kaplan–Meier curves of driver genes and methylated sites in survival models, joint survival analysis

To further explore the prognostic assessment of methylation-driven genes, we extracted relevant loci for driving gene methylation from downloaded LUAD methylation data. Then, based on the survival R package, the prognostic survival analysis of the driving genes and related methylated sites was performed by combining the clinical data and prognostic information of LUAD in TCGA, and the Kaplan–Meier curve was produced. In addition, we performed a joint survival analysis of methylation levels and gene expression levels of driver genes to further identify key genes associated with prognosis in patients with LUAD, and the joint survival curve was also obtained by the survival R package.

Results

TCGA data analysis and acquisition of methylation-driven genes

The methylation data from the analysis of this study were from 507 sample data, including 475 LUAD samples and 32 paracancerous control samples. Gene expression data were gathered from 594 samples, including 535 cancer samples and 59 paracancerous control samples. Abnormal methylation expression and gene expression data of LUAD from the TCGA were extracted and analyzed based on the LIMMA software package. The relevant data were then integrated under the same sample based on the MethylMix package for correlation analysis. A mixed model construction and Wilcoxon rank test for differential methylation were calculated, where |logFC| > 0, P < 0.05, |Cor| > 0.3. As a screening condition, 118 methylation-driven genes were obtained (Fig. 1; Additional files 1, 2, 3).

Fig. 1

Heat maps of LUAD-related aberrant methylation-driven genes. The color from green to red shows a trend from hypomethylation to hypermethylation

Functional enrichment and pathway analyses of methylation-driven genes

To further explore the molecular mechanism of methylation-driven genes in the development of LUAD, we performed a functional enrichment analysis and pathway analysis of these genes using DAVID and ConsensusPathDB online software. The results showed that the methylation-driven genes were enriched not only in multiple pathways but also in molecular functions and biological processes. Pathway analysis showed that these genes were mainly enriched in the generic transcription, RNA polymerase II transcription, gene expression (transcription) and platinum pathways, as well as pharmacokinetics/pharmacodynamics (Fig. 2; Additional file 4). Functional analysis revealed that, in the biological processes (BP) group, these genes were mainly involved in transcriptional regulation, DNA-templates, multicellular organism development, angiogenesis and blood vessel maturation, among others. The molecular function (MF) was mainly enriched in transcription factor activity, sequence-specific DNA binding, nucleic acid binding, and DNA binding. In addition, the cellular component (CC) group was mainly involved in tracellularly (Table 1).

Fig. 2

The significant enriched pathways of methylation-driven genes. Only the pathways which P < 0.01 were shown here. Node size: the number of genes; Node color: P-value; Edge width: percentage of shared genes; Edge color: genes from input

Table 1

Functional enrichment analysis of methylation-driven genes associated with LUAD

Category	Term	Count	P value
GOTERM_BP_DIRECT	GO:0006355~regulation of transcription, DNA-templated	26	1.49E−08
	GO:0006351~transcription, DNA-templated	23	8.28E−05
	GO:0001525~angiogenesis	6	0.004235
	GO:0072272~proximal/distal pattern formation involved in metanephric nephron development	2	0.009506
	GO:0032526~response to retinoic acid	3	0.016299
GOTERM_MF_DIRECT	GO:0003700~transcription factor activity, sequence-specific DNA binding	17	1.77E−05
	GO:0003676~nucleic acid binding	16	9.29E−05
	GO:0046872~metal ion binding	21	0.002409
	GO:0003677~DNA binding	17	0.007796
GOTERM_CC_DIRECT	GO:0005622~intracellular	16	0.001325

If there were more than five terms in this category, selected the first five terms based on the P value. Count: the number of enriched genes in each term

Construction and analysis of the prognosis risk assessment model of LUAD methylation-driven genes

We performed univariate and multivariate Cox regression analyses of the obtained LUAD methylation-driven genes. The results showed that the assessment model constructed by five genes (CCDC181, PLAU, S1PR1, ELF3, KLHDC9) can be used as an independent indicator to predict the prognosis of the disease. The prognostic index = (1.571 * expression level of CCDC181) + (− 1.170 * expression level of PLAU) + (3.674 * expression level of S1PR1) + (3.467 * expression level of ELF3) + (3.028 * expression level of KLHDC9). In addition, with a median PI value (value = 0.928) as a group cutoff condition, 488 samples in the methylated data sample that matched the clinical follow-up sample were divided into a high-risk group (n = 244) and a low-risk group (n = 244). Kaplan–Meier survival curve analysis of patients in the high- and low-risk groups showed that the overall survival rate was lower in the high-risk group, and the difference between the two groups was statistically significant (Fig. 3a). Using the time-dependent ROC curve to estimate the predictive performance of the risk scoring model, the AUC of the prognostic risk assessment model for the five methylation-driven genes was 0.66 at 3 years of OS (Fig. 3b). In the prognostic assessment test for patients with stage I LUAD, the Kaplan–Meier curve showed that the survival rate of the high-risk group was significantly lower than that of the low-risk group (Fig. 4). It can be seen that the model based on methylation-driven genes has certain reliability and practicability in evaluating the prognosis of LUAD patients.

Fig. 3

Kaplan–Meier and ROC curves for a linear risk model based on five methylation-driven genes. a The differences between the high-risk (n = 244) and low-risk (n = 244) groups were determined by the log-rank test. b Time-dependent ROC curves analysis for 3-year survival prediction by methylation-driven genes

Fig. 4

Kaplan–Meier curve for patients with stage I LUAD using the prognostic risk model

Survival analysis of driving gene methylated sites and joint survival analysis of methylation and expression

Prognostic survival analysis of driver gene-related methylated sites in a risk assessment model was performed based on the surviving R package. P < 0.05 was used as a meaningful cut-off condition for predicting prognosis, and we found specific methylated sites associated with the prognosis of these genes (Additional file 5). Among them, 12 methylated sites of the gene CCDC181, 1 methylated site of S1PR1, and 4 methylated sites of KLHDC9 were significantly associated with the prognosis of LUAD (Fig. 5). In addition, joint survival analysis revealed that the combination of methylation and expression of the genes CCDC181, S1PR1, KLHDC9, and PLAU had a significant correlation with the prognosis of the patient, including the gene CCDC181. The high-methylation low-expression survival rate of S1PR1 and KLHDC9 were higher, while the low-methylation high-expression survival rate of the gene PLAU was higher (Fig. 6).

Fig. 5

Kaplan–Meier survival curves of the related methylated sites. a–l Methylated sites of the gene CCDC181; m–p methylated sites of the gene KLHDC9; q methylated site of the gene S1PR1

Fig. 6

Kaplan–Meier survival curves for the joint survival analysis. a The combination of gene CCDC181 methylation and expression; b the combination of gene PLAU methylation and expression; c the combination of gene S1PR1 methylation and expression; d the combination of gene KLHDC9 methylation and expression

Discussion

LUAD is the most common type of lung cancer, and it is highly invasive and has a poor prognosis [28, 29]. In-depth studies on the molecular pathogenesis of LUAD and the early detection of the prognostic markers of this disease as well as specific driving genes are urgent and could prove useful for improving patients’ quality of life and the prognosis of patients. In recent years, it has been found that the decreased expression of genes caused by the hypermethylation and the enhanced expression of genes caused by hypomethylation play an important role in the regulation and development of malignant tumors [9, 30]. Gene methylation could lead to transcriptional disorders, causing some gene expression disorders and cell differentiation disorders [11, 31]. Studies have shown that, unlike DNA aberrations, epigenetic changes are reversible, making them attractive therapeutic targets [32]; thus, detecting and altering DNA methylation can provide new insights for the further treatment and risk assessment of cancer. For instance, Robles et al. found that methylation of the HOXA9 promoter was associated with a high risk in patients with stage I LUAD [33], whereas Sugimoto et al. found that the aberrant methylation of GRWD1 may be a protective factor in tumor development. Because the high expression of GRWD1 in tumor cells can promote tumor cell growth, the aberrant methylation of the GRWD1 gene can inhibit its activity [34]. From previous studies, the prognosis-related methylation genes and specific methylation-driven genes could serve as the new markers for further clinical applications.

In the present study, our aim was to screen for methylation-driven genes by analyzing data on LUAD methylation in TCGA and to identify new prognostic biomarkers associated with methylation. Firstly, we obtained 118 methylation-driven genes using the LIMMA software package. And then, to further investigate the cellular mechanisms involved in these methylation-driven genes, functional enrichment and pathway analyses were performed. Functional aggregation analysis revealed that these genes are involved intracellularly and that major enrichment is associated with the regulation of transcription, DNA-templates, multicellular organism development, and a variety of other functions, such as transcription factor activity and sequence-specific DNA binding. Enrichment pathway analyses indicated that the methylation-driven genes were mainly involved in the generic transcription, RNA polymerase II transcription, gene expression (transcription) and platinum pathways. These functional enrichment and pathway analyses demonstrate not only the cellular mechanisms by which DNA methylation abnormalities lead to disease development but also the functional level interaction of these genes.

In addition, according to univariate and multivariate Cox analyses of methylation-driven genes, five methylation-driven genes (CCDC181, PLAU, S1PR1, ELF3, KLHDC9) were identified and used to construct a prognostic risk model for LUAD. Using this risk model, we successfully divided the LUAD samples into high-risk and low-risk groups. Survival analysis showed significant differences in overall survival between the two groups. The results suggest that the risk model, consisting of five methylation-driven gene profiles, can effectively predict the prognosis of patients with LUAD. The results also show that the AUC of the ROC curve predicting 3-year survival by the five-gene feature was 0.66. Our study found that the risk assessment models constructed by these five gene signatures performed well in survival predictions for patients with LUAD, but further studies are needed to validate these findings. These methylation-driven genes can be used as effective biomarkers or drug targets for early diagnosis, and prognosis of patients with LUAD. C1orf114, also known as CCDC181, is a protein of unknown function and encodes a coiled-coil domain containing 181. Haldrup et al. showed that CCDC181 (C1orf114) and HAPLN3 can be used as novel biomarkers for the hypermethylation of prostate cancer (PC), proving to be a novel diagnostic and/or prognostic promoter. Moreover, the combination of C1orf114 and two other genes in the dichotomy methylation profile (AOX1 and HAPLN3) can increase the prognostic value of PC [35]. Our study found that CCDC181 is one of the methylation-driven genes of LUAD, and its hypermethylation leads to poor prognosis. The molecular mechanism of the gene CCDC18 in cancer requires further research to discover its driving genetic role in LUAD patients.

The gene PLAU is a urokinase plasminogen activator. Based on the established role of the plasminogen system in cancer invasiveness and as a key player in cancer metastasis and cancer cell invasion behavior (adhesion, migration and invasion) [36, 37], PLAU has different expression levels in individual cells of lung cancer patients and can be used as a marker to predict the prognosis of lung cancer patients [38]. Methylation and gene expression combined survival analyses found that the gene PLAU were closely associated with poor prognosis in patients with LUAD. Sphingosine 1-phosphate receptor 1 (S1PR1) is a G protein-coupled receptor involved in the regulation of physiological processes such as cell growth, differentiation, and migration as well as the immune response [39]. Previous studies have shown that S1PR1 can participate in the proliferation and invasion of cancer cells by activating the ERK signaling pathway [40]. Moreover, S1PR1 can stably activate STAT3 in tumor and bone marrow cells [41, 42], which is essential for breast cancer cell proliferation, invasion and metastasis. In addition, Yoshida et al. showed that S1PR1 expression was associated with prognosis in patients with glioblastoma and that high expression was positively correlated with favorable survival [43]. This study found that the high methylation level of S1PR1 was associated with a lower survival rate in patients with LUAD and may be a new target for the treatment and improvement of the prognosis of patients with LUAD.

In addition, as a potential oncogenic transcription factor, ELF3 promotes tumor cell growth and metastasis by modulating the PI3K/Akt and ERK pathways in NSCLC and may serve as a promising new target for the treatment of NSCLC patients [44, 45]. Other studies have shown that ELF3 overexpression is significantly associated with poor prognosis in patients with hepatocellular carcinoma (HCC) [46]. In our study, the gene KLHDC9 was used as one of the prognostic-related methylation-driven genes in LUAD, and the KM curve showed that KLHDC9 hypermethylation patients had shorter survival times and shorter methylation survival times. The combined survival analysis also found that hypermethylation and low gene expression of KLHDC9 suggested poor prognosis in patients with LUAD. There are few studies that have examined the gene KLHDC9, but this gene could be used as a new target for the diagnosis and prognosis of LUAD. However, based on the fact that DNA methylation mainly occurs at the CpG site of the gene, this study found that multiple CpG sites were closely related to the prognosis of patients with LUAD by analyzing the correlation of CpG site survival in methylation-driven genes. These results are based on bioinformatic analysis and provide a theoretical basis for further experimental validation.

Conclusion

In this study, based on the genomic methylation data provided by TCGA for LUAD patients, we obtained 118 methylation-driven genes associated with LUAD using the MethylMix algorithm. Univariate and multivariate Cox regression analyses showed that the prognostic survival model constructed from five aberrant methylation-driven genes, CCDC181, PLAU, S1PR1, ELF3, and KLHDC9, was an independent predictor of disease prognosis including. Based on the risk model of these five methylation-driven genes, LUAD patients can be divided into a high-risk group and low-risk group, and this approach provides a basis for prognosis prediction and personalized treatment plans for LUAD patients. And the predictive ability of the model was further demonstrated in the test of stage I LUAD samples. Since DNA methylation is one of the most important modifications in epigenetics, gene expression can be controlled in a variety of ways. Joint survival analysis revealed that the expression levels of the genes CCDC181, PLAU, S1PR1 and KLHDC9 can be used as independent prognostic markers or drug targets for early diagnosis and prognosis assessment of LUAD. Although further experimental verification is needed, our findings provide an important bioinformatic basis and relevant theoretical basis for guiding the subsequent in-depth study of LUAD.

Additional file 1. Methylation levels of 118 driver genes.

Additional file 2. Expression levels of 118 methylation-driven genes.

Additional file 3. The protein–protein interaction network of 118 methylation-driven genes and the related genes. Red: the methylation-driven genes; Green: the related genes. The methylation-driven genes which not associated with other genes are not shown here.

Additional file 4. Pathway analysis of methylation-driven genes associated with LUAD (P < 0.05).

Additional file 5. All relevant methylated sites of the five methylation-driven genes obtained from the TCGA database. (1–16) methylated sites of the gene CCDC181; (17–33) methylated sites of the gene ELF3; (34–47) methylated sites of the gene KLHDC9; (48) methylated site of the gene PLAU; (49–57) methylated sites of the gene S1PR1.