An analysis about heterogeneity among cancers based on the DNA methylation patterns

Abstract

Background

It is generally believed that DNA methylation, as one of the most important epigenetic modifications, participates in the regulation of gene expression and plays an important role in the development of cancer, and there exits epigenetic heterogeneity among cancers. Therefore, this study tried to screen for reliable prognostic markers for different cancers, providing further explanation for the heterogeneity of cancers, and more targets for clinical transformation studies of cancer from epigenetic perspective.

Methods

This article discusses the epigenetic heterogeneity of cancer in detail. Firstly, DNA methylation data of seven cancer types were obtained from Illumina Infinium HumanMethylation 450 K platform of TCGA database. Then, differential methylation analysis was performed in the promotor region. Secondly, pivotal gene markers were obtained by constructing the DNA methylation correlation network and the gene interaction network in the KEGG pathway, and 317 marker genes obtained from two networks were integrated as candidate markers for the prognosis model. Finally, we used the univariate and multivariate COX regression models to select specific independent prognostic markers for each cancer, and studied the risk factor of these genes by doing survival analysis.

Results

First, the cancer type-specific gene markers were obtained by differential methylation analysis and they were found to be involved in different biological functions by enrichment analysis. Moreover, specific and common diagnostic markers for each type of cancer was sorted out and Kaplan-Meier survival analysis showed that there was significant difference in survival between the two risk groups.

Conclusions

This study screened out reliable prognostic markers for different cancers, providing a further explanation for the heterogeneity of cancer at the DNA methylation level and more targets for clinical conversion studies of cancer.

Background

Recently, cancers are found to have become serious threats to human health. Through epidemiological study, experiments and clinical observations, researchers found that environment and behavior have significant effects on the occurrence of human malignant tumors. All kinds of environmental and hereditary carcinogenic factor can work in a synergistic or orderly manner in the induction of non-lethal DNA damage in cells, which leads to the activation of oncogenes and/or the inactivation of tumor suppressing genes. Moreover, substantial omics heterogeneity has been revealed for histologically homogeneous tumors in terms of genomics [1, 2], epigenomics [3], transcriptomics [4,5,6] and proteomics [7]. Actually, epigenetic modification plays an important role in the development of cancers.

Previous study has proved that epigenetic modification stands for the intersections of genes and environment [8,9,10]. Epigenetic modification can regulate the expression of genes without altering basic DNA sequence [8]. Despite increasing evidence which shows that epigenetic modifications are sensitive to environmental exposure (such as nutritional factors), the influence on epigenetic markers cast by genetic mutation has been spotted [11]. One of the most common epigenetic modifications is DNA methylation. It occurs when methyl is added to specific DNA base pairs, primarily in the background of cytosine dinucleotide (CpG). DNA methylation has been well explored and demonstrated to play essential roles in cellular processes such as regulation of gene expression [12]. According to the place where methylation takes place (such as genome and CpG islands) [13] and the level of DNA methylation, two classes are created, hypomethylation and hypermethylation. There are several most common used ways to analyze the patterns of DNA methylation: global, epigenetic genome range and candidate gene DNA methylation analysis.

Cancer is a type of disease with great genetic and epigenetic heterogeneity. So far, there have been lots of studies that confirm the feasibility of analyzing the epigenetic heterogeneity of cancers using DNA methylation patterns. For instance, it has been proved that DNA methylation heterogeneity is related to Prostatic Carcinoma [14], Low-stage Glioma [15], Esophageal Squamous Cell Carcinoma [16], and the clone of Hepatocellular Carcinoma [17]. In addition, new indicators of DNA methylation heterogeneity, such as epiallele load, Inconsistent Methylated Read Ratio and DNA Methylation Inference Regulatory Activity, are related to the clinical variables of Acute Myeloid Leukemia [18], Chronic Lymphoblastic Leukemia [19] and Sarcoma [20]. However, these researches are all based on the heterogeneity analysis of a single type of cancer, it is also required for a pan-cancer heterogeneity analysis from the global perspective.

