Pathway based analysis of genotypes in relation to alcohol dependence.

We introduce a method for detecting variants in several genes of related function with small effect on a phenotype of interest. Our method uses logistic regression to test whether multiple alleles within a functional set have significantly higher than expected predictive value, even though none individually may have strong individual effects. We illustrate this method by testing seven gene sets (including 48 genes), from a study with 1350 single nucleotide polymorphisms in 130 addiction candidate genes studied in a sample of 575 alcohol dependence (AD) cases and 530 controls. We conclude that AD is related to variation in genes participating in Glutamate and γ-amino butyric acid signaling, as has been reported elsewhere, and in stress response pathways, but not with genes in several other systems implicated in other drugs of abuse.

Introduction

Modern efforts to identify the genetic risk factors for complex psychiatric and drug abuse disorders have focused mostly on genomic regions (in linkage analysis) or individual sequence variants (or associated haplotypes) in or near specific genes (in candidate gene association studies). A few studies 1–5 have examined epistatic models of disease etiology; these studies have, typically examined only two markers, loci or linked regions at a time. Recently genome-wide association studies (GWAS) for psychiatric disorders have appeared 6–9.

In this article, we apply a different and systematic approach towards uncovering the genetic basis of complex diseases – one based on the analysis of sets of genes whose protein products collaborate in specific functions. This may be thought of as in-between a candidate gene study and a GWAS. For the purposes of this paper, we define a functional gene set as one whose protein products are involved in a closely coordinated biological function. While network approaches are common in the analysis of expression arrays, we are aware of only three prior steps in this direction in linkage or association studies in relation to complex disease 10–12. These papers have used methods devised for continuous gene expression data and adapted them to discrete genotype data. By contrast our approach is grounded in the methodology of statistical genetics extended in a natural manner to gene networks.

We conceptualize the gene set approach as follows. Consider a particular set of genes that produce proteins that act in a pathway to produce a given metabolite, or that convey a signal from the cell surface to the cell nucleus, or that come together to form a multimeric receptor. Functional variants in any one of the genes might produce abnormal proteins, many of which would impair the function accomplished by that set. Many variants in that gene set would result in similar perturbations in biological function, which could then affect propensity to a complex disease. If one variant in one of these genes conferred substantially increased disease risk, then that variant would be relatively easy to detect by traditional association methods. However natural selection tends to remove such genes, but is less efficient at removing variants which increase risk only slightly. If each of several variants of genes within that set confer modest increases in disease risk, those variants might be quite difficult to detect when examined one at a time. However greater statistical power might be achieved toward detecting a disease association if all these variants were analyzed jointly as an statistical unit.

Recent comprehensive studies of genetic contributions to disease in other fields have also found many fewer candidates than would be needed to explain the documented heritability of the disorders 13. Current discussions emphasize the role of multiple genes and posit a ‘long tail’ – many genes with very small effects on risk could account for a large fraction of the variance 14, 15. Such loci might be extremely difficult to detect by classical genetic methods that focus on a single locus at a time.

In this paper, we attempt to provide a “proof of principle” for our method by applying it to genotype data from 370 SNPs in 48 putative candidate genes obtained from a sample of individuals with alcohol dependence (AD) and population controls. A number of prior studies have examined risk loci for AD and their results may serve as a test of our approach. An unresolved issue is whether common complex disorders reflect the effects of a few common variants of several genes or the dramatic effects of many rare variants. Here we propose a method to address, from a gene set perspective, the common disease/common variant (CDCV) model.

Materials and Methods

Subjects and phenotype measurement

Participants in the Irish Affected Sibpair Study of Alcohol Dependence (IASPSAD) were recruited in Ireland and Northern Ireland between 1998 and 2002. Further details of the study design, sample ascertainment, and clinical characteristics of this sample are described elsewhere16. In brief, ascertainment of probands was mainly conducted in community alcoholism treatment facilities and public and private hospitals. Probands were eligible for study inclusion if they met the current DSM-IV criteria (American Psychiatric Association, 1994) for AD and if all four grandparents had been born in Ireland, Northern Ireland, Scotland, Wales or England. After a prospective family was identified through a proband, parents and potentially affected siblings whom the probands provided permission to contact were recruited.

