Genotype-based gene signature of glioma risk

Abstract

Background.

Glioma accounts for 80% of malignant brain tumors, but its etiologic determinants remain elusive. Despite genetic susceptibility loci identified by genome-wide association study (GWAS), the agnostic approach leaves open the possibility that other susceptibility genes remain to be discovered. Here we conduct a gene-centric integrative GWAS (iGWAS) of glioma risk that combines transcriptomics and genetics.

Methods.

We synthesized a brain transcriptomics dataset (n = 354), a GWAS dataset (n = 4203), and an advanced glioma tumor transcriptomic dataset (n = 483) to conduct an iGWAS. Using the expression quantitative trait loci (eQTL) dataset, we built models to predict gene expression for the GWAS data, based on eQTL genotypes. With the predicted gene expression, iGWAS analyses were performed using a novel statistical method. Gene signature risk score was constructed using a penalized logistic regression model.

Results.

A total of 30527 transcripts were analyzed using the iGWAS approach. Four novel glioma susceptibility genes were identified with internal and external validation, including DRD5 (P = 3.0 × 10^–79), WDR1 (P = 8.4 × 10^–77), NOMO1 (P = 1.3 × 10^–25), and PDXDC1 (P = 8.3 × 10^–24). The genotype-predicted transcription pattern between cases and controls is consistent with that between tumor and its matched normal tissue. The genotype-based 4-gene signature improved the classification between glioma cases and controls based on age, gender, and population stratification, with area under the receiver operating characteristic curve increasing from 0.77 to 0.85 (P = 8.1 × 10^–23).

Conclusion.

A new genotype-based gene signature of glioma was identified using a novel iGWAS approach, which integrates multiplatform genomic data as well as different genetic association studies.

Gliomas account for 32% of all brain tumors and 80% of all malignant brain tumors.^1,2 Mortality is high in advanced gliomas (60% of gliomas), and patients with glioblastoma multiforme (GBM) have survival rates of 4.7% at 5 years.² Other than ionizing radiation, there are no established environmental risk factors for gliomas. Various environmental factors have been implicated in the etiology of gliomas, but most findings have been inconsistent, which suggests the importance of genetics and genomics in the susceptibility of this devastating disease. Several rare genetic syndromes, such as neurofibromatosis type I, have been associated with the glioma risk,³ and a positive family history is associated with a 2-fold elevated risk of glioma,^4,5 supporting the role for a genetic component to glioma. A number of genetic regions have been identified in genome-wide association studies (GWAS).^6–11 While recent findings from GWAS revealed important regions that are associated with the glioma risk, the agnostic approach inherent to GWAS and the strict adjustment for multiple comparisons leave open the possibility that other susceptibility regions remain to be discovered. Given the lack of well-established causes for glioma, understanding genetic susceptibility will provide new insights and opportunities for progress in unraveling the biological mechanisms behind this fatal cancer.

In GWAS, a large number of single nucleotide polymorphism (SNP) markers are tested across the genome. As multiple comparison adjustments are needed in GWAS, there has been a substantial interest in improving the statistical power of testing SNP effects by borrowing additional biological information. A major criticism of GWAS lies in its agnostic style¹²: no biological knowledge is encoded in the standard GWAS analyses. To address such limitations, SNP-set analyses have been advocated to integrate biological information into statistical analyses and to decrease the number of tests.^13,14 Analyses using SNP sets grouped by physical locations have shown a better performance than the standard single SNP analyses in re-analyzing the breast cancer GWAS dataset.¹⁴ SNPs can also be grouped into a set according to biological functions and have been utilized in studying skin cancer,¹⁵ bladder cancer,¹⁶ and lung cancer.¹⁷ By decreasing the number of tests and incorporating biological knowledge in the analysis, the SNP-set approach has provided a biologically relevant alternative to pursue genetic association analyses. However, this approach has not yet been applied to glioma.

