Prognostic factor identification by analysis of the gene expression and DNA methylation data in glioma

1.   Introduction

Gliomas are the most malignant and aggressive cancers of the central nervous system originating from glial cells ^[1], the vast majority of which are characterized by diffuse infiltrative growth into the surrounding central nervous system parenchyma ^[2]. Gliomas include oligodendrogliomas, ependymomas and astrocytomas. Malignant gliomas are classified by the World Health Organization as either grade Ⅲ and Ⅳ tumors or grade Ⅳ and Ⅳ tumors ^[3], which are treated with radiation and temozolomide, with only a minor benefit in survival time ^[4]. Therefore, understanding the molecule mechanisms of malignant gliomas and translating this understanding into treatment is crucially needed.

Numerous genes are demonstrated to be genetically altered in gliomas ^[5]. For instance, EGFR, RTEL1, TP53, and TERT have been reported to increase the risk for all types of glioma. CCDC26 and PHLDB1 are implicated in gliomas with isocitrate dehydrogenase mutations ^[6,7]. Recent study has found that deregulation of cell functions in cancer is encoded in both the genome and epigenome ^[8]. DNA methylation has nowadays emerged a key regulator of gene transcription ^[9]. Many studies have reported alterations of DNA methylation in gliomas ^[10,11,12]. Kroes et al. ^[13] found that epigenetic modulation of ST6Gal1 expression played an important role in the glioma phenotype. Majchrzak-Celińska et al. ^[14] reported that SFRP1 promoter methylation could serve as a potential indicator of the survival of glioma patients. We speculated that the genes that are both differentially expressed and methylated in gliomas may play a key role in glioma progression.

In this study, we downloaded the RNAseq and DNA methylation data associated with glioma from GEO and TCGA databases to analyze the differentially expressed genes (DEGs) and methylated genes between tumor and normal tissues. These genes may play a causal role in gliomagenesis, and also have clinical and prognostic importance.

2.   Materials and methods

2.1.   Data acquisition