This study analyzes the heterogeneity of seven TCGA cancers based on DNA methylation level. We first define specific differentially methylated genes in these cancers. Then, we build methylation-correlation network and KEGG pathway network to sort out pivotal genes and find out cancer-specific methylation markers and prognostic markers. This research can provide clinicians and researchers with more therapeutic and experimental targets, and deeper understandings on cancer heterogeneity.

Methods

Acquisition and preprocessing of DNA methylation data

DNA methylation data of seven cancer types, including 337 COAD (colon adenocarcinoma) samples, 492 LUAD (lung adenocarcinoma) samples, 415 LUSC (lung squamous cell carcinoma) samples, 195 PAAD (pancreatic cancer) samples, 202 ESCA (esophageal cancer) samples, 888 BRCA (Breast invasive carcinoma) samples, 478 UCEC (Uterine Corpus Endometrial Carcinoma) samples, were downloaded from the TCGA (The Cancer Genome Atlas) database, Illumina Infinium HumanMethylation450 BeadChip platform. Specific sample information for each cancer type was shown in Table 1.

Table 1 The sample size for each cancer type

Some pre-processing is conducted on the DNA methylation data. We have removed samples with multiple missing values and recalculated missing values of remaining samples with the function impute. Knn (), R package. We also removed the unstable loci in genome, including CpG loci on sex chromosome, single nucleotide polymorphisms, and CpG loci corresponding to multiple genes. Since the methylation of CpG loci on the promotor region has a strong regulatory effect on gene expression, we only select the CpG loci in the promoter region of genes for further analysis. Here, the promoter region of the gene is defined as the upstream 2 kb region of the transcription initiation site to the downstream 0.5 kb region.

The chip HM450K checks the methylation level of over 480,000 CpG loci in the whole genome. Therefore, chances are that multiple CpG loci are tested in a single gene. Sometimes, differences are huge among those CpG loci which correspond to the same gene, so it’s not reasonable for all the genes we study, to use the average methylation level of those CpG loci to represent the methylation level of the gene. Zhang et al. propose that most of the CpG loci are hypermethylated or hypomethylated (β > 0.5 or β < 0.5) [21], hence in a single sample, we believe that the CpG loci on a gene (gene A) are of the same pattern if all of their β values are greater than or equal to (or less than) 0.5. It is reasonable for us to use the average methylation level of all the CpG loci on gene A to represent the methylation level of this gene if the ratio of samples of the same pattern reach a specific threshold. For genes don’t meet the condition, we remove them from the subsequent analysis. Finally, we use the average methylation level of all the CpG loci on a gene to represent the methylation level of it.

Differentially methylated genes identification per cancer

DNA methylation is the most extensively documented epigenetic modification that can influence cell fate and gene expression [22], which finally leads to the inhibition of gene expression through formation of heterochromatin in the gene regulatory region [23]. In this study, identification of differentially methylated genes in cancer samples and adjacent control samples for all seven types of cancer are our first task.

We use user-defined R script, the bilateral t-test, to recognize the differentially methylated genes among sample pairs. Benjamini-Hochberg method is used in multiple tests to adjust the P value. The gene whose adjusted P value is less than 0.05 and the difference of the average of β is more than 10% is considered distinctly differentially methylated gene among sample pairs.

Biological functions and pathways enrichment analysis of differentially methylated genes

In this study, using DAVID [24, 25], we conduct a GO (Gene Ontology) biological functions enrichment analysis and a KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways enrichment analysis towards the list of differentially methylated genes from the seven cancer types (hypermethylated and differentially hypomethylated genes are also included), with p controlled within 0.05, which could find out the biological characteristics and senses related.

Construction of correlation network of differentially methylated genes

In this study, Pearson correlation coefficient is used to measure the correlation of DNA methylation level of differentially methylated genes of each cancer type quantitatively. The formula is as follows:

$$ \mathrm{r}=\frac{1}{n-1}\sum \limits_{i=1}^n\left(\frac{X_i-\overline{X}}{\delta_X}\right)\left(\frac{Y_i-\overline{Y}}{\delta_Y}\right) $$