Expression quantitative trait loci (eQTL) are the SNPs that are associated with gene expression. Several studies have incorporated eQTL data into GWAS using different approaches. Some studies used eQTL to prefilter or prioritize the SNPs: a GWAS of basal cell carcinoma that focused on eQTL SNPs was reported¹⁵; an osteoporosis GWAS used eQTL to re-prioritize the ranking from the result of genome-wide scans.¹⁸ Other studies focused on the overlap between GWAS and eQTL analyses: an asthma GWAS showed that the susceptibility loci at 17q21 were also eQTL for the ORMDL3 gene.¹⁹ Still others found that eQTL were enriched in the trait-associated SNPs; eQTL SNPs were also found enriched among common susceptibility loci of type 2 diabetes,²⁰ bipolar disorder,²¹ lymphocyte count,²² and the published GWAS SNPs in the online database collected by the National Human Genome Research Institute.²³ These studies have shown the promise of integrating eQTL into GWAS analyses.

Currently the existing approaches aim at the overlap of SNP–disease (GWAS) and SNP–expression (eQTL) associations using different strategies, but whether the association of SNPs with expression really contributes to the disease risk has not been directly examined in previous studies. Statistical methods have been developed to jointly analyze SNP and gene expression data provided that both data types are collected from the same individuals.^24–26 However, how to analytically integrate the genetic and transcriptomic data when the GWAS and the eQTL study are conducted in different subjects remains a challenge. Here we introduce a new analytic framework that utilizes a novel statistical methodology²⁷ to jointly analyze SNP and gene expression data by integrating a GWAS with an eQTL study.

Methods

Study Population

Three population datasets were included in this analysis: a brain eQTL dataset, a glioma GWAS dataset, and a GBM dataset from The Cancer Genome Atlas (TCGA). The brain eQTL study was conducted at the National Institute on Aging and consisted of genomic data obtained from fresh frozen tissue samples from the frontal lobe of cerebral cortex in 354 neurologically normal Caucasian subjects. Genome-wide genotyping was performed using Illumina Infinium HumanHap 550K, 610Q, or 660W BeadChips, and mRNA expression profiling was measured using Illumina HumanRef-8 Expression Beadchips. Additional details on this dataset can be found in Gibbs et al.²⁸

The glioma GWAS dataset was archived in the Database of Genotypes and Phenotypes (study accession phs000652.v1.p1). The genome-wide genotype data were collected on 556 glioma cases and 3647 controls using Illumina 550K, 610Q, or 660W BeadChip. We randomly divided the GWAS data into a Discovery Set consisting of 370 cases and 2430 controls and Validation Set consisting of 186 cases and 1217 controls. Missing genotypes were imputed using IMPUTE2 software (see Supplement Section I).

To validate the differentially expressed genes estimated for cases and controls in the GWAS data, we collected the array-based transcriptomic data from 473 GBM tumors or stage IV astrocytomas and 10 organ-specific normal control tissues from individuals without cancer who donated tissue for other reasons, both archived in TCGA.²⁹ The transcriptomic data were preprocessed level 3 data measured with University of North Carolina AgilentG4502A_07 array. The study is based on de-identified data that have been made publicly available, and thus does not involve human subjects. The demographics of the 3 datasets are summarized in Supplementary Table 1.

eQTL Analyses and Estimation of Gene Expression for GWAS Data

Three hundred fifty-four subjects of the eQTL dataset were analyzed to identify eQTL. We first defined cis-eQTL of a gene as SNPs locating within 0.5 Mb of the gene. We then constructed univariate eQTL models using linear regression models for each transcript expression value, regressing log2-transformed expression level on an SNP genotype under additive mode (0, 1, 2 as the number of the minor allele), adjusting for age and gender. For the cis-eQTL SNPs with P-value smaller than .05, we considered them as potential eQTL and built a multivariate eQTL model based on these eQTL SNPs. The distribution for the number of potential eQTL for 30527 transcripts is shown in Supplementary Figure 1 (median = 9).

We assumed a multiple linear regression model that log2-transformed expression level of each transcript $G$ was determined by covariates $X$ (1 [for the intercept], age, and gender), and SNPs $S=(S1, …, Sp)T$ from its corresponding eQTL set with size of p SNPs:

Due to the large number of eQTL SNPs, we employed a ridge estimator to estimate $α̂X$ and $α̂S$ that minimized $∑i(Gi−XiTαX−SiTαS)+λαS2$; the tuning parameter $λ$ was chosen from 10-fold cross-validations. We then predicted the unmeasured expression levels of probes for all 4349 samples in the GWAS study using the estimated model coefficients of SNPs, $μGi=SiTα̂S$. We obtained the predicted expression values of 30527 unique probes.

Integrative Genome-wide Association Study

We have developed a novel statistical method to integrate an eQTL study into a GWAS where the 2 studies were conducted in different study subjects.²⁷ We built eQTL models in (1) to estimate the expression value for each subject in the GWAS data and then combined the cis-eQTL SNPs mapped to the gene (or transcript) and its estimated expression level to assess their joint effect on the glioma risk by the following efficient testing procedure.

We first assumed a disease risk model, one transcript at a time: for subject $i$ ($i=1,…,n$; $n$ is the sample size of GWAS data), the glioma outcome $Yi$ ($Yi=1$ and $0$ for case and control, respectively) is determined by $p$ eQTL SNPs $Si=(S1i, …, Spi)T$ identified from the above, one mRNA expression of a transcript/gene $Gi$, their possible cross-product interactions as well as covariates $Xi$: 1 for the intercept, age, gender and 4 principal components for population stratification, through a logistic model:

We then obtain the marginal model for the glioma risk under the rare disease assumption that only depends on the eQTL SNPs $Si$ and the covariates $Xi$ by taking an integral with respect to gene expression $G$:

The null hypothesis of no genetic effect is specified as:

It has been shown that $H0:Δ=0↔H0:β=(βST, βG, βSGT)T=0$, provided $S$ are eQTL SNPs.²⁷ Because of the low power in conventional multivariate tests, we further assumed $βS$, $βG$, and $βSG$ followed arbitrary working distributions with mean zeros and variances $τS$, $τG$, and $τSG$, respectively. Thus the null hypothesis $H0:Δ=0$ became equivalent to:

We constructed the test statistic $Q$ as a weighted sum of noncentered scores $UτS$, $UτG$, and $UτSG$ for $τS$, $τG$, and $τSG$ under the null (2):

where $μ0=(μ01, …, μ0n)T$, $K=w1SST+w2μGμGT+w3CSGCSGT$, $S=(S1, …, Sn)T$, $μG=(μG1, …, μGn)T$, $CSG=(CSG1, …, CSGn)T$, and $CSGi=μGiSi$. We chose the weights $w1$, $w2$, and $w3$ to be the square root of the information for $UτS$, $UτG$, and $UτSG$, respectively. $μ0i$, the glioma risk for subject $i$ under the null (2), was estimated by $eβ̂XTXi(1+eβ̂XTXi)−1$, where $β̂X$ was estimated from the logistic model under the null. We calculated P-value for the test statistic $Q$ by comparing with its underlying distribution, a mixture of chi-square distributions.²⁷

We consider the following 3 hypotheses: (1) the SNPs-only model, $H0:τS=0$; (2) the model with SNPs and gene expression effects but no SNPs-by-expression interactions, $H0:τS=τG=0$; and (3) the model with SNPs and gene expression effects as well as their interactions, $H0:τS=τG=τSG=0$. Note the models (1) and (2) are parsimonious models nested within (3). In order to synthesize the information from the 3 candidate models, we conducted an omnibus test. Specifically, we calculated P-values for the 3 models $pS$, $pSG$, and $pSGC$ and used the smallest P-value $pmin=min(ps, pSG, pSGC)$ as the test statistic of the omnibus test (see Supplement Section II for further discussion). We were able to obtain the realization of the approximated distribution for the minimum P-value using a perturbation procedure, ${pmin(b) }={pmin(1) , …, pmin(B)}$, where $B$ is the number of perturbation.²⁷ By comparing $pmin$ with ${pmin(b) }$, we calculated the tail probability as the omnibus P-value. As the proposed method protects type I error for each transcript,²⁷ standard methods of adjusting multiple comparisons on P-values can be used (ie, Bonferroni’s correction, false discovery rate, etc). We used Bonferroni’s correction in the analysis of glioma risk. Sensitivity analyses by adjusting for population stratification in eQTL analyses are present in Supplement Section III (Supplementary Tables 4 and 5; Supplementary Figures 2 and 3).