Probands, siblings and parents were interviewed by clinically-trained research interviewers, most of whom had extensive clinical experience with alcoholism. The assessment included demographic characteristics, lifetime history of AD and other comorbid conditions, alcohol-related traits, personality features, and clinical records. The DSM-IV AD diagnosis was assessed in probands and siblings using the modified Semi-Structured Assessment of the Genetics of Alcoholism (SSAGA) version 11 17 interview to reduce assessment time.

All participants provided informed consent. There were 1238 individuals who completed the SSAGA interview and met DSM-IV AD diagnosis, including 591 probands, 620 affected siblings, and 27 affected individuals from 10 complex families. Controls were recruited in the Northern Ireland from volunteers donating at the Northern Ireland Blood Transfusion Service (N=554) and in the Republic from the Garda Siochana (the national police force, N=38) and the Forsa Cosanta Aituil (the army reserve, N=34). Controls were screened and their samples excluded if they reported a history of heavy drinking or problem alcohol use. In the present case-control study design, we included 575 independent AD cases (399 probands and 176 sibs) that were drawn one each from the IASPSAD families and 530 controls, with high yield of high quality DNA for genotyping.

SNP selection for the Addiction Array

For this study, we tested 48 genes which had previously been tested individually on a custom addiction array. The array, developed by Dr. David Goldman, tested a total of 130 addiction candidate genes, selected from multiple functional systems and implicated in substance dependence phenotypes (primarily alcohol, cocaine and opiates). For the array genes, a genomic region containing 5 kb upstream and 1 kb downstream of each candidate gene was retrieved from NCBI (Human Genome Build 35.1) and genotype data from the African population, which is the most diverse, were obtained from HapMap Phase I Rel18 to re-construct haplotypes for each gene using SNPHAP (http://www.gene.cimr.cam.ac.uk/clayton/software/snphap.txt). A double classification tree search algorithm was applied to select minimum index SNPs representing maximum haplotype information and with frequency >0.6% for each gene. Probable functional SNPs (non-synonymous, splice site and putative functional SNPs from the literature) were forced in during the selection process. The performance of the initially selected SNP set was validated by the manufacturer and replacements made where necessary. In total, there were 1536 SNPs on this Illumina array, including 1350 SNPs from the 130 candidate genes, and 186 ancestry informative markers (AIMs). Our Irish population showed more linkage disequilibrium than the African population for which the chip was designed: the median r2 between neighboring variants was 0.21; 145 pairs of neighbors had r2 > 0.8. By construction of the algorithm SNPs in high LD could not both count toward the significance of a gene set.

Genotyping using an Illumina custom SNP array

Genotyping for the addiction array was performed using the Illumina GoldenGate genotyping protocols on 96-well format SentrixR arrays. All pre-PCR processing was performed using a TECAN liquid handling robot and arrays were imaged using an Illumina Beadstation GX500. The data were analyzed using GenCall v6.2.0.4 and GTS Reports software v5.1.2.0 (Illumina). Full details of array genotyping and post-experiment QC are described elsewhere 18. In brief, a total of 1105 unique samples were genotyped. After screening for poor performing SNPs based on technical metrics from the GX500 and assessment of Hardy-Weinberg equilibrium, the data available for this sample included 1286/1350 (95.3%) of gene SNPs and 176/186 (94.6%) AIM markers. Markers which passed QC metrics were then assessed for Hardy-Weinberg equilibrium, and those with p < 0.001 were discarded. Those SNPs or individuals in which more than 10% values were missing were dropped from further analysis. Any SNPs with minor allele frequency less than 5% were also discarded.

Identification of signaling gene sets

In this study, we focused on a subset of signaling systems which were well represented in the available data. We divided up the 130 assayed genes into families based on the signaling systems in which they participate. In most cases these were sets of receptors for a particular signaling molecule, which annotations were provided with the chip; the genes for CRH were selected by DK who has specialized expertise. This process yielded seven gene sets which had at least four member genes; several other signaling systems (e.g. opioids) had fewer than four genes represented on the chip; these were not considered further for a systems-analysis approach. We concentrated on those seven, which included receptors and other molecules involved in neurotransmission via norepinephrine, acetylcholine, dopamine, CRH, gamma-amino butyric acid (GABA), glutamate, and serotonin. The genes in each group were as comprehensive as possible limited only by the coverage on the chip. The genes in the gene sets are shown in Table 1.

Table 1

Genes considered within gene sets

Norepinephrine system	ADRA1A, ADRA2A, ADRA2B, ADRA2C, ADRB2
Cholinergic system	CHRM1, CHRM2, CHRM3, CHRM4, CHRM5
Dopaminergic system	DRD1, DRD2, DRD3, DRD4, DRD5, DBH, DDC, TH
GABAergic system	GABRA2, GABRA3, GABRA4, GABRA6, GABRB1, GABRB2, GABRB3, GABRD, GABRE, GABRG2, GABRG3, GABRQ, SLC6A7
Glutaminergic system	GRIK1, GRIN1, GRIN2A, GRIN2B, GRIN2C, GRM1
Serotonergic system	HTR1A, HTR1B, HTR2A, HTR2B, HTR2C, HTR3A, HTR3B
CRH system	CRH, CRHBP, CRHR1, CRHR2

Statistical approach

Model

Our model is directed toward the Common Disease Common Variant hypothesis. We posit that each variant contributes independently and additively to risk. This can be simply expressed in a logistic regression model: the logit of relative risk is predicted by a weighted sum of counts of risk-bearing alleles as in equation (1): where n represents the counts of allele i and β reflects the contribution of allele i to the phenotype. We are anticipating many small contributions of common alleles to AD risk and therefore anticipate that many of the true coefficients β will be non-zero but not large. Therefore when we fit model (1) to our data, we expect to see that many of the fitted estimates b, of β, will be significantly different from 0, i.e. more fitted coefficients b will be associated with large increases in deviance explained in model (1), and therefore small p-values, than we would see in similarly powered studies if there no underlying effects (if all β were 0). In order to test whether there are more non-zero coefficients than we would expect by chance, we permute the phenotype labels to generate many data sets with the same LD structure among the genotypes, but no consistent relation to the AD phenotype, and ask how unusual is it to find as many significant estimated coefficients b as in fact occur with the true phenotype labels.

SNP selection

We expect that most typed variants will not be even weakly associated with AD risk, and statistical testing would be weakened by including all variants. We therefore selected subsets of SNPs across gene sets in two ways. First we followed standard practice for building a multivariate linear model based on increase in log-likelihood. We set a relatively modest threshold of a log likelihood increase of 2.5 for inclusion in a model. Second we implemented a faster procedure by selecting a subset of alleles that were potentially individually associated with AD using a liberal threshold for inclusion of p < 0.20. Both these inclusion criteria are fairly weak, and statistical significance cannot be assessed by conventional means. For both procedures we prevented from entering the model those variants in LD (r2 > .25) with variants already in the model.

Statistical significance

In order to properly assess significance of the AD risk contributions of selected SNPs, therefore, we generated a null distribution by repeating the whole procedure, including SNP selection, 1,000 times with randomized phenotype data (for the LR inclusion criterion) or 10,000 times (for the faster p-value threshold criterion); for each permutation we recorded the total deviance. We then compared the total deviance increase for the SNPs chosen using the true outcomes to the distribution of total deviance increases for SNPs chosen using the identical procedure for randomized outcomes; the p-value associated with each gene set was the proportion of permutation fits whose total explained deviance exceeded that of the true fit. We found that the sets of SNPs chosen and the significance of gene sets differ only modestly depending on which criterion for SNP inclusion is employed or on which log likelihood threshold is used for the likelihood criterion. We report results for the likelihood criterion, since this is more traditional; the results based on the other inclusion criterion were slightly more significant (not shown).

Results

The individual genes selected for each of the 7 gene sets are shown in Table 1 and main results of our analyses are shown in Table 2.

Table 2

Summary of gene sets and significance of permutation testing

	ACh	CRH	Dopamine	GABA	Glutamate	NE	Serotonin
Genes	5	4	8	13	6	5	7
SNPs assessed	40	22	58	121	65	20	44
SNPs selected	2	3	4	13	7	3	3
Nominal p-value	4.2 ×10−2	1.1 ×10−2	2 ×10−3	1.4 ×10−8	3.2 ×10−5	6.7 ×10−3	1.3 ×10−2
Permutation p-value	0.78	0.01	0.32	0.04	0.05	0.02	0.59
Q-value	0.78	0.07	.45	0.09	0.09	0.07	0.69

Using the standard Benjamini-Hochberg procedure for FDR, at an expected proportion of 10% of false positives, we identified the following four systems as potential contributors to genetic risk for AD: norepinephrine, glutamate, GABA, and CRH. The 27 SNPs in these systems together accounted for 6.8 % of the total phenotype variance, although this is likely an overestimate due to selection bias.

If we had done a standard chip-wide association analysis on this data, none of the individual variants in this study would be significantly associated with AD after multiple comparisons corrections. In fact, the smallest Q-value would be 0.7 using the standard Benjamini-Hochberg procedure. However among our seven gene sets, four gene sets (norepinephrine, glutamate, GABA, and CRH) achieved a Q-value < 0.1.

In order to examine the distribution of effect sizes within the positive gene sets, we plotted the odds-ratios for all of the SNPs selected by the algorithm in the four systems selected by our procedure (Fig 1). Modest to moderate association signals were distributed widely within both gene sets and were not confined to a few markers. To give some sense of the statistical significance of the odds-ratios (ORs) associated with these SNPs, for associations with a randomized phenotype, only 1 in 100 SNPs would have an OR ≥ 1.25 We examined 165 SNPs in the GABA system and 90 in the Glutamate system. Thus, one would expect one or two of the ORs ratios from GABA receptors SNPs to exceed 1.25, whereas we observe nine. For glutamate receptor SNPs, we would expect 1 and we observe 4. For the NE system we would expect to see 0 or 1 SNP whose OR exceeds 1.25; we see three. For the CRH system we expect on average to see one SNP whose OR exceeds 1.2 one half of the time; in fact we see two.

Figure 1

Upper Left Effect sizes (as odds-ratios) for SNPs selected in the GABA gene set. Three SNPs for GABA occur in the same gene (GABRG3), but the markers are not in LD (r2 < 0.25 by construction). Upper Right. Odds ratios for SNPs in the glutamate receptor gene set. Lower Left Odds ratios for selected SNP’s in the Norepinephrine gene set. Lower Right. Odds ratios for SNP’s selected in the CRH gene set.

Discussion

The goal of this report was to test a new method for the detection of small effect risk variants for complex neuropsychiatric disorders. We proposed 7 gene sets based on a priori biological hypotheses and tested these sets for their impact on risk for AD. Using the standard Benjamini-Hochberg FDR procedure, we identified four gene sets significantly associated with AD with at an FDR of 10%, indicating a low likelihood of their being false discoveries. In contrast, no individual SNP in any of the selected gene sets would have been selected in a chip-wide study at an FDR less than 40%. If we had hypothesized a priori that our candidate genes would belong to one of the seven groups tested and restricted our single locus tests to those SNPs only, then we would have had only one SNP significant at an FDR of 0.2, (which would also be a Bonferroni-corrected p-value of 0.2), but none at an FDR of 0.15.

Interestingly, the distribution of the association signals in the two largest sets was among several genes rather than being concentrated in one gene. Such a wide distribution of effects would be predicted by the original motivation for this gene set approach.

The validity of our findings is also supported by the work of many previous researchers who have identified the contribution of GABA and glutamate receptors to AD 19. Some of the genes contained in the gene sets tested here have been previously tested as candidate genes albeit with inconsistent results. Some of the others have yet to reveal association with AD; however since the effect sizes are small, they would not likely have yielded association when tested individually, even if the association is real. For example, previous studies have shown that a number of individual genes for components of the GABAA receptor are associated with AD. Investigations in the COGA sample showed that variations in GABRA2 were highly associated with AD as well as the beta frequency of the electroencephalogram 20. Another study 21 produced evidence showing gene-environment correlation and interaction in the impact of GABRA2 variants and risk for AD. Other studies have reported association between AD and GABRA6 22–24, GABRB3 and severe alcoholism 25 and Song et al 24 found evidence of association with GABRA5 and GABRB3. The most frequent combination of GABAA subunits is α1, β1, and δ 26. Of these three common subunits, only β1 and δ were sampled on the array, but when SNPs for these the genes encoding these subunits were considered together as a group they also reached a significance level of 0.025.

There is also increasing evidence that glutamatergic neurotransmission is involved in alcohol tolerance and withdrawal through its role in synaptic plasticity 27, 28. A limited number of studies testing NMDA receptor subunit genes in human samples have reported association with AD and related phenotypes 29–33. Genetic association results for genes coding subunits of the NMDA, kainate and AMPA receptors have produced mixed results 30 despite the substantial functional evidence implicating this system in AD. Results from the current study suggest that other glutamatergic genes also represent putative risk factors for AD.

Although genes in the Norepinephrine and CRH systems have not been tested for AD to the extent that GABA and glutamate have, there is evidence to suggest that the Norepinephrine and CRH systems have a functional role in the etiology of AD. Previous findings have shown that stress system disturbances are observed in AD and are involved in alcohol-induced neurotoxicity and cognitive deficits 34. The Hypothalamic-pituitary-adrenal (HPA) axis activation influences reward pathways 35 and GABAA receptor function involved in ethanol sensitivity 36. Evidence for genetic factors is supported by family studies showing altered HPA response to ethanol 37, 38, animal studies linking a CRH haplotype in alcohol consumption 39, and clinical studies associating CRHR1 variants with drinking patterns 40, 41. Findings of HPA axis disturbance in drinking or abstinent AD individuals, however, are inconsistent and difficult to interpret, with both heightened and reduced activity reported 34, 42. Moreover, acute and chronic alcohol use differ in their effects on HPA axis activity, with disturbed activity in AD potentially due in part to the impact of chronic alcohol consumption on brain regions involved in negative feedback 43. HPA disturbance in AD may also be due to the high comorbidity of AD with major depression (MD). Support for an etiologic role for altered stress system activity in the pathophysiology of MD is growing, with some evidence that genetic variation may be involved 44–46. Epidemiological research suggests MD may be the primary pathology in some AD patients 47. Results of AD studies selecting for 48 or screening out 49 a history of MD suggest disturbance of the stress system in some AD patients may be due to the well-documented disturbance of this system in MD.

The fraction of variance accounted for by the subset of SNPs in the GABA receptor genes selected was 3.8 %, while the fraction of variance accounted for by the subset of SNPs in the six selected glutamate receptor genes was 3.3 %. We note that these estimates are comparable to previously published estimates for the contributions of multiple variants in distinct neurotransmitter systems to a different psychiatric phenotype (ADHD) 50. For comparison purposes, the fraction of variance that would be explained by comparably sized subsets of SNPs, randomly selected from those of comparable minor allele frequencies, would be 1.1% and 0.6% respectively. We think that a systematic search of other gene sets (including those related to metabolism of alcohol and taste) along the lines performed here may yield further useful results.

One concern about this study may be that the genes were a priori selected to be likely candidates for association with addictive behavior. We agree that this means it is more likely that we have higher levels of association than in a GWAS. However this caveat applies equally well to the individual tests, which as noted above, do not result in experiment-wide significance. And in any event a priori selection does not invalidate the statistical evidence of an independent study.

This method is likely to prove most useful in complex biomedical phenotypes where information exists about plausible pathophysiological mechanisms. Particularly when applied to genome wide data, the task will be to develop informed sets of candidate pathways that can be tested as plausible a priori hypotheses. Defining the boundaries of individual pathways may prove challenging as biological pathways are rarely discrete. Nonetheless, if the CDCV hypothesis is correct, genetic variation for many common diseases may be scattered throughout a number of genes within particular biological pathways. Testing a limited number of such pre-defined pathways is likely to lead to a substantial gain in statistical power if it is possible to nominate the likely candidate pathways. Such a gain in power could significantly reduce the need for very large samples for GWAS needed to compensate for the multiple testing burden.

In summary, we see the most important results here as methodological. We have developed and tested a network based approach to clarifying the genetic substrate of complex diseases. Applying this to a pre-selected set of candidate genes genotypes in cases with AD and controls, we found evidence that variants in four neurotransmitter pathways (norepinephrine, glutamate, GABA, and CRH) significantly contributed to risk for AD. These results are broadly consistent with the prior literature and suggest that network analyses may be a useful addition to the more standard approaches to clarifying the genetic basis of complex disorders, especially when prior biological hypothesis of disease etiology can lead to predicted gene pathways.