Comprehensive analysis of genes based on chr1p/19q co-deletion reveals a robust 4-gene prognostic signature for lower grade glioma

Introduction

Diffuse low-grade and intermediate-grade gliomas (which together make up the lower-grade gliomas, World Health Organization grades II and III) is an infiltrative neoplasm of young adults that mostly arise from cerebral hemispheres¹ Despite the recent progress in neurosurgery, radiotherapy and chemotherapy, no great improvement of the surveillance, epidemiology and end results has been reported in the past 10 years²

The prognosis of lower grade gliomas has traditionally been determined by histological type and WHO grade. More recently, molecular markers have received more and more attention³ Chromosomal 1p and 19q (chr1p/19q) co-deletion is considered to be a good prognostic factor in lower grade glioma. About half of patients have chr1p/19q co-deletion, and such patients are often sensitive to radiotherapy and chemotherapy, especially alkylating agents⁴ However, we found that predicting with prognosis using only the chr1p/19q co-deletion was not accurate enough. The survival of patients without chr1p/19q co-deletion ranged from 5 years to less than 1 year. Therefore, the establishment of a more accurate model to predict the prognosis of lower grade glioma patients is very important.

It is welll known that chr1p and chr19q contain genes associated with DNA damage repairing, spindle checkpoint function, apoptosis, WNT signaling pathways, TGF-βsignaling pathways and tumor suppression.^5,6 Chromosome 1p deletion is generally associated with the initiation of carcinogenesis⁷ However, in lower grade glioma, chr1p/19q co-deletion is a beneficial marker to the prognosis, which is contrary to previous knowledge⁸ We hypothesize that the co-deletion of chr1p/19q leads to the loss of important genes during tumor development. The primary objective of this study was to identify “driver” genes from “passengers” on chr1p/19q and to establish a more accurate prognostic model to account for survival prolongation caused by chr1p/19q co-deletion. We used multiple bioinformatics techniques to identify and validate the gene signature based on deletion gene on chr1p/19q. Through this study, we hoped to find more prognostic biomarkers and therapeutic targets for lower grade glioma patients.

Materials and methods

Patients and samples information

RNASeqV2 level3 data from 510 lower grade glioma patients with complete clinical data in The Cancer Genome Atlas (TCGA) cohort were downloaded from TCGA data portal (http://cancergenome.nih.gov/). Fragments per kilobase million (FPKM) normalized expression level was used to quantify gene expressions in this data set. The FPKM sequencing data for 181 lower grade glioma samples in WHO grade II-III from Chinese Glioma Genome Atlas (CGGA) was downloaded from http://www.cgga.org.cn/and used as an independent validation cohort. Microarray dataset GSE16011 with 284 glioma (WHO grade II-III, 125; IV, 159) and eight normal brain samples measured by using Affymetrix U133 plus 2.0 array (GPL8542) were downloaded from Gene Expression Omnibus (GEO) database. The UCSC xena browser (https://xena.ucsc.edu) was used to obtain the chr1p/19q deletion and IDH1 mutation status of TCGA patients.

Establishment of 4-gene-signature

GSE16011 was used for differential expression gene (DEG) analysis. The TCGA-LGG cohort was selected as the training set and the CGGA cohort was selected as the validation set (Table 1). The establishment of 4-gene-signature included the following steps (Figure S1): 1) by performing differential expression gene analysis on the expression profile chip of GSE16011, 1346 genes with significantly higher expression in cancer samples compared to normal samples were obtained. |log₂FC| >1 and adjusted P-value <0.05 served as the screening conditions for significantly differential expression genes. A total of 76 genes encoded by chr1p/19q were selected based on annotation information of the GSE16011. 2) The univariate COX analysis was performed on the 76 expressed genes in the TCGA cohort. A total of 66 significant genes were obtained with P<0.05 as a threshold. 3) The robustness test was performed on the 66 genes using robust likelihood-based survival modeling. The histological type, chr1p/19q deletion status, IDH1 mutation status and WHO grade were included as prognostic risk factors to exclude their effects. After 100 iterations, a total of 23 prognostic-related robust genes were selected. 4) Multivariate COX analysis was performed on 23 genes using the coxph package in R software. Analysis was performed using the maximum likelihood ratio method. Four genes were finally obtained at a threshold of P<0.05.

Table 1 Three lower grade glioma datasets used in this study

Figure S1 Flowchart of the present study. Abbreviations: DEGs, differentially expressed genes; TCGA, The Cancer Genome Atlas; LGG, lower grade glioma; CGGA, Chinese Glioma Genome Altas; GSEA, gene set enrichment analysis.