Gene Signature

The transcripts identified in the iGWAS analysis were further combined to construct a risk prediction model that used transcript signature to predict glioma risk with the 2-stage discovery-validation process. We fit the following logistic regression model using the Discovery Set:

where $m$ is the number of genes or transcripts. Due to the large degrees of freedom and the correlation among $μG1$, . . . ,$μGm$, we introduced an $L2$ penalty into the maximum likelihood estimation to stabilize the estimator of $γ̂train=(γ̂1, …, γ̂m)T$ and $γ̂X, train$. The tuning parameter for the penalty was chosen to minimize mean squared error from 10-fold cross-validation. The risk score of the gene signature was derived as $ψi=μGSiTγ̂train$, where $μGSi=(μG1i,…,μGmi)T$. We compared the classification performance between the model with covariates and the m-transcript/gene risk score ($XiTγ̂X,train+ψi$) and that with only covariates ($XiTγ̂X,train$), using receiver operating characteristic (ROC) curve. To avoid overfitting in ROC, we calculated $ψi$ using the expression values $μGSi$ and the covariates $Xi$ in the Validation Set and $γ̂train$ and $γ̂X, train$ estimated from the Discovery Set (Figure 4A and C). We also swapped the two sets, ie, using the Validation Set to estimate $γ̂train$ and $γ̂X,train$ and calculate $ψi$ with the expression values $μGSi$ and the covariates $Xi$ in the Discovery Set (Figure 4B and D). Note that the risk score of the gene signature is the weighted sum of the genotype:

where $Sji$ is the genotype of eQTL for gene (or transcript) $j$ and $α̂Sj$ is the association of the eQTL SNPs with the gene (or transcript) expression.

Results

iGWAS

The analysis strategy is depicted in Figure 1. The glioma GWAS dataset was randomly divided into Discovery and Validation Sets. We conducted iGWAS to analyze 30527 transcripts using the Discovery Set. The distribution of the 30527 omnibus P-values is very close to a uniform distribution with a spike at the low end (Supplementary Figure 4). The quantile-quantile plot is shown in Supplementary Figure 5 with a genomic inflation factor of 1.061. Note that the 30527 transcripts may not be independent due to the fact that eQTL could be mapped to multiple transcripts and that multiple transcripts could be derived from the same gene. iGWAS omnibus P-values of the 30527 transcripts in the Discovery Set were presented in a Manhattan plot (Figure 2). Note that each dot represents a statistical significance level of a transcript where we jointly analyzed its expression value and the genotypes of its eQTL. The black dots are the 55 transcripts with P < .001 in both Discovery and Validation Sets.

Among the 30527 transcripts, 371 (1.30%), 397 (1.39%), 400 (1.40%), and 467 (1.63%) transcripts were significant at discovery P < .01 in tests for the SNP-only model, tests for SNP and gene expression main effect model, tests for interaction model, and omnibus tests, respectively; 98 (0.34%), 96 (0.34%), 101 (0.35%) and 107 (0.37%) transcripts were significant at discovery P < .001. It suggests that better statistical power was achieved by incorporating genotype-estimated gene expression, which is consistent with the numerical studies.²⁷