T test is used to perform a hypothesis test towards correlation coefficient. In addition to that, we also use permutation test to examine the correlation between DNA methylation levels in each pair of genes. Script of python and R are used to complete the process, and we use the function cor.test() in R for calculation and test of correlation coefficient. We build a methylation correlation net. This net is built, analyzed and visualized using Cytoscape 2.8.2 [17] (http://www.cytoscape.org/). The statistical and functional significance of the network, is proposed to be measured using various statistical parameters, namely in the proposed case, degree (the number of edges per node) and average clustering co-efficient C(k), the ratio of the number of edges E of the node having a k degree with neighbors to the total possible number of such edges.

In the DNA methylation correlation network, different nodes are of different importance, for those whose degrees are large, they often are pivots of the network with lots of genes related to them. If they go abnormal, vertexes adjacent to them will be affected, leading to dysfunction of the pathway and causing cancer. We assume that those key nodes may be associated with the prognosis of cancer patients, thus we pick the top 20% nodes in the network as candidate genes for further analysis. It is also necessary to analyze the interaction information in the pathways of these differentially methylated genes from a functional perspective. DAVID online bioinformatics tools are used in the enrichment analysis of the pathways and functions involved those genes. The result is visualized using EnrichmentMap function in Cytoscape.

Construction of KEGG pathway network of differentially methylated genes

In this study, XML format files of pathways enriched by differentially methylated genes in each cancer type are obtained from KEGG database. User-defined Perl script is used, <relation></relation> block is used to find the molecular interaction pairs within each pathway, <entry></entry> block is used to obtain information about the specific genes or compounds of each pair. Among all those interaction pairs, the interactions of real proteins are our only concern, therefore only ‘PPreal’ type of interaction pairs in relation remain undeleted. Then, the resulting interaction id is then converted into a gene symbol to facilitate visualization and analysis. The network of KEGG pathway is also built by Cytoscape 2.8.2 [26] (http://www.cytoscape.org/) to analyze and visualize the network.

We also pick the top 20% nodes whose degrees are the biggest as candidate genes for further analysis and have a discussion on the functions of those genes. DAVID online bioinformatics tools are used to conduct an enrichment analysis on the pathways and functions in which those genes are involve, the results are visualized using EnrichmentMap [27] in Cytoscape.

The construction of prognostic model and survival analysis

In order to be accurate, all cancer patients in each cancer type were divided into two data sets on average in this study, a training set and a test set. The training set is used for establishment of models and screening of prognostic markers while the test set is for follow-up validation of screened prognostic markers. The division of two sets should meet the following criteria: (1) All samples are divided into training set and test set randomly. (2) There were no significant differences in age distribution, staging, follow-up time and mortality between the two sets (Use Fisher’s exact test or t test). That is to say, patients of all types were randomly assigned to the training and test sets, including patients with missing clinical information. Then, we use the samples of each cancer in the training set and the differentially methylated candidate prognostic markers in each type of cancer obtained from correlation network and the KEGG pathway network to construct a model to screen for specific prognostic markers in cancers.

In the first step, we find out DNA methylation spectrum of candidate markers for each cancer type, as well as clinical phenotype information and follow-up information of the samples and establish a univariate COX proportional risk regression model, so as to assess the association between patient survival and DNA methylation levels. Additionally, we also construct univariate COX proportional risk regression models to determine the clinical factors that significantly affect patient survival. In the next step, significant genes in each cancer type and the clinical factors that significantly affect survival in this cancer type are introduced into the multivariate COX proportional risk regression model to find independent prognostic factors (genes). For each gene i, the formulas of univariate and multivariate COX proportional risk regression models are defined as follows:

$$ {\displaystyle \begin{array}{c}\mathrm{h}{\left(t,x\right)}_i=h0(t)\exp \left({\beta}_{methy} meth{y}_i\right)\\ {}\mathrm{h}{\left(t,x\right)}_i=h0(t)\exp \left({\beta}_{methy} meth{y}_i+\sum {\beta}_{clinical} clinical\right)\end{array}} $$

In the formula, methy_i is the DNA methylation level vector of Gene i in all Samples, clinical represents clinical attribute information, β_methy, β_clinical are the coefficients of the regression model. The positive regression coefficient indicates that the increase of methylation level is related to the increase of death risk (risk gene), while the negative regression coefficient indicates that the increase of methylation level is related to the decrease of death risk (protective gene). Univariate and multivariate COX proportional risk regression models are constructed using function coxph() in survival R package.

After univariate and multivariate COX proportional risk regression analysis, independent prognostic markers that are still significant are used to calculate risk scores in the training set. Risk score is a linear combination of DNA methylation level and regression coefficient of these markers, representing different risk levels of patients. The formula is as follows:

$$ \mathrm{Risk}\kern0.17em \mathrm{Score}={\sum}_{i=1}^{\mathrm{n}}{\beta}_i{X}_i $$

In the formula, β_i is the COX regression coefficient of Gene i in the training set, X_i is the methylation level of Gene i, n is the number of genes that have a significant impact on survival. Next, taking the median risk score as the threshold, the patients in the training set are divided into high-risk group and low-risk group. The survival difference between the two groups is analyzed, the overall survival status of patients is estimated by Kaplan-Meier method and the statistical significance of the difference is determined by log-rank test. Functions survfit() and survdiff() in survival R package are used in the process.

Then, the regression coefficients and the threshold of risk score from the training set are directly applied to the test set, and the patients in the test set are also divided into high-risk group and low-risk group. The prognostic differences between the two risk groups were assessed using the same method as in the training set.

Results

Heterogeneity of differentially methylated genes per cancer

In this study, we have compared the number of the genes obtained and the proportion of the rest of the genes at different ratio threshold (Fig. 1). We hope that we could find a ratio threshold which retains as many genes as possible meanwhile improves the accuracy of calculation of gene methylation level. Eventually, we select 70% as the ratio threshold, which guarantees that about 50% of the original genes remain. At last, we use the average methylation of the CpG loci as the methylation level of the gene in further analysis.

Fig. 1

Comparison of the number of genes and the proportion of remaining genes obtained at different ratio thresholds. a. Comparison of the number of genes obtained at different ratio thresholds. b. Comparison of the proportion of remaining genes when different ratio thresholds are used

Through the process mentioned above, we identified 2214 differentially methylated genes in the total seven cancer types. The numbers of hypomethylated and hypermethylated genes are shown in Table 2. The differentially methylated genes are shown in the volcano plot (Fig. 2), which is drawn using ggplot2 R package.

Table 2 The numbers of differential methylated genes in 7 cancer types

Fig. 2

Volcano plot of differentially methylated genes in seven cancers. a Volcano plot of differentially methylated genes in BRCA. b Volcano plot of differentially methylated genes in COAD. c Volcano plot of differentially methylated genes in ESCA. d Volcano plot of differentially methylated genes in LUAD. e Volcano plot of differentially methylated genes in LUSC. f Volcano plot of differentially methylated genes in PAAD. g Volcano plot of differentially methylated genes in UCEC

All those differentially methylated genes are shown in Additional file 1: Figure S1, which indicates the great heterogeneity of differential methylation markers among the cancer types. Besides, we also use heat map to display the methylation level of differentially methylated genes in cancer samples and adjacent control samples (Fig. 3). We utilize the function pheatmap() in pheatmap package of R to create these graphs. It is from those graphs that we can see that each and every one of the differentially methylated genes of all the cancer types is able to separate cancer samples and adjacent control samples clearly.

Fig. 3

Heat map of differentially methylated genes in seven cancers. a Heat map of differentially methylated genes in BRCA. b Heat map of differentially methylated genes in COAD. c Heat map of differentially methylated genes in ESCA. d Heat map of differentially methylated genes in LUAD. e Heat map of differentially methylated genes in LUSC. f Heat map of differentially methylated genes in PAAD. g Heat map of differentially methylated genes in UCEC

Heterogeneity of pathways and biological functions differentially methylated genes involved

From the result of enrichment, we can see that differentially methylated genes in every cancer type are involved in various biological pathways and functions (Additional file 2: Figure S2, Additional file 3: Figure S3, Additional file 4: Figure S4, Additional file 5: Figure S5, Additional file 6: Figure S6, Additional file 7: Figure S7, Additional file 8: Figure S8). It was found that the most enriched gene ontology and KEGG pathways of these seven cancers are olfactory receptor activity, G-protein coupled receptor activity, odorant binding and Olfactory transduction, which have been reported to have association with cancers in previous studies [28,29,30]. At the same time, the distribution shown in Additional file 9: Figure S9 shows that the heterogeneity of biological pathways and functions enriched from differentially methylated genes among various cancer types are great. Specifically, 28 GO functions and 1 KEGG pathways are enriched from differentially methylated genes in two cancer types, 10 GO functions and 2 KEGG pathways are enriched from differentially methylated genes in three cancer types, 5 GO functions are enriched from differentially methylated genes in four cancer types, 2 GO functions are enriched from differentially methylated genes in five cancer types, 6 GO functions are enriched from differentially methylated genes in 6 cancer types, only 8 GO functions and 1 KEGG pathway are enriched from differentially methylated genes in all seven cancer types. The other 93 GO functions and 8 KEGG pathways are cancer specific, which shows that the heterogeneity of biological pathways and functions enriched from differentially methylated genes among various cancer types are great. Even within the same cancer type, differentially hypomethylated genes and hypermethylated could be involved in different pathways and functions. Enrichment pathways and top GO functions are shown in the graph (Fig. 4, Attached Additional file 10: Figure S10, Additional file 11: Figure S11, Additional file 12: Figure S12, Additional file 13: Figure S13, Additional file 14: Figure S14, Additional file 15: Figure S15).

Fig. 4

The enrichment analysis of differential methylated genes in BRCA. a The enrichment analysis of hypermethylated genes in BRCA. b The enrichment analysis of hypomethylated genes in BRCA

Identification and functional analysis of key genes in correlation network

We get 48,816 pairs of gene pairs whose DNA methylation levels are of strong correlation evidently, there are 7345 pairs in BRCA, 5477 in COAD, 5074 in ESCA, 24818 in LUAD, 4587 in LUSC, 9538 in PAAD, 1488 in UCEC.

The net contains a total number of 48,816 edges (Fig. 5). To assess biological significance of the pathway network, topological properties of the network is studied, the average degree of the nodes is 70.953 and the average clustering coefficient is 0.597, and above all, the degree of the network obeys power law distribution (Additional file 16: Figure S16), which indicates that this network conforms to the characteristics of scale-free biomolecular networks, that is, most of the nodes in the net have small degrees, only a small number of nodes have large degrees.

Fig. 5.

DNA methylation correlation network of differentially methylated genes. The nodes in network represent genes, and the edges represent a strong correlation between the two genes. The nodes marked as colors in the legend represent differential methylation of the gene in the cancer type, and a node with multiple color annotations indicates that the gene is differentially methylated in various cancers

According to the degree ranking of the nodes, the first 274 genes are selected, with a maximal degree of 342 and a minimal degree of 137.Then, we discuss the function of those 274 genes, DAVID online bioinformatics tools are used in the enrichment analysis of the pathways and functions involved those genes. The result is visualized using EnrichmentMap function in Cytoscape (Additional file 17: Figure S17A). We can learn from the graph that these genes are significantly enriched in the biological processes related to G-protein-coupled receptor activity and signal pathway, ion channel-related biological processes and the regulation of cell proliferation and differentiation.

Identification and analysis of key genes in KEGG pathway network

We obtain 6120 pairs of gene interactions in BRCA, 6934 in COAD, 4550 in ESCA, 5329 in LUAD, 6968 in LUSC, 2934 in PAAD, 7996 in UCEC.

The network of KEGG pathway is built (Fig. 6, Cytoscape 2.8.2 [17] (http://www.cytoscape.org/). The nodes in the network represent the genes in the pathways enriched by the differentially methylated genes in each type of cancer, and the edges represent the interaction between the two genes in the pathways. The colored nodes represent the gene is differentially methylated for this type of cancer, the gray nodes represent the non-differentially methylated genes extracted from the pathways but the genes that interact with differentially methylated genes. The size of the nodes is marked by the degree of the node, but the colored nodes are larger because different colors are required to be displayed. There are 1628 nodes and 12,765 edges in the network (Fig. 6). To assess biological significance of the pathway network, topological properties of the network is studied, the average degree of the nodes is 15.682 and the average clustering coefficient is 0.131, and above all, the degree of the network obeys power law distribution (Additional file 18: Figure S18), which indicates that this network conforms to the characteristics of scale-free biomolecular networks, that is, most of the nodes in the net have small degrees, only a small number of nodes have large degrees.

Fig. 6.

The KEGG pathway network. The nodes in network represent genes, and the edges represent the interaction of the two genes in the pathways. The nodes marked as colors in the legend represent differential methylation of the gene in the cancer type or a non-differentiated methylation gene obtained from the pathway. C. Enrichment analysis of prognostic marker candidate gene sets

325 genes are selected with a maximal degree of 510 and a minimal degree of 18. Among those genes, 44 are genes differentially methylated in cancers, 281 are acquired from expansion of the pathways.

We also have a discussion on the functions of those genes. DAVID online bioinformatics tools are used to conduct an enrichment analysis on the pathways and functions in which those genes are involve, the results are visualized using EnrichmentMap in Cytoscape (Additional file 17: Figure S17B). Only the most significant enrichment (FDR < 1E-30) entries are shown in the figure, nodes in the graph represent biological functions or pathways where genes are significantly enriched, and the thickness of edge represents the correlation between these functions and pathways, which are measured by the number of shared genes. We can learn from the graph that those genes are significantly enriched in cancer and multiple signaling pathways, as well as metabolic and biosynthetic pathways.

Integration and functional analysis of cancer-specific prognostic candidate marker sets

In this study, we first obtained the key candidate genes in various cancer types at the epigenetic modification level by DNA methylation correlation between genes, and further obtained more candidate genes from the perspective of functional interaction by pathway enrichment analysis.. The candidate gene obtained by these two methods has only one intersection gene (ADCYAP1R1), which is a common differential methylation gene among three cancers, COAD, PAAD and ESCA. The screening of these two complementary modes avoids the omission of the marker gene, and the candidate marker genes obtained by the two methods are integrated together as a basis for screening and analysis of the next specific cancer type prognostic marker. This study only performed a prognostic efficacy analysis of differentially methylated genes in each cancer type, thus removing 281 genes from the pathway that interacted with the differential genes. Finally, 317 differentially methylated genes in these cancers were obtained as prognostic marker candidate gene set.

Functional analysis of these candidate gene sets revealed significant enrichment of genes in sensory organ-related biological processes, many drug metabolisms, and biological processes and pathways for multiple enzyme synthesis (Additional file 17: Figure S17C). Therefore, it is speculated that abnormalities in these genes may lead to dysregulation of related biological processes and pathways, thus inducing cancer.

Identification and analysis of specific prognostic markers per cancer

After the process mentioned above, we described sample information from two datasets for each cancer type in detail in Table 3, and we identify, from the univariate COX regression model, 4 prognostic risk markers for BRCA, 14 for COAD, 10 for ESCA, 7 for LUAD, 5 for LUSC, 16 for PAAD and 31 for UCEC, clinical factors are included as well as gene methylation. You can find information in detail in the attached table below. In the further analysis of multivariate COX regression, in all seven types of cancer, 3 risk genes that independently affecting prognosis of patients are found in BRCA, 6 in CPAD, 5 in ESCA, 2 in LUAD, 3 in LUSC, 11 in PAAD and 19 in UCEC. You can find information in detail in Table 4.

Table 3 Clinical characteristics of patients in the training set and testing set

Table 4 Results of multivariate COX regression analysis

Survival analysis of the two groups of patients of each type of cancer shows that there are significant differences in survival between the two risk groups in all types of cancer (Fig. 7, attached Additional file 19: Figure S19). Further validation based on the reserved test set using the method stated above shows that there are significant differences in survival between the two groups in all the seven types of cancer except ESCA whose p value of significance is 0.0563 (higher than 0.05) (Fig. 7 attached Additional file 19: Figure S19). Although the significance of ESCA does not reach below 0.05, as we can tell from the figure, the two groups of patients can be separated using the prognostic marker genes sifted out. This suggests that the prognostic markers screened out in this study are reliable and can be used to distinguish the high and low risk of patients. And it’s also worth noting that, prognostic markers, in most types of cancer, are specific to this type of cancer. A few exceptions are the one common prognostic marker in BRCA and UCEC (SNORD114.16), SULT1E1 in BRCA and PAAD, SNORD113.5 in COAD and UCEC. SULT1E1 is a protective factor in both BRCA and PAAD, however, The other two markers play opposite roles in the two types of cancer (risk factor and protective factor).

Fig. 7.

Kaplan-Meier survival curve. a Survival curve of BRCA training set. b Survival curve of BRCA test set. c Survival curve of COAD training set. d. Survival curve of COAD test set

After looking through papers, only 4 genes of these prognostic markers have been verified to be relavant with according cancers, including CCL4 [31, 32] in COAD, CACNA2D3 [33, 34] and SMO [35,36,37] in ESCA, and IL23R [38] in LUAD. Other genes have not been tested to be efficient in treating cancer, which may be potential targets for scientists and doctors to further research on them.

Discussion

The heterogeneity of cancers is one of the reasons why cancers are so hard to be cured clinically, therefore, molecular analysis of the mechanism of cancer heterogeneity and screening of cancer-specific diagnostic and prognostic molecular markers are of great importance for clinical treatment. In addition to genetic mutations, DNA methylation is an important epigenetic alteration that can modify gene expression and is commonly perturbed in cancers [39]. So far, DNA methylation is proposed as a molecular biomarker for cancer detection [40] but also as a biomarker for prediction and stratification of patients with risk of distinct clinical outcome and response to therapies [41], which are found abnormal in the early stage of cancer generation which is a stable marker in cancers. It is a severer change in that it affects the transcriptional regulation of genes, which makes it a potentially important marker for early detection, precise treatment and prognosis assessment of cancer. In cancer detection, DNA methylation also has several advantages over somatic mutation analysis, such as high clinical sensitivity and dynamic range. Moreover, the change of DNA methylation pattern is one of the first detectable tumor-specific changes associated with tumorigenesis. Therefore, it is an important research direction to interpret the heterogeneity of cancer from the perspective of epigenetic abnormality.

Yang et al. provides a comprehensive investigation and reveals meaningful cancer common and specific DNA methylation patterns, contributing to a deeper understanding of pan-cancer studies [42]. They discovered a potential tumorigenesis mechanism that involved of three pan-cancer differentially methylated CpG sites (PDMCs) and 62 PDMCs that are significantly associated with patient survival. They also found that cancer-specific DMCs are enriched in known cancer genes and cell-type-specific super-enhancers.

We also conducted a research on pan-cancer analysis from epigenetic perspective. Compared to the study conducted by Yang at al, we first performed a differential methylation analysis of genes (DMGs) and aimed to find reliable prognostic markers for each cancer from gene levels, and made a supplementation of their survival analysis. In this study, the heterogeneity of DNA methylation markers among cancers is discussed in detail by using the large sample DNA methylation data of seven cancers in TCGA database detected by the open available HM450K chip platform. Differential methylation analysis identifies specific and common tumor markers in each type of cancer, which provides more potential targets for cancer diagnosis and experimental researchers. These cancer type-specific tumor markers are also involved in different biological functions and pathways. In the next step, through using two biological molecular networks, DNA methylation correlation network and KEGG pathway network, the marker sets are further optimized and integrated from the perspective of correlation and functional interaction. At last, the specific prognostic markers for each type of cancer are screened out by using the establishment of prognostic model. These markers can classify the risk of patients ideally, and are verified in the test set. The searching of prognostic markers for cancer provides important reference for clinicians to monitor conditions of patients and to alter regimens of treatment in time.

Conclusions

In this study, DNA methylation markers of only 7 cancer types in TCGA are screened out and analyzed, but the method in this study is also applicable to other cancer types. Also, though the preliminary verification of these markers is realized by the compute in this study, which lays a solid theoretical foundation for the reliability of these markers, further experimental confirmation is still a necessity to promote the process in which those molecular markers are put into clinical use.