Risk score calculation

The following linear equation was established by a multivariate COX proportional hazard model to calculate the risk score for each patient: Risk score= expression of IFI44×0.012 + expression of KIF2C×0.068 + expression of GNG12×0.035 + expression of EMP3×0.0200.

Validation of 4-gene-signature

The area under the curve (AUC) values of 1 year, 3 years, and 5 years were calculated in the training set (TCGA) and the validation set (CGGA), respectively. The median of the risk score in TCGA was used as the cut off value to define the high-risk group and the low risk group, and a Kaplan-Meier curve was drawn to visualize the conclusion. Kaplan-Meier curves were plotted in different subgroups to verify the independence of the 4-gene-signature. To clarify whether the proposed 4-gene-signature was more accurate, we compared the index of concordance (C-index) of the published model with the signature model using the time-dependent AUC curve.

KEGG of DEGs on chr1p/19q

Kyoto Gene and Genomic Encyclopedia (KEGG) analysis was used to analyze the biological pathways where there was an obvious enrichment of DEGs on chr1p/19q. At the P-value of <0.05, the KEGG pathway was significantly enriched by these genes.

Gene set enrichment analysis (GSEA)

GSEA was performed using GSEA v2.2.2 (www.broadinstitute.org/gsea) to identify risk score-related gene sets⁹ First, the risk score was divided into four equal parts from high to low, and the highest 1/4 and the lowest 1/4 were selected as the high and low groups. The risk score-related pathway was annotated by using the C2:CP KEGG gene set from MSigDB¹⁰ Six genomes (FDR <0.05) associated with high expression of risk score were selected as enriched genomes.

Statistical analysis

The R package limma was used for differential gene analysis¹¹ The R-package survival and survminer were used for survival analysis and COX regression analysis. The R package pROC was used to plot the ROC curve¹² The R package rbsurv was used for robustness testing. The comparison of different models was performed with the R package pec::cindex¹³ The R package clusterProfiler was used for KEGG analysis. JavaGSEA was used for GSEA analysis.

Results

Identification of cancer-associated genes on the chr1p/19q of lower grade glioma

In order to identify cancer-associated genes on chr1p/19q, we first used the GSE16011 to screen for differential expressed genes (DEGs) between lower grade glioma patients and normal brain. A total of 1,346 up-regulated genes in cancer were screened (Figure 1A). Among them, 56 genes were located on chr1p and 20 genes were located on chr19q (Figure 1B). These 76 genes were significantly higher in chr1p/19q non-co-deletion lower grade glioma compared to normal tissues and chr1p/19q co-deletion lower grade glioma, which suggests they may be associated with prolonged overall survival (OS) in patients with chr1p/19q co-deletion.

Figure 1 Differentially expressed genes between lower grade glioma samples and normal samples in GSE16011. (A) Heat maps show hierarchical clustering of differentially expressed gene expression levels. Each row represents a patient; each column represents a gene. Red dots represent up-regulated genes, and green dots represent down-regulated genes. (B) The volcano plot shows 1,346 differential genes using the limma package with |log2FC| >1 and adj-P<0.05. Red dot indicates a gene that is up-regulated in cancer compared to normal, and green dot indicates a gene that is down-regulated in cancer. The 76 genes encoded by chr1p/19q are labeled.

Kyoto encyclopedia of genes and genomes (KEGG) analysis of DEGs in chr1p/19q

To elucidate the mechanisms of longer survival in patients with chr1p/19q deletion, we performed KEGG analysis on 76 DEGs encoded by chr1p/19q. KEGG analysis showed that the up-regulated gene encoded by chr1p/19q was most significantly enriched in the “cell cycle,” and were also associated with “breast cancer,” “chronic myeloid leukemia,” and “Notch signaling pathway” (Figure 2A). Genes specifically enriched in each KEGG term were visualized (Figure 2B and C).

Figure 2 KEGG results of 76 up-regulated DEGs encoded by chr1p/19q in GSE16011. (A) DEGs enriched in KEGG pathways “cell cycle,” “breast cancer,” “chronic myeloid leukemia,” and “Notch signaling pathway,” etc. Fold enrichment of each KEGG term is indicated by the x-axis and bar color. (B) Hierarchical clustering of gene expression profiles in each KEGG pathway. (C) Chord plots show the relationship between genes and the KEGG pathway.

Establishment of the 4-gene prognosis signature