The 55 validated transcripts are summarized in Table 1 and Supplementary Table 2. Fifty-four of them had discovery P-values lower than the Bonferroni-adjusted genome-wide significance level, and all 55 combined P-values reached the genome-wide significance. Results of the SNP-only model, the SNP and gene expression main effect model, and the interaction model are summarized in Supplementary Table 3. The estimated gene expression provided advantage in detecting statistical significance in transcripts of OR7E85P, SLC2A9, DRD5, WDR1, VNN3, SNORD100, SNORA33, NOMO1, and PDXD1, but not the 40 transcripts of the ubiquitin specific peptidase 17-like family, NTAN1, STMN3, and LIME1, which had higher statistical significance in the SNP-only model. The most significant genes included DRD5 (dopamine receptor D5; P = 3.0 × 10^–79), WDR1 (WD repeat domain 1; P = 8.4 × 10^–77), SLC2A9 (solute carrier family 2 [facilitated glucose transporter] member 9; P = 1.4 × 10^–27), NOMO1 (NODAL modulator 1; P = 1.3 × 10^–25), and PDXDC1 (pyridoxal-dependent decarboxylase domain containing 1; P = 8.3 × 10^–24).

In addition to the iGWAS approach, we investigated candidate regions and previous GWAS loci using univariate genetic association analyses. For the 20 previously reported glioma (or GBM) susceptibility loci that can be mapped in our data, we confirmed 10 SNPs: rs2736100 (5p15.33; TERT), rs2853676 (5p15.33; TERT), rs2252586 (7p11.2; EGFR), rs4977756 (9p21.3; CDKN2BAS), rs1412829 (9p21.3; CDKN2BAS), rs1063192 (9p21.3; CDKN2A/B), rs2157719 (9p21.3; CDKN2A/B), rs3851634 (12q23.3; POLR3B), rs2297440 (20q13.33; RTEL1), rs6010620 (20q13.33; RTEL1, TNFRSF6B) (Supplementary Table 6). Based on our eQTL analyses, rs2297440 and rs6010620 were associated with the expression of STMN3 and LIME3.

As shown in the Manhattan plot, there are 4 hot spots of glioma susceptibility, including 4p16.1, 6q23.2, 16p33.11, and 20q13.33 (Figure 2). In 4p16.1, rs6824806 had highly significant genetic association with glioma risk, and the SNP also had very strong association with transcription levels of DRD5, WDR1, SLC2A9, OR7E85P, and ubiquitin specific pepetidase 17-like family. In 6q23.2, rs4458717 and r12200377 had significant association with SNORD100, SNORA33, and VNN3 expressions as well as the glioma risk. In 16p33.11, rs11075260 was associated with PDXDC1, NOMO1, NTAN1 expressions, and the glioma risk. In 20q13.33, rs2297440, rs6010620, rs6089953 had strong association with the glioma risk and moderate association with STMN3 and LIME3 transcription level. We note that 20q13.33 is also the region that harbors significant SNPs (rs2297440 and rs6010620) in previous GWAS.^6,9

Genotype-based Gene Signature

Hierarchical clustering was performed based on Euclidean distance among the estimated expression values of the 55 validated transcripts (Figure 3). Glioma cases are enriched in the red cluster, compared with the other clusters (20.1% vs 8.1%, P = 1.2 × 10^–29). With the 55 transcripts, we developed a gene signature with use of $L2$ penalized logistic regression. Compared with the model with only covariates: age, gender, and 4 principal components of population stratification, the model with the additional 55 transcripts had better discrimination between glioma cases and controls with area under the ROC curve (AUC) increasing from 0.76 to 0.83 (P = 1.3 × 10^–5) in the Validation Set with models built by the Discovery Set (Figure 4A), from 0.77 to 0.88 (P = 1.1 × 10^–35) in the Discovery Set with models built by the Validation Set (Figure 4B), and from 0.77 to 0.88 (P = 1.0 × 10^–35) with model building and testing using both sets combined. ROC curves for a single transcript of all 55 transcripts are shown in Supplementary Figure 6.

