Screening of key methylation-driven genes CDO1 in breast cancer based on WGCNA

2.1 Data acquisition

The RNA-seq transcriptome data, DNA methylation data, and clinical information of breast cancer were obtained from the TCGA database (https://portal.gdc.cancer.gov/) [12]. There were 1222 samples with transcriptome data, including 113 normal samples and 1109 tumour samples, there were 892 samples with DNA methylation data, including 96 normal samples and 796 tumour samples; There were 1097 samples with clinical information.

2.2 Identification of methylation-driven gene

Sorting out the transcriptome data and methylation data, the differentially expressed genes (DEGs) between tumour and normal samples were screened by the limma (version 3.46.0) R package [13]. The cut-off criteria were set as | log2-fold change (FC) | $>$ 1, 0.05, and 5483 DEGs were finally obtained. Of the DEGs, Methylation driven-genes in breast cancer were careened by the MethylMix R package [14] and the cut-off criteria were set as | logFC | $>$ 1, 0.05, and cor 0.3 [15]. Then the obtained methylation-driven gene expression and methylation pattern in tumour and normal samples were visualized by heat map, and the transcriptome heat map and methylation heat map of methylation-driven genes were obtained by the heatmap R package.

2.3 Enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) signal pathway

First, methylation-driven genes were converted into IDs by the org.Hs.eg.db (version 3.12.0) R package. GO functional enrichment analysis of the methylation-driven genes was carried out by using clusterProfiler (version 3.18.1) R package ( 0.05) [16], and the output results were visualized by using the GOplot R package. In the same way, KEGG analysis of the methylation-driven genes was carried out by the clusterProfiler (version 3.18.1) R package ( 0.05) [17], and the output results were visualized by GOplot R package [18].

2.4 Construction of a prognostic models of the methylation-driven genes

First, the expression degree and methylation degree of methylation-driven genes were integrated with survival data, and then the linear risk assessment model of the methylation-driven genes was constructed by univariate and multivariate Cox regression analysis with survival R package ( 0.05) [19]. The risk score(rs) of each patient was calculated by the following formula: r ⁢ s = gene.exp ⁢ 1 × β ⁢ 1 + gene.exp ⁢ 1 × β ⁢ 2 + … ⁢ gene.expn × β ⁢ n , where β is the coefficient of each methylation-driven gene from multivariate Cox regression analysis, gene.exp is the expression of each gene. Then the median value of the risk score was set as the critical value, patients above the median value were included in the high-risk group, and patients below the median value were included in the low-risk groups. The survival and survminer R packages were used to analyse the survival of the high-risk and low-risk groups. The accuracy of the prognostic model in predicting patient survival was verified by drawing multi-index ROC curve with survival ROC R package [20]. ROC can reflect the trend of the model’s sensitivity (FPR) and its specificity (TPR) when selecting different thresholds. When the ROC curve is closer to the upper left corner, the area under the ROC curve (AUC) is larger. The higher the FPR and the lower the FPR of the test, that is, the higher the sensitivity and the lower the false positive rate, indicating that the test is higher diagnostic value [21].

2.5 Identification of key coexpression modules by WGCNA

WGCNA mines genes with similar expression patterns based on different gene expression patterns, which are, included in a module. The module features and hub genes within modules were identified, and each module was associated with other modules or external sample features to screen core genes in modules that can serve as biological markers [22]. In this paper, the weighted gene coexpression network of the breast cancer transcriptome dataset was constructed by using the WGCNA package in R software (version 4.0.3) [23]. To further select the functional modules in the coexpression network, we calculated the correlation coefficient between the modules and the phenotype of interest, obtained the correlation degree between the modules and clinical features, and selected the module with a high correlation coefficient for subsequent analysis, which was considered as a functional module related to clinical features.

The overlapping genes between the co-expressed genes of functional modules with high correlation coefficients obtained by WGCNA and the prognosis-related methylation-driven genes were regarded as key methylation-driven genes.

2.7 Transcriptome and methylation combined survival analysis and expression level verification