The model was tested on 510 patients in TCGA-LGG with complete clinical pathology parameters. The chr1p/19q co-deletion status and the IDH1 mutation status were predicted by the UCSC xena website. First, we used COX proportional hazard regression model to compute the univariate analysis of 76 genes, where 66 genes appeared to be associated with OS (P<0.05). We used the rbsurv package in R to further test the robustness of these 66 genes. As shown in Table 2, 23 genes were considered to be robust prognosis-related genes. They were included in the COX multivariate proportional hazard model, and finally four genes were obtained as independent prognostic factors for OS, including IFI44, KIF2C, GNG12, and EMP3 (Figure 3).

Table 2 Survival-associated gene signature screening using forward selection

Figure 3 Forest plot of Cox PH ratios of the 4-gene-signature, with a 95% CI indicated around the values. n=510.

The 4-gene-signature can stratify the prognosis of lower grade glioma

Based on the multivariate Cox model, the prognostic risk score for each patient was calculated based on the coefficient value and the expression level. We divided patients into high-risk and low-risk groups by the median risk score. Patients in the TCGA had a risk score ranging from 0 to 5 and the median risk score of 0.528 (Figure 4A, top panel). More patients survived in the low-risk group than in the high-risk group (Figure 4A, middle panel). The expression levels of the four genes in the high-risk group were significantly higher than those in the low risk group, thus suggesting that they may be the biomarkers of poor prognosis (Figure 4A, bottom panel). These conclusions were further confirmed in the CGGA dataset including 181 patients with lower grade glioma (Figure 4B). The Log-rank text and Kaplan-Meier curves showed that lower grade glioma patients were significantly stratified in two groups by using the median risk score in TCGA as the cut-off value (Figure 4C and D). Since our signature was obtained from the gene on chr1p/19q, we further examined the stratification ability of this signature in patients without chr1p/19q co-deletion. The results showed that in patients with chr1p/19q non-co-deletion, the stratification ability of this signature was still prominent (Figure 4E).

Figure 4 Survival analysis based on 4-gene-signature. (A) The 4-gene-signature-based risk score in TCGA cohort. Upper panel: the risk score plot based on 4-gene-signature. Middle panel: the live/dead state corresponding to the risk score in the upper panel. The X axis represents the patient; the Y axis represents the number of days of survival. Bottom panel: The z-score transformed expression value of each gene in the 4-gene-signature. (B) The 4-gene-signature-based risk score in CGGA cohort. (C-E) Kaplan–Meier plot for overall survival based on the 4-gene-signature in TCGA (C), CGGA (D) and non-chr1p/19q patients in TCGA (E).

Four-gene-signature risk score is an independent predictor of overall survival

Since IDH1 mutation status, histological type, and WHO classification are important for the prognosis of lower grade glioma, it was necessary to determine whether our risk score was an independent prognostic factor for OS. Our 4-gene-signature can separate patients into high-risk group and low risk groups according to different gender (Figure 5A) and age (Figure 5B). Similarly, in the subgroups of different IDH1 status (Figure 5C), histological type (Figure 5D) and WHO level (Figure 5E), the overall survival of the high-risk group was significantly shortened. The above results indicated that the 4-gene-signature was an independent prognostic predictor for OS.

Figure 5 Stratification analysis by different clinical pathological parameters. The Kaplan-Meier curve showed significant statistical differences in overall survival between the high-risk group and the low-risk group in different gender (A), age (B), WHO grade (C), IDH status (D), and histological type (E).

Comparison of the predictive accuracy for OS between the 4-gene-signature and current staging or prognostic factor

Receiver operating characteristic (ROC) curves were used to assess the prognostic ability of our 4-gene signature. The area under the receiver operating characteristic (ROC) curve (AUC) of 1-year, 3-year and 5-year survival rate were 0.858, 0.853 and 0.837 in TCGA, respectively, which was significantly higher compared to other clinical pathological parameters such as chr1p/19q deletion, histology type, WHO classification and IDH mutations (Figure 6A–C). The same conclusion was reached for the CGGA where the AUC of 4-gene-signature in 1-year, 3-year and 5-year survival rate were 0.834, 0.879 and 0.876, respectively. The AUC values of other clinical parameters were below 0.8 (Figure 6D–F). In the chr1p/19q non-co-deletion cohort of TCGA, the existing parameters had limited predictive power for prognosis in this population (Figure 6G–I). However, our 4-gene-signature showed an accurate predictive ability with an AUC value of 0.860, 0.884 and 0.850 in 1-year, 3-year and 5-year survival rate. These results suggested that the signature is a more accurate and useful tool for the prediction of OS in lower grade glioma patients, especially for chr1p/19q non-co-deletion patients.