Among the 55 transcripts, 8 (DRD5, WDR1, SLC2A9, VNN3, NOMO1, PDXDC1, STMN3, and LIME1) can be identified in TCGA gene expression data of GBM. We investigated the estimated gene expression by comparing differential expression patterns between cases and controls with those in TCGA GBM tumors versus organ-specific control tissues. Four genes (DRD5, WDR1, NOMO1, and PDXDC1) were validated: expression of DRD5, NOMO1, and PDXDC1 was significantly lower and that of WDR1 was significantly higher in GBM tumor tissue than normal brain tissue, which is consistent with findings based on the estimated gene expression for the glioma GWAS data (Figure 5). With the 4 genes further validated by TCGA data, we developed a parsimonious 4-gene signature. Its performance is similar to the 55-transcript signature, with AUC increasing from 0.76 (the model with only covariates) to 0.82 (the model with covariates and the 4 genes) (P = 2.3 × 10^–4) in the Validation Set (Figure 4C), from 0.77 to 0.86 (P = 1.4 × 10^–27) in the Discovery Set (Figure 4C), and from 0.77 to 0.85 (P = 8.1 × 10^–23) in both sets combined. The other 4 genes that were not validated by TCGA data were SLC2A9, VNN3, STMN3, and LIME1. Interestingly, the signature constructed by these 4 genes had only modest increase in AUC (from 0.76 to 0.81 in the Discovery Set; from 0.76 to 0.79 in the Validation Set) (Supplementary Figure 7). The genotype-estimated expressions of the 55 transcripts between glioma cases and controls are shown in Supplementary Figure 8.

Discussion

Here we present a new approach to conduct a genome-wide association study, which identifies novel susceptibility genes of glioma risk: DRD5, NOMO1, PDXDC1, and WDR1. Gene-centric GWAS analysis decreases the number of tests and has become popular. Genotype-based methods such as the SNP-set Sequence Kernel Association Test (SKAT) analyze all genetic variants within a gene and have shown better power than single-SNP analyses.^14,30 A special case of our method focusing on an SNP-only model can be construed as SKAT analyses with the linear kernel. Similar to our approach, an eQTL-based method PrediXcan proposed to utilize reference data to impute transcriptome in GWAS.³¹ The difference is that PrediXcan focused only on the association of imputed gene expression with phenotype, and we jointly analyze the genotype data, the imputed gene expression data, as well as the genotype-by-transcription interaction using a score-type variance component test.²⁷ It is critical to include the genetic effect, which captures the trans-eQTL effect through other genes or effects of other transcriptional regulation and biological pathways not via the imputed transcription. It has also been shown that analysis with single platform loses statistical power compared with the proposed joint analyses.²⁷ By integrating both genetic and transcriptomic data, we present a unified framework where SKAT and PrediXcan are 2 special cases, and our approach further provides omnibus tests to synthesize the optimal information.

A challenge of jointly analyzing SNP and gene expression data in genome-wide association studies is its difficulty in gaining access to the target tissue. In most GWAS, investigators collect peripheral blood samples and genotype DNA extracted from blood cells. While mRNA can also be extracted and profiled from blood cells to study immune-related diseases,¹⁹ it is subject to serious limitations in interpretation when generalizing to non-immune related diseases such as glioma, since brain tissue is the target tissue, not blood cells. For glioma, it is not impossible to obtain target tissue from cases; however, it is often difficult and unethical if not impossible to obtain target tissue from controls. To integrate genomic data collected from different genetic association studies, we have developed a statistical method²⁷ under the framework of causal mediation modeling.^32,33 Utilizing the newly developed algorithm, our proposed iGWAS approach provides an analytic paradigm to exploit the information from an external eQTL study to infer the missing transcriptomic profile for GWAS subjects. We note that our eQTL analyses focused on the expression profile from only the cerebral tissue that is only limited to the frontal lobe, which may attenuate the signals or even diminish the likelihood of identifying signals. The information of anatomical sites can be easily incorporated in our analytic framework by developing lobe-specific eQTL models for prediction, which may increase the power of detecting susceptibility genes.

