Abstract
Lower grade gliomas are invasive brain tumors that are difficult to completely resect neurosurgically. They often recur following resection and progress, resulting in death. Although previous studies have shown that specific germline variants increase the risk of tumor formation, no previous study has screened many germline variants to identify variants predictive of survival in patients with glioma. In this study, we present an approach to identify the small fraction of prognostic germline variants from the pool of over four million variants that we variant called in The Cancer Genome Atlas whole-exome sequencing and RNA sequencing datasets. We identified two germline variants that are predictive of poor patient outcomes by Cox regression, controlling for eleven covariates. rs61757955 is a germline variant found in the 3′ UTR of GRB2 associated with increased KRAS signaling, CIC mutations, and 1p/19q codeletion. rs34988193 is a germline variant found in the tumor suppressor gene ANKDD1a that causes an amino acid change from lysine to glutamate. This variant was found to be predictive of poor prognosis in two independent lower grade glioma datasets and is predicted to be within the top 0.06% of deleterious mutations across the human genome. The wild-type residue is conserved in all 22 other species with a homologous protein.
Implications: This is the first study presenting an approach to screening many germline variants to identify variants predictive of survival and our application of this methodology revealed the germline variants rs61757955 and rs34988193 as being predictive of survival in patients with lower grade glioma.
Introduction
Grade II and grade III (lower grade) gliomas are primary brain tumors that are derived from glial cells and include astrocytomas and oligodendrogliomas. They are most commonly found in the cerebral hemispheres. They are highly invasive and therefore difficult to completely resect neurosurgically without significant patient morbidity. Following surgery, patients are typically treated with chemotherapy and radiation, although these tumors typically recur or progress to grade IV gliomas and are fatal (1). The median survival following lower grade glioma diagnosis is around 7 years (2).
While the 2007 World Health Organization's (WHO) classification of central nervous system neoplasms differentiated between neoplasms primarily based on histologic features, the updated 2016 WHO classification system now utilizes both molecular and histologic parameters (1). Isocitrate dehydrogenase mutation (IDH) status, 1p/19q codeletion status, telomerase reverse transcriptase (TERT) promoter mutation status, MGMT promoter methylation, TP53 mutation status, and ATRX mutation status may be used to molecularly characterize gliomas (1, 3). The availability of genomic data from patient glioma samples from groups such as The Cancer Genome Atlas (TCGA), the Chinese Glioma Genome Atlas (CGGA), and the Ivy Glioblastoma Atlas Project has substantially contributed to our understanding of these tumors (4, 5).
Many studies have utilized these datasets to identify gene expression signatures, miRNA expression patterns, somatic mutation status, and imaging characteristics that are predictive of survival in lower grade gliomas (6–8). While studies have shown that germline mutations can increase an individual's susceptibility for specific cancers (9–12), including a recent study that identified 853 pathogenic or likely pathogenic germline variants found in 8% of 10,389 patients with cancer (13), no study has comprehensively screened all of the germline variants in a given cancer type to discover the prognostic variants in that cancer type. Although germline mutations have been shown to be prognostic in breast cancer (14) and medulloblastoma (9) in genes that have been well-characterized in the context of these cancers, these variants were not identified using an unbiased approach that screened a large number of germline variants. Identifying prognostic germline variants is challenging due to the limited effect size of germline variants, the large number of germline variants, and confounding clinical factors that may be associated with germline variants. Here we present a novel methodology for identifying prognostic germline variants and report two germline variants that we have found to be associated with survival in patients with lower grade glioma.
Materials and Methods
Glioma datasets
A total of 491 whole-exome–sequenced normal blood samples (WXS normal), 503 whole-exome–sequenced tumor samples (WXS tumor), and 501 RNA-sequenced tumor samples (RNA tumor) from TCGA lower grade glioma (4) patients available on the Cancer Genomics Cloud (CGC; ref. 15) platform were used as part of this analysis. The clinical information was downloaded directly from the TCGA data portal using the GenomicDataCommons (https://bioconductor.org/packages/release/bioc/html/GenomicDataCommons.html) R package available through Bioconductor. Additional molecular characteristics about these TCGA patients were acquired by downloading the supplement from Ceccarelli and colleagues (16). The raw sequencing data from the Chinese Glioma Genome Atlas patients was downloaded using accession number SRP027383 from the Sequence Read Archive. Clinical information for these patients was downloaded directly from the project's website (http://www.cgga.org.cn/).
Variant calling
Variant calling was performed on the TCGA lower grade glioma whole-exome–sequenced normal blood samples (WXS normal), whole-exome–sequenced tumor samples (WXS tumor), and RNA-sequenced tumor samples (RNA tumor) using VarDict (17) on CGC. The VarDict settings were set at default except for requiring mapping quality greater than 30, base quality greater than 25, a minimum of 3 variant reads, minimum allele frequency of 5%, and the removal of duplicate reads. We compiled a list of all of the unique variants and ran “samtools (18) depth” on all sequencing files requiring a mapping quality greater than 30. We determined the status of each variant in each patient from the three datasets (WXS normal sample, WXS tumor sample, and RNA tumor sample). The variant status at positions with fewer than ten reads for a given patient was changed to unknown. We used the WXS tumor samples to insert variant calls into the WXS normal samples at positions at which a variant status was listed as unknown in the WXS normal samples. If the variant status was still missing in a given patient, we then used the RNA tumor sample to insert variant calls into the combined WXS variant call set, allowing us to create the combined set of variant calls.
The same program parameters and approach were used to variant call and process the CGGA RNA sequencing dataset. All computation on the CGGA dataset was performed locally and not on CGC.
Quality control
We used Annovar (19) to determine the allele frequencies of the variants called by VarDict as listed in gnomAD (http://gnomad.broadinstitute.org/). We calculated the allele frequency of the variants in our study using the following formula:
The R package GGally (https://cran.r-project.org/web/packages/GGally/index.html) was used to calculate the correlation between the four variant call sets and to display their correlations with each other. Only variants with an allele frequency of greater than 5% in gnomAD and found in 15 or more of the TCGA lower grade glioma patients were tested for an association with survival by Cox regression.
Because we used the WXS tumor and RNA tumor samples to fill in missing variant calls, we evaluated whether somatic mutations were affecting the validity of our results. We first determined the percentage of variants called in the WXS tumor sample, which were somatic mutations. To do this, we downloaded the set of somatic mutations generated by the TCGA Research Network (20). We then calculated the number of somatic mutations called in each patient in this variant call set and divided that number by the total number of variants called in that patient's WXS normal sample. To assess whether somatic mutations were affecting the integrity of our results, we counted the number of times that a somatic mutation called by the TCGA Research Network overlapped with the set of germline variants that we were testing for an association with survival.
Because we used the RNA tumor sample to fill in missing variant calls, we evaluated whether RNA editing was having a significant impact on our analysis. To do this, we downloaded the set of over 2.5 million known RNA editing sites from a rigorously annotated database of RNA editing sites, RADAR (21). We counted the number of times that the germline variants that we were testing for an association with survival overlapped with any of the known 2.5 million RNA editing sites.
Principal component analysis
To calculate principal components that could separate patients on the basis of race, we used PLINK (22) to create a pruned set of germline variants to avoid bias from variants in linkage disequilibrium. Pruning was performed using a window size of 50 variants and a variance inflation factor of 2. These variants were used to calculate principal components using base R.
Cox regression and ROC curves
Lasso in the R package glmnet (23) was run on 17 covariates (Table 1). Information about patient age, gender, tumor location, grade, treatment site, and TP53 mutation status was acquired from the TCGA data portal, while data for patient somatic mutation count, percent aneuploidy, TERT expression, IDH mutation status, 1p/19q codeletion status, MGMT promoter methylation status, and chromosome 7 gain with chromosome 10 loss status was acquired from Ceccarelli and colleagues (16). The principal components were calculated as described above. Eleven of these 17 covariates were selected for inclusion in the final model for survival prediction. The R packages survival (24) and survminer (25) were used to run Cox regression and create Kaplan–Meier curves. For each minor allele, we tested with our model whether the minor allele was associated with a difference in survival outcomes with respect to the reference allele. False discovery rate correction was performed through Bonferroni correction.
List of variables that are known to be associated with differences in survival in lower grade glioma patients
ROC curves were created and evaluated using the survivalROC (https://cran.r-project.org/web/packages/survivalROC/survivalROC.pdf) and pROC (https://cran.r-project.org/web/packages/pROC/pROC.pdf) R packages. To test whether rs61757955 significantly improves the survival model consisting of the 11 covariates selected by Lasso, we compared the two ROC curves using the bootstrap method with 1,000 iterations. We also used this bootstrapping approach to determine whether ANKDD1a expression levels, GRB2 expression levels, rs61757955, and rs34988193 together improve the survival model with respect to the eleven covariates selected by Lasso.
RNA-sequencing data processing
We downloaded the HTSeq FPKM quantification files for each patient from the Genomic Data Commons data portal. We only used gene quantification files from primary tumor samples as part of this analysis. Replicate samples from a single patient were averaged.
Variant correlation to covariates and somatic mutations
To test for associations between the germline variants and genomic and histologic tumor characteristics, we divided patients based on their germline variant status. We used the Wilcoxon rank-sum test to test for significant differences in each of the continuous variables between patients with and without a given variant. We used Fisher exact test to test for differences in each of the discrete variables using a similar approach. Somatic mutation calls were downloaded from Ellrott and colleagues (20).
Gene set enrichment analysis
Gene set enrichment analysis (GSEA) of mRNA changes associated with rs61757955 and rs34988193 was performed by dividing the patients into two groups for each variant based on whether or not they had the reference allele at the position of the variant. For each germline variant, we calculated the log fold change for all genes expressed greater than one fragment per kilobase per million mapped reads (FPKM) between patients with the variant and without the variant. For each gene, fold change was calculated by dividing the median expression of the gene in patients with the variant by the median expression of the gene in patients without the variant. We used the log fold change to rank the genes from greatest log fold change to smallest log fold change. This file was used as input for GSEA (26).
Variant annotation
To identify deleterious mutations, we annotated all variants by combined annotation–dependent depletion (CADD) scores and only analyzed the variants predicted to be within the top 0.1% of all deleterious variants (CADD > 30; ref. 27). This led us to identify rs34988193 in ANKDD1a as a potentially deleterious variant predictive of survival. Because rs34988193 causes an amino acid change from positively charged lysine to negatively charged glutamate, we ran a BLASTp (https://blast.ncbi.nlm.nih.gov/Blast.cgi) search so that we could determine how many species have a protein homologous to ANKDD1a and how consistently the wild-type lysine residue was conserved. We identified homologous sequences in 22 other species. These sequences were aligned using ClustalW in MEGA (28). We also annotated this variant with its PhyloP score (29). Because the crystal structure for ANKDD1a was not available, we downloaded the predicted model for this protein from Modbase (https://modbase.compbio.ucsf.edu/modbase-cgi/index.cgi) and calculated the Gribskov score using prophecy on EMBOSS (30). We retrieved linked variants from Ensembl using the population of Utah residents with Northern and Western European ancestry which is demographically similar to the TCGA lower grade glioma patient population.
Results
Identification of high quality germline variants
Our variant calling pipeline is shown in Fig. 1. Briefly, we used the variant caller VarDict on Cancer Genomics Cloud to identify variants from whole exome sequencing (WXS) and RNA sequencing samples in about 500 lower grade glioma patients. In total, we found 4,453,701 unique variants. We used “samtools depth” to determine the sequencing depth at each of these variants for each patient and changed the variant status to “unknown” for patients with sequencing coverage less than 10 reads at a given position. We created a set of combined variant calls by using the WXS and RNA tumor samples to fill in unknown values in the whole-exome–sequenced normal samples that resulted from having a sequencing coverage of less than 10 reads at a given position. This approach increased our sample size and enabled us to include many more variants in our analysis than if we had solely used variant calls from the whole-exome–sequenced normal blood samples. Ultimately, this left us with four sets of variants – WXS normal, WXS tumor, RNA tumor, and a combined set that resulted from merging the other three variant call sets, giving preference to the WXS normal and then WXS tumor variant calls. We used the combined variant call set when testing variants for an association with survival. We only tested variants found in 15 or more lower grade glioma TCGA patients and listed in gnomAD as having an allele frequency of greater than 5%.
A flowchart describing the steps involved in identifying prognostic germline variants.
Tumor variant calls are not significantly affected by somatic mutations or RNA editing after filtering
Because we used sequencing data from the WXS tumor and RNA tumor samples to fill in missing calls in the WXS normal samples, we evaluated our variant calls for contributions from somatic mutations and RNA editing. We first showed that the majority of variant calls in the tumor sample are germline variant calls. To do this, we counted the number of somatic mutations called by the TCGA Research Network's analysis in each patient and divided that number by the number of variants that we called in the WXS normal sample (20). The median number of somatic mutations called per patient was 39. The median number of variants called in the WXS normal sample was 95,794. We therefore estimated that over 99.9% of variants called in the WXS tumor sample consisted of germline variants and that the percentage of somatic mutations in the WXS tumor sample across all patients was quite small (Supplementary Fig. S1). Because somatic mutations rarely occur at the same position, we suspected that the number of somatic mutations included in our study was extremely small because we limited our analysis to variants found in 15 or more of the patients with lower grade glioma found in gnomAD with an allele frequency of greater than 5%. Indeed, only one of the 196,022 variants that we tested overlapped with a somatic mutation. This somatic mutation occurred in only a single patient (Supplementary Table S1). Ultimately, we did not find any evidence to suggest that somatic mutations were impacting the quality of our analysis.
We next determined whether RNA editing was affecting our analysis by downloading the 2.5 million known RNA editing sites from the rigorously annotated RNA editing database, RADAR (21). Only 215 of the 196,022 variants that we tested were located at a position that overlapped with a known RNA editing site. We did not find any of these variants to be prognostic as part of our analysis. We therefore did not find any empirical evidence to suggest that somatic mutations or RNA editing impacted our findings (Supplementary Table S1).
Finally, we established that our four variant call sets (WXS normal, WXS tumor, RNA tumor, and combined) were concordant with each other by calculating the allele frequency of each variant called in the four sets and demonstrating a very strong correlation between all pairs of variants (r > 0.98 for all pairs, Supplementary Fig. S2). To further evaluate the quality of our variants calls, we calculated the frequency of each allele and compared it with the frequency of these alleles as listed in gnomAD. Our allele frequencies were well correlated with gnomAD (r > 0.93 for all four variant sets, Supplementary Table S2). As expected, the distribution of allele frequencies is negatively skewed as the majority of the identified variants are rare (Supplementary Fig. S2). We used the variants from the WXS normal samples to determine the principal components. As expected, these principal components effectively separate patients on the basis of reported race (Supplementary Fig. S3).
Identification of 271 prognostic germline variants that are independent of the clinical covariates
To identify clinically relevant germline variants, we restricted our analysis to variants found in at least 15 patients in the TCGA dataset and found in gnomAD with an allele frequency of greater than 5%. This restricted our analysis to 196,022 testable variants (Fig. 2A). To reduce the risk of identifying variants that are prognostic because they are confounded by other covariates known to be associated with survival, we used the machine learning algorithm Lasso to determine which of the 17 covariates should be controlled for in our Cox regression model. Lasso regression was useful in the screening of these 17 covariates because it penalizes models based on the number of coefficients, allowing for the elimination of less predictive coefficients from the model. The algorithm selected 10 covariates known to be associated with differences in survival in lower grade glioma (age, somatic mutation count, percent aneuploidy, histologic subtype of astrocytoma, tumor grade, treatment site, IDH mutation status, 1p/19q codeletion status, MGMT promoter methylation status, chromosome 7 gain/chromosome 10 loss status) along with the third principal component that we calculated (Table 1). Although the first two principal components are more effective in stratifying patients on the basis of race than the third principal component, the selection of the third principal component over the first two suggests that the third principal component contributes more information to the survival model than the first two principal components. This third principal component primarily separates African Americans from each other, suggesting that a subpopulation of African Americans experienced worse clinical outcomes in this dataset compared with other groups. We ran Cox regression on all 196,022 variants one at a time, controlling for these 11 covariates, to identify germline variants predictive of survival.
Prognostic germline variants in the TCGA dataset. A, Of the 4.4 million unique variants called in the TCGA dataset, we ran Cox regression on the 196,022 germline variants found in gnomAD with an allele frequency greater than 5% and found in 15 or more of the TCGA lower grade glioma patients. B–E, Similar to the 196,022 germline variants, the 271 prognostic variants are mostly found in protein-coding genes (B), are located in introns (C), and are single nucleotide polymorphisms (SNP; D). E, Most single-nucleotide polymorphisms cause transitions.
We identified 271 germline variants that are predictive of survival (P < 0.001) (Fig. 2A). As is the case with germline variants in general, the majority of these germline variants are found in protein-coding genes (Fig. 2B), are located in introns (Fig. 2C), and are single-nucleotide polymorphisms (Fig. 2D). Most single-nucleotide polymorphisms are transitions (Fig. 2E).
The germline variant rs61757955 in GRB2 is associated with poor prognosis
We identified two germline variants that are highly predictive of survival after false discovery rate correction (FDR < 0.10; Fig. 3A; Table 2A). rs61757955 results in a mutation in the 3′ UTR of Growth Factor Receptor Bound Protein 2 (GRB2) and is associated with a poor prognosis [P = 7.08E-10, HR = 20.4, Fig. 3B; Table 2A]. To determine whether rs61757955 enhances the survival model compared with the 11 clinical covariates alone, we calculated a risk score for each patient using a Cox regression model with rs61757955 and the other 11 covariates and a risk score using the 11 covariates alone. Using these risk scores, we determined the rate at which a patient would be correctly labeled as alive or dead at 7 years with a given false positive rate to create a ROC curve. The increased AUC suggests that rs61757955 enhances the survival model compared with the 11 clinical covariates alone (P = 0.0489, Fig. 3C). The allele frequency of rs61757955 is close to 0% according to the 1000 Genomes Project (31) in the Chinese population and, as expected, did not show up in the Chinese Glioma Genome Atlas. We also found rs28672782, a germline variant found in the intron of BRSK2, to be associated with a favorable prognosis, although the testable sample size for this variant was small and the maximum follow-up for patients with this variant was only 3 years. Therefore, we did not investigate this variant further (Supplementary Fig. S4; Table 2A).
rs61757955 is a highly prognostic germline variant identified in the TCGA dataset. A, Manhattan plot showing the P values resulting from testing each germline variant by Cox regression, controlling for the 11 variables in bolded in Table 1. Two variants passed the FDR threshold in the TCGA dataset. B, A Kaplan–Meier plot depicting the deleterious outcome associated with rs61757955, adjusting for the eleven covariates. C, ROC curve at 7 years. rs61757955 increases the area under the curve compared with the 11 covariates alone, suggesting that it improves the clinical model. D, Separation of patients on the basis of whether or not they have this germline variant to determine which genes are induced or repressed in patients with rs61757955. Subsequent gene set enrichment analysis reveals that patients with this germline variant exhibit upregulation of the genes involved with KRAS signaling.
Description of the prognostic germline variants identified in this study
To test whether rs61757955 in GRB2 is associated with an increased risk of other genomic abnormalities, we separated patients on the basis of this variant to see whether there was a difference in the incidence of the genomic or histologic variables (Table 3). We found this variant to be associated with an increased incidence of 1p/19q codeletions (P = 0.038). Because 1p/19q codeletions are frequently seen in Capicua transcriptional repressor (CIC)-mutated gliomas (32) and CIC aberrations are known to be a driver in lower grade glioma tumorigenesis (33), we tested whether there was a difference in the incidence of CIC mutations in patients with this variant. Thirty-eight percent of patients with this variant had CIC-mutated gliomas, whereas only 16% of patients without the variant had a CIC mutation (P = 0.0168; Table 3). Although the incidence of oligodendrogliomas was elevated in patients with the variant compared with patients without the variant, consistent with reports from the literature that 1p/19q codeletions and CIC mutations are enriched in oligodendrogliomas (32), this difference was not statistically significant (P = 0.475). Because rs61757955 is in a noncoding region, we also tested whether this variant is associated with differences in gene expression. We separated patients based on their variant status and calculated the log fold change of each gene between patients with the variant and patients without the variant. These data were used as the input for gene set enrichment analysis (GSEA). We found rs61757955 to be associated with increased KRAS signaling (FDR = 0.015; Fig. 3D).
The association between the germline variant rs61757955 and genomic and histological variables
Because we only have whole-exome sequencing and RNA-sequencing data from The Cancer Genome Atlas, we do not know whether the upregulation of genes in the KRAS signaling pathway and the increased incidence of CIC mutations and 1p/19q deletions are due to this variant or a linked variant in a regulatory region that we would be able to analyze with whole-genome–sequencing data. Therefore, we identified the four other variants that are genetically linked to rs61757955 in the European population, the population that is most similar to the TCGA lower grade glioma patient population (Supplementary Table S3). These variants did not pass the criteria to be included within the 196,022 testable variants that we had identified at the beginning of this study but could become useful in the future.
rs34988193 is a deleterious germline variant present in ANKDD1a associated with poor outcomes
To identify prognostic variants that are predicted to be deleterious due to effects on the encoded protein, we repeated our analysis, but restricted it to only variants with a CADD score greater than 30 and expression greater than one FPKM on average. Eighty-one variants met these criteria. These variants correspond to the top 0.1% of deleterious mutations as predicted by this scoring system. We found the germline variant rs34988193 in the tumor suppressor gene ANKDD1a to be associated with poor prognosis in the TCGA dataset (P = 0.001, HR = 1.73, FDR < 0.10, Fig. 4A and B; Table 2B). Because this variant is found in both the European and Asian populations, we were able to test whether this variant is also predictive of survival in the independent CGGA dataset. We found this variant to be predictive of survival in the CGGA dataset and we found the HR that we calculated in CGGA to be very similar to the HR calculated in the TCGA dataset (P = 0.0743; HR = 1.79, Fig. 4C; Table 2B). rs34988193 is not linked with any other variant in the European population. We did not find any enriched pathways after performing gene-set enrichment analysis and this variant was not associated with differences of any of the genomic or histologic variables (Supplementary Table S4).
rs34988193 is a prognostic variant predicted to be highly deleterious. A, A Manhattan plot with the P values resulting from testing each germline variant by Cox regression, controlling for the 11 covariates in bolded in Table 1. rs34988193 is prognostic (FDR < 0.10) in the TCGA when restricting the analysis to the top 0.1% most deleterious variants by CADD. B and C, Kaplan–Meier plots depicting the deleterious outcome associated with rs34988193 in the TCGA (C) and CGGA datasets (D), adjusted for the 11 covariates. D, A schematic showing that this variant is located in the ninth ankyrin repeat of ANKDD1a. The predicted protein structure of ANKDD1a reveals that this variant leads to an amino acid change from lysine to glutamate on the loop of an ankyrin repeat. E, Multiple sequence alignment of ANKDD1a in 22 species showing that lysine is conserved at this position in all of the species with this protein. F, ROC curves comparing the ability of two survival models to label patients as alive or dead after 7 years of follow up. The inclusion of rs61757955 variant status, rs34988193 variant status, GRB2 expression, and ANKDD1a expression significantly improves the survival prediction compared with the 11 covariates bolded in Table 1 alone (P = 0.0279).
ANKDD1a contains ten ankyrin repeat domains and one death-like domain. This variant causes a nonsynonymous mutation in the last codon of the ninth ankyrin repeat domain. The AAG to GAG codon change results in the incorporation of negatively charged glutamate instead of the wild-type positively charged lysine residue in the loop between ankyrin repeats nine and ten (Fig. 4D). This variant has a CADD score of 32 and is therefore predicted to be in the top 0.06% of deleterious mutations across the human genome. We performed a BLASTp search using the ANKDD1a protein sequence to identify homologous sequences in 22 other species. We aligned these sequences using ClustalW and found that this lysine residue is conserved in all 22 of these species (Fig. 4E). The PhyloP score at this position is 8.42, suggesting that evolution is occurring much more slowly than expected at this residue assuming no selection pressure. We determined the position-specific profile Gribskov's score for a lysine to glutamate amino acid change at this position using the multiple sequencing alignment from 23 species to be 15 to 3, suggesting that this variant is highly unfavorable.
Combined model predicts survival better than clinical covariates alone
As a result of this analysis, we found the germline variants rs61757955 in the 3′ UTR of GRB2 and rs34988193 in the protein-coding region of ANKDD1a to be predictive of survival in patients with lower grade glioma. We constructed a survival model consisting of the 11 clinical covariates, rs61757955, rs34988193, GRB2 expression, and ANKDD1a expression and generated a ROC curve by using this model to categorize patients as alive or dead after 7 years of follow up. This combined model is significantly better at predicting survival compared with the 11 clinical covariates alone (P = 0.0279; Fig. 4F).
Discussion
Up until this point, the identification of prognostic features in gliomas has been limited to clinical factors, somatic mutations, gene expression changes, and methylation pattern changes (6–8). Although many studies have commented on how germline variants could enable physicians to better individualize patient care by being able to better predict how a patient might respond to chemotherapeutic treatment (34–36), most large-scale studies have focused on identifying germline variants that predispose or protect an individual to a disease (13, 37). These studies have not focused on understanding how germline variants can be used to individualize patient care following diagnosis. Identifying prognostic germline variants is difficult due to the large number of germline variants, the limited effect of any single germline variant, and clinical factors that may confound the effect of germline variants. In this study, we have developed a novel method that can be used to identify prognostic germline variants and we have used that method to identify two variants that are predictive of survival in the TCGA dataset. The germline variant rs61757955 in GRB2 is not found in the Asian population and so could not be confirmed in an independent dataset. In contrast, the germline variant rs34988193 in ANKDD1a is found in both the European and Asian populations, and remarkably, was found to be prognostic with very similar HRs in both the TCGA and CGGA datasets.
Studies of germline variants using TCGA datasets typically solely utilize the WXS normal blood samples (13, 38). One major disadvantage to this approach is that it limits the analysis to genes within the capture regions of the whole-exome sequencing kits used by the study (4). In this study, we combined the information from both the whole-exome sequencing and RNA-sequencing datasets for a given patient to identify germline variants outside of the whole-exome sequencing capture region. Our approach had the added benefit of providing us with more information for a given variant for variants with low sequencing depth in the whole-exome sequencing datasets. We do not believe that this approach significantly affected the accuracy of our variant calls because the allele frequencies calculated from the RNA-sequencing dataset were well correlated with the allele frequencies from gnomAD and with the allele frequencies calculated from the whole-exome sequencing datasets. We showed that somatic mutations and RNA editing did not affect the integrity of our finding. Only one somatic mutation in a single patient overlapped with the 196,022 variants that we tested in our analysis and only 215 of the 196,022 variants that we tested overlapped with the 2.5 million known RNA editing sites. We did not find any of these variants to be predictive of survival. Instead, we feel that the increased sample size resulting from the additional sequencing coverage greatly outweighs any effect that somatic mutations or RNA editing had on our results.
We next needed to devise an approach to using these germline variants in a Cox regression model. We first had to decide how to deal with the absence of a variant in the variant call file. The variant could be absent because the patient was wild type for that allele or because the sequencing depth at that position was too low to make the variant call. We therefore determined the sequencing depth of each variant at each position so that we could exclude patients with low sequencing depths for the testing of specific variants. Testing a large number of variants increased the probability of a variant being significant solely because it was confounded with another significant variable. To avoid this issue, we tested each variant while controlling for 11 other covariates that we found to be predictive of survival. In this study, we found rs61757955 to be associated with differences in 1p/19q codeletion status. By including the 1p/19q codeletion as a covariate in our model, we were able to estimate the effect of rs61757955 independent from the 1p/19q codeletion status and the other 10 covariates.
GRB2 is a signal transduction adaptor protein that plays an oncogenic role in a variety of cancers (39–42). GRB2 plays an important role in the RAS/RAF/ERK pathway. Its SH2 domain binds the phosphotyrosine of activated growth factor receptor, while its two SH3 domains bind the guanine nucleotide exchange factor son of sevenless (SOS) protein, resulting in SOS recruitment to the plasma membrane and subsequent RAS activation. RAS binds and activates the kinase RAF, which phosphorylates the kinase MEK. MEK phosphorylates and activates extracellular signal-regulated kinase (ERK) which transmits the signal to transcription factors in the nucleus. This results in cell proliferation (43). We found the variant rs61757955 located in the 3′ UTR of GRB2 to be associated with poor prognosis in patients with glioma. Separating patients on the basis of this variant revealed that the KRAS signaling pathway is upregulated in patients with this variant. As described above, GRB2 plays a well-characterized role in this pathway (43). We also found this variant to be associated with an increased incidence of CIC mutations and 1p/19q codeletions. CIC is a known driver of lower grade glioma pathogenesis (33). Mutations in CIC are common in oligodendrogliomas and are associated with poor prognosis (4, 32). Although patients with rs61757955 variant exhibited an elevation in the incidence of oligodendrogliomas, which we expected given the increased incidence of CIC mutations and 1p/19q codeletions (32), this difference was not statistically significant. It is possible that this germline variant or the four other germline variants that it is linked with increase a patient's risk for oligodendrogliomas with the CIC mutation and 1p/19q codeletion.
In this study, we were only able to study variants in the whole-exome or RNA-sequencing data. Although it is possible that the 3′ UTR of GRB2 has regulatory activity or affects GRB2 protein translation efficiency, it is also possible that one of the variants that rs61757955 is linked to regulates the KRAS signaling pathway. None of the four linked variants are in the protein-coding sequence of GRB2 so that if they upregulate RAS activity, like the rs61757955, they likely do so by regulating the expression of GRB2. While recent large-scale sequencing studies have published patient whole-genome sequences (44), these data are not yet available for gliomas. We will be able to apply our approach to variants in regulatory regions in the future to specifically identify these prognostic variants when whole-genome sequencing data for gliomas is available. Our inability to further study this variant in the CGGA dataset due to this variant being rare in Asian populations is a limitation of this study, which could be addressed in the future with the availability of additional glioma-sequencing datasets. This result also suggests that the clinical usefulness of specific germline variants is dependent on the frequency of that germline variant in the population.
ANKDD1a is a tumor suppressor gene that has been shown to inhibit cell autophagy and induce apoptosis in glioblastoma multiforme (GBM). It directly interacts with and upregulates FIH1, resulting in inhibition of HIF1α activity and decreased HIF1α half-life. This induces apoptosis in GBM cell lines in hypoxic microenvironments. Hypermethylation of this gene is common in GBM and leads to decreased ANKDD1a expression and increased cell proliferation (45). We found the germline variant rs34988193, located at the end of the ninth of ten ankyrin repeat domains in this protein, to be associated with a poor prognosis in patients with lower grade glioma in both the TCGA and CGGA datasets. The HR independently calculated using the two datasets is remarkably similar. The wild-type lysine residue is conserved in all 22 species with a homolog to ANKDD1a and this position has a high PhyloP score. This variant is predicted to be within the top 0.06% of deleterious mutations in the human genome by CADD score (27), because it causes a change from a positively charged lysine residue to a negatively charged glutamic acid residue in the loop of this ankyrin repeat. Ankyrin repeats are common domains known for their involvement with protein–protein interactions (46, 47). Previous studies have suggested that mutations in the loops of ankyrin repeats may disrupt protein–protein interactions (48–50). The change from a positively to negatively charged amino acid resulting from the germline variant rs34988193 in the loop of ANKDD1a may disrupt ANKDD1a's protein interactions and could explain the poor prognosis associated with this variant seen in two independent datasets. Given the amino acid change, further studies involving rs34988193 in ANKDD1a could be directed toward experimentally determining whether or not this variant alters ANKDD1a's protein–protein interactions.
rs61757955 in GRB2 and rs34988193 in ANKDD1a could also be used to enhance predictions made by survival models clinically, as we found that these variants are significant predictors of prognosis even after controlling for 11 covariates. The prognostic effect of rs34988193 in ANKDD1a seems to be fairly reliable, as we found that this variant had a similar HR in both the TCGA and CGGA datasets. Our approach could be used in the future to identify sets of germline variants that together enhance the predictions made by survival models, although the current number of lower grade glioma sequencing samples is small relative to the large number of possible combinations of germline variants. Focused studies on particular sets of genes or pathways could potentially get around this low sample size problem by drastically limiting the number of variants studied. We believe that this study provides researchers with an effective approach to identifying biologically significant germline variants and provides clinicians with germline variants that could enhance currently existing survival models.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: A. Chatrath, A. Dutta
Development of methodology: A. Chatrath, P. Kumar, A. Ratan
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Chatrath
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Chatrath, M. Kiran, A. Ratan, A. Dutta
Writing, review, and/or revision of the manuscript: A. Chatrath, P. Kumar, A. Ratan, A. Dutta
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Chatrath, M. Kiran, P. Kumar, A. Dutta
Study supervision: P. Kumar, A. Ratan, A. Dutta
Acknowledgments
The authors thank Dr. Wei Min Chen at the University of Virginia for his statistical input on our Cox regression models. We thank Drs. Ana Damljanovic, Liz Williams, and Manisha Ray for their assistance with computation on the Cancer Genomics Cloud Platform. We thank dbGAP for providing us with access to The Cancer Genome Atlas data. We thank all of the patients and their families that participated in The Cancer Genome Atlas and Chinese Glioma Genome Atlas studies. We thank the Dutta lab members for the valuable feedback during the drafting of this manuscript. This work was supported by grants from the NIH R01 CA166054 and CA60499 and a Cancer Genomics Cloud Collaborative Support grant. The Seven Bridges Cancer Genomics Cloud has been funded in whole or in part with Federal funds from the National Cancer Institute, NIH, contract no. HHSN261201400008C and ID/IQ agreement no. 17 × 146 under contract no. HHSN261201500003I.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Footnotes
Note: Supplementary data for this article are available at Molecular Cancer Research Online (http://mcr.aacrjournals.org/).
- Received September 18, 2018.
- Revision received November 22, 2018.
- Accepted January 7, 2019.
- Published first January 16, 2019.
- ©2019 American Association for Cancer Research.
References
Volume 17, Issue 5