Figure 6 Receiver Operator Characteristic (ROC) curve for the 4-gene-signature. (A–C) The ROC curve of the 1-year (A), 3-year (B), and 5-year (C) survival rates of 4-gene-signature in TCGA (training set). (D–F) The ROC curve of the 1-year (D), 3-year (E), and 5-year (F) survival rates of 4-gene-signature in CGGA (validation set). (G–I) The ROC curve of the 1-year (G), 3-year (H), and 5-year (I) survival rates of 4-gene-signature in chr1p/19q non-co-deletion patients in TCGA.

Comparing the prediction accuracy of OS between 4-gene-signature and other published models

Next, we compared the predictive power of this 4-gene-signature and other published models in the TCGA and CGGA datasets. The time-dependent ROC curve showed that our 4-gene signature had significantly superior predictive capability compared to previously published models^14–16 (Figure 7A and B). Whether in the TCGA or CGGA dataset, our 4-gene-signature C-index was greater than 0.8, thus indicating a high consistency of the model with patient survival. For patients with a patient survival over 5 years, this model has the strongest predictive power.

Figure 7 Concordance index showing measure of concordance of 4-gene-signature with overall survival of lower grade glioma patients. (A) TCGA-LGG; (B) CGGA.

GSEA

To investigate the possible functions of the 4-gene signature, we did a GSEA in TCGA-LGG based on the risk group calculated by the signature. In TCGA-LGG and CGGA datasets, six pathways, namely, “KEGG_JNK_STAT_SIGNALING_PATHWAY,” “KEGG_MISMATCH_REPAIR,” “KEGG_APOPTOSIS,” “KEGG_CELL_CYCLE,” “KEGG_P53_SIGNALING_PATHWAY,” and “KEGG_DNA_REPLICATION” were significantly enriched in the high risk group (Figure 8A and B). These results suggest that the difference in prognosis between the high-risk group and the low-risk group may be due to the different proliferation and apoptosis of tumor cells, and this different biological behavior may be mediated by the p53 pathway and the JNK-STAT pathway.

Figure 8 Gene Set Enrichment Analysis (GSEA) of risk score results. KEGG pathway, named “KEGG_CELL_CYCLE,” “KEGG_APOPTOSIS,” “KEGG_DNA_REPLICATION,” “KEGG_MISMATCH_REPAIR,” “KEGG_NOTCH_SIGNALING_PATHWAY,” and “KEGG_P53_SIGNALING_PATHWAY” was significantly enriched in TCGA (A) and CGGA (B).

Discussion

The prognosis of lower grade gliomas has traditionally been determined by histologic type and histologic grade, while over recent years the molecular markers have gained increasing attention³ According to Brat et al,³ patients with IDH mutations and 1p/19q co-deficiency have the best prognosis, while patients with only IDH mutations have the second highest prognosis, and IDH wild-type and 1p/19q wild type patients have the worst prognosis. Other prognostic factors include ATRX mutation,^3,17 TERT promoter mutations,¹⁸ CIC loss,^19,20 FUBP1 loss^19,20 and PTEN loss²¹ The 1p/19q co-deletion is a well-recognized prognostic factor in glioma that often leads to high sensitivity to chemotherapy.^22–25 Yet, the specific mechanisms still remain unclear. The TCGA and CGGA datasets have large sample size, detailed clinical pathological information and good data quality, which facilitates the study of the molecular characteristics of chr1p/19q co-deletion in lower grade glioma at the genomic and transcriptome levels.

In this study, we obtained a total of 1,346 DEGs by comparing the cancer with normal brain tissue in GSE16011. We further screened the genes encoded by chromosome 1p/19q, resulting in 76 DEGs. Chromosome 1 is gene-dense, and it contains 3,141 genes and 991 pseudogenes, and many overlapping coding sequences. Chromosomal 1p contains genes associated with DNA damage repairing, spindle checkpoint function, apoptosis, Wnt signaling pathways, and tumor suppression⁵ Chromosome 19 has approximately 64 million base pairs, accounting for more than 2% of the human genome²⁶ Chromosome 1p deletion is generally associated with the initiation of carcinogenesis⁷ Abnormal chromosome 19 may cause tumorigenesis through tumor suppressor gene deletion,²⁷ DNA repair dysfunction,²⁸ abnormal regulation of TGF-β pathway⁶ and NF-κB pathway²⁹ However, in lower grade glioma, chr1p/19q co-deletion is a beneficial prognosis marker, which is not consistent with existing research. We hypothesized that during the tumorigenesis, certain “driver genes” on the chr1p/19q become activated, and their high expression leads to a worse phenotype of tumor cells, thus negatively affecting the patient’s prognosis. Because chr1p/19q contains a large number of genes, a part of these genes may lead to tumorigenesis, while another part may be related to the prognosis of patients.