The glioma associated RNAseq and DNA methylation data ^[1] were download from both GEO database ^[15] (http://www.ncbi.nlm.nih.gov/geo/) and TCGA (https://xenabrowser.net/datapages/) dataset in UCSC Xena database ^[16]. For the RNAseq data (GSE15222) in GEO, there were 187 normal brain tissue samples (normal), and the gene expression was quantified using Sentrix humanref-8 BeadChip (GPL2700). The data in UCSC Xene database were normalized gene expression matrix, which contained 530 tumor samples. For the DNA methylation samples (GSE74486) in GEO, the corresponding platform was Illumina Infinium HumanMethylation450, and 58 normal samples were selected. The DNA methylation dataset downloaded from the UCSC Xena database contained 530 tumor samples. Additionally, the clinical phenotypic data of samples were downloaded.

2.2.   Differential expression analysis

RNAseq data matrixes were normalized with betaqn method in R. Based on the normalized dataset, limma package ^[17] (version 3.10.3) was used for differential analysis of tumor vs. normal, and paired t-test was used for significance test. All genes were tested to obtain the corresponding p values. Benjamini & Hochberg (BH) method was used for multiple test correction, and the corrected p value was adj.P.Value. The threshold values were set as adj.P.Value < 0.01 and |log fold change (FC)| > 2.

2.3.   Functional and pathway enrichment analyses of DEGs

DAVID ^[18] (version 6.7, https://david-d.ncifcrf.gov/) was used for gene ontology (GO) analysis to understand the functions of DEGs in terms of biological process (BP), cellular component (CC) and molecular function (MF). Results were visualized using GOplot ^[19]. Based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database ^[20], pathway analysis was performed using Gene set enrichment analysis (GSEA, version 3.0) ^[21]. Enrichment results with significance threshold of adj.P.Value < 0.05 were screened.

2.4.   Co-expression analysis of DEGs

The Pearson correlation coefficients between the DEGs pairs were calculated and the correlation test was carried out by using corr.test method of psych ^[22] (parameter: ci = F, adjust = "BH"). Multiple test correction was performed by BH method. Gene pairs of |R| > 0.9 and p value < 0.05 were screened, and the co-expression network was constructed using Cytoscape software ^[9] (version 3.7.0).

2.5.   Prognostic gene screening

The prognosis information of the corresponding patients, including overall survival (OS) and OS status, was collected from the downloaded clinical data. The DEGs was used as candidate genes, and their expression medians were used as the boundary value to divide the samples into high expression and low expression groups. Kaplan-Meier (K-M) survival analysis was conducted by combining with the gene expression values and prognostic information, and the K-M curves were drawn. Significance of p value was calculated by Logrank test, and mRNA with p value < 0.05 was preliminarily screened.

2.6.   Differential methylation sites analysis

The methylation data matrixes were normalized with betaqn method in R. Based on the normalized dataset, differential analysis of tumor vs. normal was performed using the limma package followed by significance test using the paired t test. The obtained p values were corrected by multiple tests using BH method, obtaining adj.P.Value. The threshold values were set as adj.P.Value < 0.05 and |log FC| > 1.

The coordinates of the above selected differential methylated chips on the genes (including the corresponding genes, chromosomes, initiation, and termination sites) were downloaded from UCSC Xena database (illuminaMethyl450_hg19_GPL16304_TCGAlegacy), and the corresponding version of the genome annotation was downloaded from GENCODE database (https://www.gencodegenes.org/human/release_19.html). The methylation sites are located on the genes by comparing the coordinates of the methylation sites with those of the genes, including exon, intron, 5’-UTR and 3’-UTR.

The Pearson correlation coefficients between methylated sites and gene expressions were calculated in R, and multiple test correction was performed by Bonferroni method. The threshold values were set as adj.P.Value ≤ 0.05 and |r| ≥ 0.7. In addition, the relationships among gene expressions, degree of methylation and histological types were analyzed by Analysis of Variance (AOV) and Bonferroni methods in R.

2.7.   Key genes identification

The genes that were significantly differentially expressed and differentially methylated between tumor and normal samples were select. We used DAVID database to analyze the functional annotation of these genes, and visualized the results with GOplot tool.

3.   Results

3.1.   Differential expression analysis

The number of annotated genes in TCGA and GEO databases was 20,530 and 16,754, respectively, with a total of 10,678 annotated genes. Among these annotated genes, there were a total of 2190 DEGs between tumor and normal groups, including 1001 up-regulated and 1189 down-regulated genes in tumor group (Figure 1A). The top 10 up- and down-regulated genes are shown in Figure 1B.

3.2.   GO enrichment analysis

GO enrichment analysis was conducted on the DEGs obtained above, and the results showed that 23 BP, 16 CC and 2 MF terms were significantly enriched (adj.P.Value < 0.05). For BP, CC and MF, the results of the top 5 are presented in Figure 2. DEGs were significantly enriched in neuron differentiation, neuron development, axonogenesis, transmission of nerve impulse, cell morphogenesis, translation, ribosome, cytosolic ribosome, cell-cell signaling, and homophilic cell adhesion.

3.3.   Pathway analysis

According to GSEA enrichment analysis of pathways, we identified a total of 54 significantly enriched pathways (adj.P.Value < 0.05), including 40 up-regulated pathways (normalized enrichment score, NES > 0) and 14 down-regulated pathways (NES < 0) (Figure 3A and 3B). KEGG ribosome had the highest NES in the significantly up-regulated pathways (Figure 3C). In addition, the significantly up-regulated pathways also included oxidative phosphorylation, protein export, and DNA replication. The most significantly down-regulated pathway was neuroactive ligand receptor interaction (Figure 3D). Moreover, many signal transduction pathways were significantly down-regulated, such as the Ras signaling pathway, sphingolipid signaling pathway, and CAMP signaling pathway.

3.4.   Co-expression network analysis

For the 2190 DEGs, the correlation coefficients and significance degree of gene pairs were calculated based on their expression levels in different samples. The gene pairs with |r| > 0.9 and adj.P.Value < 0.05 were screened. The visualized network is shown in Figure 4. There were 440 nodes (genes) in this co-expression network, among which 435 genes were up-regulated and only 5 genes were down-regulated. As can be seen from Figure 4, there were 5 subnetworks, among which 98.9% of the genes (435/440) were located in the same subnetwork. In this network, there were some high-connectivity genes (hub nodes), such as dipeptidyl peptidase like 6 (DPP6), mitogen-activated protein kinase 10 (MAPK10), RAB39B, member RAS oncogene family (RAB39B) and ribosomal protein L3 (RPL3), whose co-expressed genes are 285,178,110 and 109, respectively.

3.5.   Gene expression and prognosis

The downloaded TCGA clinical data contained the information of 515 patients, where the survival status and the corresponding survival time were extracted, and the mislabeled cases were deleted. Finally, the clinical information of 514 patients was retained. According to the log FC values of DEGs, the first 10 genes that were up-regulated and down-regulated were selected respectively, and the KM survival curves were plotted by combining the expression values of these genes in different samples and the clinical information (Figure 5). Based on the threshold of p value < 0.05, 18 genes (except DPP6 and VSIG1) were significantly correlated with the survival time of the patients.

3.6.   Differential methylation sites analysis

The methylation sites of TCGA data set and GEO data set were 485,577 and 482,421, respectively, and the number of sites shared by the two datasets was 482,421. The differential analysis of the methylation levels of these sites showed that there were 6186 methylation sites (2834 up-regulated and 3352 down-regulated) with significant differences in tumor vs. normal (Figure 6A). These methylation sites involved a total of 4898 genes, and the intersection (632 genes) of these genes with the DEGs is shown in Figure 6B.

The above 632 simultaneous differentially expressed and DNA differentially methylated genes were involved in 1010 methylation sites and constituted a total of 1116 pairs of methylation site-gene. The results showed that most methylation sites (77%) were located in noncoding regions (UTR and intron), and only 23% methylation sites were located in exon region (Figure 6C). In addition, the number of methylation sites located in the 5 'UTR zone (30%) was significantly higher than that in the 3' UTR zone (7%). In addition, based on the adj.P.Value ≤ 0.05 and |r| ≥ 0.7, a total of 25 pairs of methylation-gene were screened, including 7 positive correlations and 18 negative correlations (Table 1).

3.7.   The correlation among gene expressions, degree of methylation and histological types

Combined with the histological types of TCGA clinical information, a total of 529 valid samples were collected, including 197 Astrocytoma, 198 Oligodendroglioma, and 134 Oligoastrocytoma. After analyzing, 256 gene expression values and methylation at 280 sites were significantly associated with histological types (adj.P.Value < 0.05, Supplementary table 1). Figure 7 showed that top four methylation sites (cg02132760, cg14243481, cg14640066, and cg25766425) and top four genes (CPLX2, MYT1L, NOG, and PSD) were closely related to Astrocytoma, Oligodendroglioma, and Oligoaastrocytoma.

3.8.   Functional analysis of genes corresponding to the differential methylation sites

Functional analysis of these 632 genes showed that a total of 21 BP, 8 CC and 6 MF terms were significantly enriched (Figure 8). These genes were mainly involved in functions of neuron differentiation, neuron development, cell morphogenesis, cell morphogenesis involved in neuron differentiation, cell-cell adhesion, cell-cell signaling, cell junction, transcription regulator activity, transcription factor activity, etc.

4.   Discussions

In this study, 2190 DEGs were identified between tumor and normal tissues, which were significantly enriched in neuron differentiation associated functions, as well as ribosome pathway. In the constructed co-expression network, DPP6, MAPK10 and RPL3 were hub genes. Survival analysis of 20 DEGs obtained 18 prognostic genes, among which 9 were differentially methylated, such as LHFPL tetraspan subfamily member 3 (LHFPL3), cadherin 20 (CDH20), complexin 2 (CPLX2), and tenascin R (TNR). In addition, it was found that most methylation sites (77%) were located in noncoding regions (UTR and intron), and only 23% methylation sites were located in exon region. A total of 25 pairs of methylation-gene were screened, including 7 positive correlations and 18 negative correlations. After analyzing, 256 gene expression values and methylation at 280 sites were significantly associated with Astrocytoma, Oligodendroglioma and Oligoastrocytoma.

DPP6 is a regulatory subunit of the voltage-gated A-type Kv4.2 potassium channel complex expressed in neuronal dendrites and soma ^[23]. DPP6 has been somatically altered in pancreatic cancers and plays a role in pancreatic cancer invasion ^[24]. MAPK10 is a member of the MAPK gene family, which has been found to be implicated in the proliferation and migration process of tumor cells ^[25,26]. A recent study reported that the expression level of MAPK10 gene changed significantly in glioma cells ^[27]. RPL3 encodes a ribosomal protein. It has been suggested that changes in the expression levels of ribosomal proteins have predictive value to distinguish cancer and normal cells ^[28]. A recent study has reported that human RPL3 can act as stress sensing molecule essential for cancer cell response to ribosomal stress in lung cancer cells lacking active p53 ^[29]. In our study, RPL3 was enriched in ribosome pathway, an up-regulated pathway with the highest NES. The enhanced growth of cancer cells requires an increase in global protein synthesis that it is correlated with increased ribosome activity ^[30]. Our result further strengthened the link between ribosome defects and glioma progression. Taken together, we speculated that DPP6, MAPK10 and RPL3 may play important roles in tumorigenesis of glioma.

Previous study had indicated that DNA methylation may be critical to the development of cancers, and have great potential to serve as biomarkers in predicting the prognosis and monitoring response to therapy ^[31]. In our reseach, it was found that most methylation sites (77%) were located in noncoding regions (UTR and intron), and only 23% methylation sites were located in exon region. A total of 25 pairs of methylation-gene were screened, including 7 positive correlations and 18 negative correlations. By analyzing, 256 gene expression values and methylation at 280 sites were closely related to histological types. Additionally, survival analysis of 20 DEGs obtained 18 prognostic genes, among which LHFPL3, CDH20, CPLX2 and TNR were differentially methylated. LHFPL3 is a member of LHFPL-like family and functions as a translocation partner of HMGIC in lipoma ^[32]. Amplification of LHFPL3 has been correlated with mesenchymal differentiation in gliosarcoma ^[33]. Additionally, alterations of LHFPL3 have been reported to be more frequent in high level of genomic instability, and in grade Ⅳ glioma ^[34]. Recently, Li et al. ^[35] detected 20 pairs of glioma tissues and revealed that the expression level of LHFPL3 was significantly higher in glioma tissues compared with in the normal tissues. Study has reported that the survival rate of glioblastoma patients with LHFPL3 mutations is significantly worse than those without the alterations ^[36]. CPLX2 is a member of the complexin/synaphin family, being involved in synaptogenesis and regulation of the neurotransmitter release from pre-synaptic terminals in brain ^[37]. Therefore, CPLX2 plays a key role in the maintainance of normal neurological function ^[38]. Komatsu et al. ^[39] have reported that CPLX2, as an important regulator of neuroendocrine function, is closely related to the growth of high grade lung high grade neuroendocrine tumors. A recent study suggested that CPLX2 was up-regulated in the glioblastoma tissues compared with normal brain tissues ^[40]. Therefore, LHFPL3 and CPLX2 may be served as potential prognostic biomarkers in glioma.

In addition, CDH20 was significantly involved in GO functions associated cell adhesion. Cell-cell adhesion determines the cell polarity and is involves in the maintenance of tissues. Generally, cell-cell adhesiveness is decreased in human cancers ^[41]. Study has reported that changes in the function of cell-adhesion molecules contribute to the progression of tumors both by affecting cell signaling and altering the adhesion status of cells ^[42]. Adaptive adhesion systems have been revealed to mediate the invasion of glioma cells in complex environments ^[43]. TNR is a brain-specific member of the tenascin family comprising tenascins C (TNC), X, and W. TNC is the most prominent member of the tenascin family, and is highly expressed in glioma tissue. Growing evidence has indicated that TNC plays a key role in cell migration or invasion of glioma ^[44]. In this study TNR was enriched in nervous system associated functions, such as neuron differentiation. Recently, a study revealed that a single factor that promotes neuron differentiation can suppress cell growth of glioma cells ^[45], which suggesting the relationship between s neuron differentiation and occurrence of glioma. Therefore, a better understanding of the cell adhesion, migration and invasion molecules, such as CDH20 and TNR, may help to develop new treatment methods.

Although these DEGs and differentially methylated genes were identified through a series of analyses, the results were not validated by animal or clinical experiments. Therefore, further experimental study is needed to confirm our findings.

In conclusion, DPP6, MAPK10 and RPL3 may play important roles in tumorigenesis of glioma. Additionally, methylation of LHFPL3, CDH20, CPLX2, and TNR may serve as prognostic factors of glioma.

Acknowledgement

This work was supported by the Natural Science Youth Foundation of China (Program No. 81400404).

Conflict of interest

The authors declare that there are no competing financial interests exist.