To further study the influence of the expression level and methylation degree of the key methylation-driven genes on the prognosis of breast cancer patients, the survival and Survminer R packages were used to analyse the survival rates related in gene expression and methylation, degree of these genes to better predict the prognosis of breast cancer patients [24, 25]. Gene expression profiling interactive analysis (GEPIA) is a public database (http://gepia.cancer-pku.cn/) for the gene expression profile analysis of cancer and normal tissues [26, 27]. The expression level of the key methylation-driven genes in tumour tissues and normal tissues were further verified by the online GEPIA tool.

3.1 Identification of methylation driven genes in breast cancer

The transcriptome sequencing data of 1222 breast cancer patients were obtained from the TCGA database, including 113 normal patients and 1109 tumour patients. A total of 892 samples had DNA methylation data, including 96 normal samples and 796 tumour samples. With | log2 FC | $>$ 1, and 0.05 as the cut-off criteria, a total of 5483 DEGs were screened for further analysis. With | log2FC | $>$ 1, 0.05, anc cor 0.3 as the cut-off criteria, we identified 17 methylation-driven genes (Table 1). The expression level and methylation degree of the 17 methylation-driven genes in tumour and normal samples are displayed by a heat map (Fig. 2A and B).

OPEN IN VIEWER

Seventeen methylation-driven genes were analysed in terms of their GO functions and KEGG signalling pathway to identify related biological functions and signalling pathways. The results showed that the methylation-driven genes were mainly enriched in oligodendrocyte development, oligodendrocyte differentiation, glial cell development, myelination, ensheathment of neurons, and axon ensheathment ( 0.05) (Fig. 3A). The results of KEGG signalling pathway enrichment analysis showed that the methylation-driven genes were mainly enriched in the taurine and hypotaurine metabolism, cysteine and methionine metabolism, malaria, and hedgehog signalling pathways ( 0.05) (Fig. 3B).

OPEN IN VIEWER

To determine the prognostic role of methylation-driven genes in breast cancer, univariate Cox regression was used to analyse the methylation-driven genes related to prognosis, and 4 risk genes with HR $>$ 1 were identified, namely, USP44, CDO1, EVC2, and ZSCAN23. The multivariate Cox regression analysis, results are shown in Table 2. Based on these results, USP44, CDO1, EVC2 and ZSCAN23 were used to construct a prognostic model. The risk score of each breast cancer patient was calculated by the following formula: r ⁢ s = ( 0.230 ) × 𝑈𝑆𝑃 ⁢ 44 ⁢ expression value + ( 0.232 ) × 𝐶𝐷𝑂 ⁢ 1 ⁢ expression value + ( 0.568 ) × 𝐸𝑉𝐶 ⁢ 2 ⁢ expression 𝑣𝑎𝑙𝑢𝑒 + ( 1.823 ) × 𝑍𝑆𝐶𝐴𝑁 ⁢ 23 ⁢ expression value . Then, the median value of the risk score was used as the critical value, and breast cancer patients with complete survival information were divided into a high-risk group ( n = 388) and a low-risk group ( n = 389). Then, the survival curves of the high-risk group and low-risk group were drawn (Fig. 3A). The results showed that the overall survival of the low-risk group was significantly higher than that of the high-risk group ( 0.001). The ROC curve showed that the risk score (AUC = 0.784) had a good predictive ability (Fig. 4B) and could be used for the prognosis of breast cancer. Based on the risk score distribution (Fig. 4C) and the related survival status (Fig. 4D) of the four methylation-driven genes, the high-risk group and the low-risk group could be better distinguished.

OPEN IN VIEWER

The breast cancer transcriptome data set was analysed by WGCNA. First, sample clustering was performed on the expression matrix, and no outliers were found (Fig. 5A). Second, the clinical data were prepared, imported, and compared with the sample cluster tree. The results showed that the sample cluster tree matched the clinical data (normal and tumour) well (Fig. 5B). Then, the appropriate soft threshold β value was screened to construct a network that was more in line with the characteristics of a nonscale network. When the scale-free topological fitting index R2 reached 0.9, the appropriate β value of 6 was selected (Fig. 5C). The minimum module size parameter was set to 50, and the cutting height was set to 0.25 by using the dynamic cutting tree method. The modules that were close to each other were merged into new modules, and 22 gene modules were obtained (Fig. 5D). Different colors represent different modules, each leaf node in the tree diagram represents a gene, and densely connected branches represent close genes. Then, a heat map of the relationship between modules and traits was drawn to evaluate the relationship between two clinical characteristics (tumour and normal) and each module. According to the results of the analysis of the relationships between modules and clinical traits, the green module (Fig. 5E) had the highest correlation with breast cancer (ME green: r = 0.73, p = 3e-203). The green module contains 2010 co-expressed genes. In addition, the coexpressed genes of breast cancer in the green module had the highest correlation coefficients with this module (Cor = 0.84, 0.01) (Fig. 5F).

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

Venn diagram of the 2010 coexpressed genes in the green module and 4 methylation-driven genes related to prognosis was drawn (Fig. 6A). It was found that only CDO1 was an overlapping gene, and it was used as a key methylation-driven gene related to prognosis. The degree of CDO1 methylation in breast cancer patients and normal patients was compared, and the correlation between the degree of methylation and gene expression was determined (Fig. 6B and C). The degree of methylation of the CDO1 gene in tumour samples was higher than that in normal samples, and its methylation degree was negatively correlated with the expression level.

3.6 Combined analysis of the methylation and expression of CDO1 and verification of its expression level

Kaplan Meier curve analysis of the methylation-driven genes showed that the combination of the methylation and expression of CDO1 was significantly correlated with the prognosis of breast cancer patients ( 0.01) (Fig. 7A). Patients with hypermethylation and low expression of CDO1 had a lower survival rate. The expression level of the CDO1 gene in normal tissues and tumour tissues was further verified by GEPIA public database (Fig. 7B). Compared with normal tissues, the CDO1 gene was significantly downregulated in breast cancer. The expression level of CDO1 in various tumours was analysed in the GEPIA public database (Fig. 7C) and was almost lower than that in normal tissues, except in glioblastoma multiforme (GBM) and lower grade glioma (LGG).