In lower grade glioma, the genes on chr1p/19q have strong collinearity. In order to obtain stable prognosis-related genes, we performed a robustness test on the genes obtained from univariate COX analysis with P<0.05. We performed 100 iterations with R package rbsurv and finally obtained 23 genes with stable prognosis in the TCGA-LGG cohort. After further incorporating these 23 genes into the multivariate COX proportional hazard regression model, only four genes remained. Regardless of the training set or the validation set, the AUC value of 1 year, 3 years or 5 years was 0.8 or more. Therefore, we hypothesize that these four genes are the “driver genes” in chr1p/19q. Kinesin family member 2C (KIF2C/MCAK, chr1p34) is the member of the motor proteins family³⁰ that has been identified as a novel marker for prognosis in human gliomas³¹ Interferon induced protein 44 (IFI44, chr1p31.1) is part of the type I IFN-inducible gene family that was previously thought to be an inflammation-related gene.^32,33 Currently, the role of IFI44 in cancer remains unclear. Epithelial membrane protein 3 (EMP3, chr19q13.13) that is often inactivated by hypermethylation-mediated transcriptional repression in many types of cancer, was previously thought to be a tumor suppressor gene.^34–36 Recent studies have also reported that EMP3 may promote the development of glioblastoma by activating the transforming growth factor (TGF)-β/Smad2/3 signaling pathway and is associated with poor prognosis of GBM³⁵ G protein subunit gamma 12(GNG12, chr1p31.3) is a member of the G protein family that belongs to the G protein g subunit37 Previous studies have shown that it may be involved in processes such as cytoskeletal function³⁸ and cell growth³⁹

The reasons why the loss of chr1p/19q lead to an increase in overall survival in lower grade glioma patients still appear to be unclear. In the current study, we aimed to reveal the molecular mechanism underlying chr1p/19q co-deletion leading to prolonged survival in patients with lower grade glioma from the bioinformatics perspective. We first performed KEGG analysis and visualization on 76 DEGs encoded by chr1p/19q. The cell cycle pathway was the most significantly enriched pathway, which can explain the prognosis of patients with chr1p/19q co-deletion to a certain extent. Furthermore, based on the risk scores obtained from the 4-gene-signature, we divided the patients into two groups and performed GSEA analysis. Similar to the results of KEGG, the cell cycle was also significantly enriched in GSEA. In addition, the p53 pathway, the DNA damage repair pathway, and the apoptotic pathway were also significantly enriched. Liu et al⁴⁰ have found that the expression of Ki-67 protein (an indicator of cell proliferation) in 1p/19q co-deficient oligodendroglioma is down-regulated Palfi et al⁴¹ have found that glioma cell apoptosis is more pronounced in ch1p/19q co-deletion patients compared with non-co-deletion patients⁴¹ This is consistent with our functional analysis results.

There are some limitations in this study. First, the prognostic risk score established was based on the Illumina RNA-seq platform. If the data comes from other platforms, it needs to be standardized before it can be included in the formula. Therefore, the application of this model is limited. Second, the genes in this model are in general less studied in lower grade glioma, which is why the specific mechanisms remain unclear, and should be further investigated and verified using the experimental approach. Third, the prognostic model only considers genes that are highly expressed in lower grade glioma relative to adjacent cancer, but does not consider the prognostic role of non-differentiated genes. Some molecules may not be highly expressed in cancer, but still affect the prognosis of patients via other means. These molecules were not considered in this study. Fourth, since the number of normal samples is small (eight cases), the differential gene screening may lead to an increase in the number of final DEGs and affect the accuracy of the results.

Conclusions

In this study, we established a 4-gene-signature based on the chr1p/19q co-deletion by employing integrated bioinformatics analysis. We demonstrated that this model has stronger predictive power than traditional grade, molecular biomarkers and published models. In addition, we found that the prolongation of survival in patients with chr1p/19q co-deletion might be related to the changes in pathways such as cell cycle and DNA mismatch repair. The four genes in this signature are promising molecular targets, especially for chr1p/19q non-co deletion patients.