In the example illustrated here, we identified 4 new genes that have not been linked to glioma before: DRD5, NOMO1, PDXDC1, and WDR1. DRD5 encodes the D5 subtype of the dopamine receptor, a G-protein coupled receptor which stimulates adenylyl cyclase.³⁴ Consistent with our findings that low expression was associated with glioma risk (Figure 5A and B), the Human Protein Atlas shows that DRD5 protein has high expression in cerebral cortex but is not detectable in glioma tumor tissue; the expression in cerebral cortex is mostly observed in neuropil but not in glial cells, endothelial cells, or neuronal cells.^35,36DRD5 has been associated with neurologic disorders such as Parkinson’s disease,³⁷ multiple sclerosis,³⁸ schizophrenia,³⁹ and attention-deficit hyperactivity disorder,⁴⁰ and is also an FDA-approved drug target.⁴¹ However, the literature is very limited on its role in carcinogenesis, which deserves more research.

NOMO1, also known as PM5, is one of the 3 highly similar genes located in a duplication region on p arm of chromosome 16⁴²; the 3 genes encode proteins that may have the same function, and one protein has been identified as part of a complex involved in the nodal signaling pathway during development.^34,43 NOMO1 protein is highly expressed in normal cerebral cortex, particularly in glial cells and neuronal cells, but not in endothelial cells or neuropil; in contrast, its expression in glioma tissue ranges from medium to nondetectable.^35,36 The protein expression pattern is consistent with our findings in the GBM data from TCGA (Figure 5F) as well as in the eQTL-estimated expression (Figure 5E). NOMO1 was found overexpressed in the cutaneous T-cell lymphoma cell line compared with normal peripheral blood monocytes,⁴⁴ but little is known about its role in glioma. Discussion of WDR1 and PDXDC1 is provided in the Supplementary material.

We note that the results from TCGA should be interpreted with caution because the set from TCGA was derived from tumor DNA, and GWAS and eQTL sets were derived from normal DNA. Because gene expression is profiled after tumors have developed, such differential expression could be from either causal genes that induce carcinogenesis or reactive genes with their expression altered secondary to cancer development.⁴⁵ The gene expression of the tumor consists of 2 components: the causal gene expression and reactive gene expression. Since both GWAS and eQTL sets were derived from the nontumor DNA, we hypothesize that the predicted transcriptome is more likely to represent the causal gene expression than the reactive one. The 4 genes, SLC2A9, LIME1, STMN3, and VNN3 with estimated expression not validated by TCGA data could be due to the difference of causal (in iGWAS data) and reactive (in TCGA data) expression, or simply false positives. Moreover, the validity of the estimated transcriptomic data for the GWAS subjects by the eQTL data relies on the assumption that the eQTL subjects are representative of the GWAS subjects, conditioning on age, gender, and ethnicity/population stratification. Exploratory analyses using low-grade glioma tumors were also conducted to validate DRD5 and WDR1 (Supplement Section IV). We expect that the set from TCGA tends to provide a more conservative validation, ie, we may miss causal signals due to its being intertwined with the reactive genes.

In conclusion, we identified 4 novel susceptibility genes of glioma: DRD5, NOMO1, PDXDC1, and WDR1, and constructed a genotype-based gene signature of glioma risk using integrative genomic analytics, iGWAS. The iGWAS approach integrates multiplatform genomic data, that is, SNP and gene expression data as well as different genetic association studies (ie, an eQTL study and a GWAS). The iGWAS has advantages of identifying biologically plausible results and improving statistical power, which make it a promising strategy to revisit existing GWAS data identifying new disease susceptibility genes that are transcriptionally functional.

Supplementary Material

Supplementary material is available at Neuro-Oncology online.

Funding

National Cancer Institute grant 5R03CA182937-02 (Integrative Modeling of Gene Expression into GWAS of Glioma, to Y-T.H., D.S.M.)

Conflict of interest statement. The authors report they have none to declare.

Acknowledgments

We thank the US National Institutes of Health/National Cancer Institute grant 5R03CA182937-02 (Integrative Modeling of Gene Expression into GWAS of Glioma, to Y-T.H., D.S.M.) for the support of the analysis project we present in this paper, and thank Dr J. Raphael Gibbs for sharing the brain eQTL data and providing the data preprocessing information.

References

Author notes