|Home | About | Journals | Submit | Contact Us | Français|
Recent genome-wide association studies (GWAS) have implicated a range of genes from discrete biological pathways in the aetiology of autism. However, despite the strong influence of genetic factors, association studies have yet to identify statistically robust, replicated major effect genes or SNPs. We apply the principle of the SNP ratio test methodology described by O'Dushlaine et al to over 2100 families from the Autism Genome Project (AGP). Using a two-stage design we examine association enrichment in 5955 unique gene-ontology classifications across four groupings based on two phenotypic and two ancestral classifications. Based on estimates from simulation we identify excess of association enrichment across all analyses. We observe enrichment in association for sets of genes involved in diverse biological processes, including pyruvate metabolism, transcription factor activation, cell-signalling and cell-cycle regulation. Both genes and processes that show enrichment have previously been examined in autistic disorders and offer biologically plausibility to these findings.
Autism is a complex neurodevelopmental disorder characterized by impairments of varying severity in the three core areas of communication, social interaction and repetitive behaviour. Population prevalence of autism is approximately 15–20 per 10000 with all autism spectrum disorders (ASD) estimated at 60 in 10000 children.1, 2 The role of genetic factors in the development of autism is undisputed. Heritability has been estimated as high as 91–93% using a multi-threshold liability model.3 However, despite the strong influence of genetic factors, autism linkage studies and association studies of common SNPs have not identified any genes of major effect. Recent genome-wide association studies (GWAS) have implicated a number of genes from discrete biological pathways in the aetiology of autism.4, 5, 6 In a recent study by the AGP using these data, we identified genome-wide significant association with MACROD2.7 However, we did not observe strong marker-wise associations within the cadherin gene region (CDH9, CDH10) or the TAS2R1, SEMA5A region that were highlighted in the work of Wang and co-workers,4 Ma and co-workers5 and Weiss and co-workers.6 In addition to identifying genome-wide significant association, it can be hypothesized that additional true vulnerability loci may exist within the nominal to modest range of statistical significance and confer risk to the disorder.8 A milieu of nominal to modestly associated risk variation fits with a polygenic model of disease and presents additional challenges for the identification of patterns of association within expected experimental noise.9
One promising approach is to examine association enrichment within ‘pathways' or groups of genes. The underlying hypothesis of association enrichment analysis is that functional polymorphisms that exist within a group of biologically inter-related genes are in essence ‘disrupting' the normal functioning of the biological process of the pathway. Consequently, one can consider the biological process, rather than the individual gene or SNP, in the development of the disease/disorder. By examining the ratio of association signals within a group of genes, we can determine whether there is enrichment of the signal above that expected by chance. This strategy decreases the multiple-testing burden that accompanies GWAS, and can have increased power to observe association.
A number of pathway-based methodologies have been developed to examine gene enrichment in association data (reviewed by Wang et al10). These include gene ranking algorithms,11 gene-enrichment algorithms, for example, ALIGATOR (Association LIst Go AnnoTatOR)9 and SNP enrichment approaches such as the SNP ratio test (SRT).12 The SRT provides a formal test of whether markers within pre-defined pathways show enrichment in association signal over that expected by chance alone. For case–control data, the basic algorithm underpinning the SRT is to first calculate the ratio of the number of nominally associated SNP markers within a pathway to the total number of markers within the pathway. Significance is assigned through a case-randomisation permutation routine, which takes account of the linkage disequilibrium between markers.
To apply the SRT to family-based data, we are unable to perform standard case randomisation; therefore, a pseudo-sibling model is generated from the alleles that are not transmitted to the proband. A proband-randomisation procedure is performed within the family, whereby the affection status of the offspring (case and pseudo-sibling) is permuted. This method allows retention of the linkage disequilibrium structure within the families and retains the advantages of the transmission disequilibrium test (TDT) design for the family-based association. In this study, we chose the SRT over other approaches for a number of reasons. Firstly, as the SRT retains all of the markers from the association analysis, it is sensitive to more than one true association signal per gene, and therefore gains information in the presence of allelic heterogeneity. Secondly, the SRT's use of multiple association signals across a gene as opposed to a single maximum signal limits potential genotyping artefact effects. Genotyping error at a single point may highlight a gene erroneously in a maximum signal design where this becomes the only observation. However, taking the ratio of all signals across a gene restricts the impact of single points of error as they are more likely to be diluted across the gene. Thirdly, the SRT also controls for gene size and linkage disequilibrium effects by permuting case-ness independently of genotype, consequently maintaining the same recombination patterns. Approaches that do not apply a gene-wise correction to GWAS data can show inflated signals for pathways containing larger genes. This is often the case in brain expressed pathways that are enriched for larger genes such as cell-surface receptors and can lead to misinterpretation of any association enrichment. Finally, as the SRT uses an SNP-wise association statistic over a gene-wise association statistic, we have sufficient observations to examine pathways that may contain fewer genes. Thereby, we are able to examine discrete ‘niche' pathways as well as larger, more diverse gene sets for enrichment in the GWAS.
For this study, we use gene-set lists derived from the gene-ontology (GO) (http://www.geneontology.org) database to examine whether association enrichment is present in a cohort of individuals from the Autism Genome Project (AGP) with a diagnosis of autistic disorder.
The individuals examined in this study were collected as part of the AGP Consortium genome analysis project. The AGP represents more than 50 centres in North America and Europe. Subjects with known karyotypic abnormalities, fragile X mutations or other known genetic disorders were excluded. Diagnostic and ancestral definitions were as previously reported by this group.7 Briefly, families are grouped into two nested diagnostic classes (Strict and Spectrum) based on proband diagnostic measures. To qualify for the Strict class, affected individuals met the criteria for autism on both primary diagnostic instruments: the Autism Diagnostic Interview-Revised (ADI-R13) and Autism Diagnostic Observation Schedule (ADOS14). ADI-R-based diagnostic classification of subjects as ASD followed the criteria published by Risi and co-workers.15 Specifically, individuals who almost met the ADI criteria for autism were classified as ASD if: (1) they met the criteria on social and either communication or repetitive behaviour domains; or (2) met the criteria on the social domain and were within two points of criteria for communication, or met the criteria on the communication domain and were within two points of social criteria, or within one point on both social and communication domains. The Spectrum class included all individuals who met the Strict criteria and those individuals who were classified as ASD or autism on both the ADI-R and ADOS or who were not evaluated on one of the instruments, but were diagnosed with autism on the other instrument. A summary of the sample sizes for the discovery and replication datasets for each DiagnosticAncestry subset is shown in Table 1.
As described elsewhere,7 ancestry for these individuals was determined for the proband by using 5239 widely spaced, independent SNPs that had a genotype completion rate of ≥99.9%. The software used was Spectral-GEM,16 which estimated five significant dimensions of ancestry. Subsequent clustering on dimensions of ancestry identified nine clusters; five clusters were used to describe European ancestry and the remaining clusters best reflect Asian, African (East/West) and Latin American origins. The All ancestry class included all individuals including those who met the European ancestry criteria.
The discovery sample were genotyped using the Illumina Infinium 1M-single SNP microarray, the replication sample were genotyped on a either the Illumina Infinium 1M-single SNP microarray or the Illumina 1M-duo microarray. All quality control (QC) procedures were maintained across datasets; in addition, QC marker sets from both the discovery and replication datasets were matched and only those markers meeting QC for both the discovery and replication datasets were carried forward to analysis. Additional QC details are described elsewhere.7 A total of 856932 SNPs passed QC on both the discovery and replication sample. TDT statistics were generated using PLINK v.1.07.17
The pedigree SRT (pedSRT) is a modification to the SRT described by O'Dushlaine and co-workers,12 which is applicable to family-based data. Briefly, the SRT tests the ratio of the number of associated SNPs to the total number of SNPs in a pre-defined set of genes. A marker is considered ‘associated' if the association statistic is observed below a given threshold. The threshold used is arbitrary, but is set by default at an unadjusted P≤0.05. The significance of the ratio is determined through permutation using an empirical P-value derived from the proportion of the ratios for the permuted datasets that are greater than or equal to the observed ratio.12 We performed 10000 permuted GWAS analyses for each of the diagnostic, ancestry strata for both the discovery and replication datasets. The pedSRT determines association using the TDT,18 as implemented in PLINK.17 In a case–control model, permutation is performed using case randomisation. In the TDT design, case randomisation is performed by creating a pseudo-sibling. The pseudo-sibling is created from the non-transmitted alleles from the parents. Within each permutation cycle either the proband or pseudo-sibling is considered the ‘case'. Alternate case randomisation for the TDT is implemented in PLINK, using the alternate phenotype routine.
It is important to note that to reduce type-I error in the SRT due to inflation of the original association signal, for each permutation ‘associated' SNPs are assigned according to their rank in the dataset.12 In short, the numbers of SNPs (T) that meet the ‘associated' threshold are calculated from the primary dataset. For each permuted dataset, the top T SNPs are termed ‘associated'.
All SNP ratio statistics were calculated using custom scripts in STATA version 10 (Stata Corp., College Station, TX, USA).
Individual SNP codes from the Illumina 1M Infinium SNP array platform were updated to reflect build 130 of dbSNP. SNPs were assigned to genes using gene criteria from the dbSNP/NCBI criteria: namely, if the SNP resides within the locus containing the gene transcript including 2kb 5′ and 500bp 3′ of the transcript. The gene assignment protocol was performed using the NCBI criteria and facilitated using the file b130_SNPContigLocusId_36_3.bcp available at ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/database/org
Gene sets were described using the GO database (http://www.geneontology.org).19 Gene lists were obtained from the OBO format 1.2 database release available from http://www.geneontology.org/GO.downloads.ontology.shtml (build release date 15 December 2009). GO terms are structured in a semihierarchical relationship within the cellular component, molecular function and biological process nodes. Daughter ontology terms are more specialized and parent ontology terms are less specialized. However, unlike a hierarchy, a term may have more than one parent term.
Parent terms were populated by their daughter terms to describe a composite list of genes for each term. SNP ratios were calculated on GO terms with greater than 20 SNPs, but less than 2000 SNPs, and greater than 1 gene, but no more than 1000 genes. A total of 6853 GO terms met these criteria. To account for identity of terms, we merged those GO terms containing identical gene lists; in total the list of unique terms is 5955.
As mentioned above, the GO terms used in this study can show considerable overlap owing to term redundancy, biological overlap and the hierarchical nature of the database. Simulations were performed to calculate the null distribution and subsequent expectancy for the total number of associated GO terms at a given threshold in a single study given the GO terms used.
We performed 1000 pedSRT permutations on a case-randomized sample derived from 1248 families from the discovery dataset. A GWAS TDT was performed on each dataset followed by pedSRT, using 10000 additional permutations on the 5995 GO terms. For each of the 1000 original permutations, the proportion of the 5955 GO terms that met a significance threshold of P≤0.05 in the subsequent 10000 was calculated. The mean proportion across the 1000 permutations was used to predict the expected number of associated GO terms in a dataset.
Visual representation of overlap in enriched GO terms was performed using the EnrichmentMap (http://baderlab.org/Software/EnrichmentMap20) plugin for Cytoscape 2.8.0 (http://www.cytoscape.org/21). Consistent with the author's recommendations for use with the GO database, nodes were joined if the overlap coefficient was ≥0.5.
Across all analysis in the discovery dataset, 1035 unique GO terms show association enrichment at SRT-P-value≤0.05. Examination of those GO terms that show strong enrichment (SRT-P-value<0.001) highlights diverse processes such as the regulation of cell division (mitosis and meiosis), ribosome processing and apoptosis. A visual representation of enriched pathways is shown in Supplementary Figure 1. A summary of the total number of GO terms that show enrichment at SRT-P-value≤0.05 is given in Table 2. On the basis of simulated data, 4.46% (SD=0.8%) of the 5995 unique but non-independent pathways are expected to be associated at SRT-P≤0.05 level. Given this level, we would expect 267 GO terms to be associated per experiment. To provide a greater distinction of potentially important GO terms, we examined the overlap of enriched GO terms in an independent replication dataset. On the basis of 4.46% of GO terms showing enrichment, we would expect to observe replication for 12 of the 5995 pathways. All individual discovery samples show more GO terms associated than would be expected by chance (see Expected 1; Table 2). Moreover, the overlap between the discovery and replication sample also show enrichment over what would be expected by chance (see Expected replication 2; Table 2). When we use a more cautious interpretation based on the total number of observed associated GO terms in the discovery data and the predicted replication of 4.46% we would expect to replicate is between 15 and 17 pathways (see Expected 3; Table 2). Under this model, we still show enriched replication for each ‘DiagnosisAncestry' groupings. Overall, compared to simulated data we observe between 1.5- and 3.2-fold enrichment in the overlap of pathways in the discovery and replication dataset above what can be expected by chance.chance.
A summary of the replicated pathways, summary statistics, gene number and genes tagged in this analysis is shown in Tables 4, ,5,5, ,6,6, ,77 (full lists of replicated pathways can be found in Supplementary Tables 1A–C). A total of 88 unique GO terms were shown to be replicated within analytic groupings (see Supplementary Table 2), 22 GO terms were replicated within two of the analytic groupings and four GO terms were replicated within three of the analytic groupings (see Table 3). Replication was only considered within strata, such that, for example, GO terms identified in the discovery StrictEuropean analyses were examined in the StrictEuropean replication dataset. The four GO terms that show enrichment across three groupings are as follows: GO:0006090, GO:0032872, GO:0032874 and GO:0042156, involved in pyruvate metabolism, regulation of the mitogen-activated protein kinase (MAPK) cascade and zinc-mediated transcriptional activation. A visual representation of replicated enriched pathways is shown in Supplementary Figure 2.
The interpretation of GWAS data purely on the strength of association data is challenging where the distribution of association is close to or barely exceeding what is expected by the number of tests. In the absence of clear association enrichment across the entire dataset, interpretation has relied on rank order or via the application of suboptimal significance thresholds that juggle type-I and type-II error. The principle of association enrichment approaches is to discover whether within this milieu of data there are underlying patterns to the association. In these approaches, we ask whether SNPs that are linked to genes of common function show greater proportion of nominal association than expected by chance. Although a modest association signal at an individual SNP within a gene may not warrant further investigation, the cumulative association of SNPs within a gene family may offer insight into the biology of the disorder.
Gene enrichment approaches have been primarily developed to aid interpretation of data from microarray expression studies. In this context, each gene is tagged by either one or a small number of probes regardless of gene size. However, when applying these technologies to SNP-based data, we do not measure gene-wise variation or gene-wise association; instead, we can potentially examine multiple points of association at any given gene using many tagging SNPs. This brings additional challenges and bias. When applying association enrichment, we must account for and correct for these potential bias in these data. Firstly, when examining larger genes, we utilize more SNP markers to tag the variation than for smaller genes. If we choose a maximum association signal approach per gene, we observe by chance an inflated signal for the larger genes. By calculating the ratio of associated to not associated SNPs, we can adjust each GO term to the total number of SNPs examined per GO term. Secondly, where multiple markers tag a gene, one might observe multiple strong association signals due to strong linkage disequilibrium between the associated markers. To reduce this effect, we calculate significance of the data through permutation. Permutation is performed by case randomisation within families where a pseudo-control sibling is created from the alleles that are not transmitted to the proband. By using the non-transmitted alleles, we retain the linkage disequilibrium structure across the genome, thereby retaining linkage-disequilibrium-related inflation in the original association signal.
We have applied the SRT to family-based data from the AGP to identify 88 gene sets from the GO database that show a replicated enrichment for association signal. Of the overlapping GO terms, we observe enrichment in sets involved in diverse biological processes, including pyruvate metabolism, transcription factor activation, cell-signalling and cell-cycle regulation.
One of the strongest findings from the discovery and replication findings was observed across the ‘Strict diagnosisAll ancestries' grouping for the GO term GO:0031146; SCF-dependent proteasomal ubiquitin-dependent protein catabolic process (Discovery SRT-P=0.0001; Replication SRT-P=0.0009). GO:0031146 is described by only two genes (FBXO31 and FBXO6). Both genes are members of the F-box protein family, which are involved in a variety of molecular and cellular functions, including protein degradation, synapse formation and circadian rhythm.22 FBXO6 has also been suggested as a putative biomarker for autism,23 as one of 13 genes highlighted in the work of Nishimura and Brown24 who show differential expression at this gene in the lymphoblastoid cell lines from individuals with both the FMR1 mutation and autism compared with typically developing controls.
Those GO terms that show replication across multiple diagnostic and ancestral groups are also worth noting as they are robust to differences in sampling used in our analyses. Four replicated GO terms were observed in three analytic groupings (see Table 3). These include GO:0006090, GO:0032872, GO:0032874 and GO:0042156. GO:0006090 (pyruvate metabolic process) describes a group of 39-tagged genes (see Supplementary Table 3) covered by 589 SNPs. These genes are involved in the biological processes connecting the chemical reactions and pathways involving pyruvate. Pyruvate metabolism is a component of the energy metabolism pathway, which has received considerable attention with respect to autism. The biological plausibility of the pyruvate metabolic process association enrichment is supported by numerous studies showing evidence of aberration in pyruvate levels in individuals with autism.25 The GO term GO:0042156 (zinc-mediated transcriptional activator activity) describes a group of three genes tagged by 37 SNPs (MTF1, RNF4 and ZNF384). One of the constituent genes, MTF1, human metal-regulatory transcription-factor-1, has previously warranted investigation as putative candidate gene for autistic disorder under an environmental exposure model of autism.26 Finally, GO:0032872 (regulation of stress-activated MAPK cascade) and GO:0032874 (positive regulation of stress-activated MAPK cascade), which differ by a single gene (see Supplementary Table 3), describe 10 and 9 genes, and 122 and 116 SNPs, respectively. These pathways are involved in increasing the signalling of the stress-related MAPK signalling pathway. Stress-activated MAPKs are thought to have a critical role in modulating inflammation, DNA damage response, apoptosis in cancer27 and negative regulation of cell-cycle progression.28, 29 Cell-cycle progression and DNA damage response are also highlighted in enriched replicated GO terms in these analyses, for example, GO:0032404 (mismatch repair complex binding) and GO:0031571 (G1/S DNA damage checkpoint).
In a recent study by this group, we explored enrichment in GO terms for rare deleted CNVs.30 Using individuals from the Discovery Group, we identified 24 enriched GO terms that show enrichment in rare CNV at FDR q<0.05 that highlighted five biological domains: namely cell proliferation, cell projection and motility, MHC-I, GTPase/RAS signalling and kinase activation/regulation. We do not observe any overlap between the 88 gene sets showing replicated enrichment in the GWAS data with the 24 significant GO terms identified for rare structural variation. However, we do observe some overlap for GO terms enriched only in the discovery dataset. These include overlap in ‘cell migration', ‘cell motility', ‘cell morphogenesis' and GO terms identified as having a role in protein kinase regulation.
We can take some encouragement that highlighted pathways are supported in the autism literature. We have emphasized biological plausibility of some of these pathways with autism and ASD. However, one major caveat when interpreting these data is whether this overlapping evidence reflects the considerable literature surrounding autism research and is therefore coincident, or is biologically meaningful concordance.
Pathway approaches, such as the SRT and pedSRT, can also be applied to research questions using candidate gene list. Candidate genes rely on the selection of genes and markers based on previous knowledge of biology, function and position of the gene or marker. The pathway approach in the form used in this manuscript applies a ‘hypothesis-free' design, in which we examine all GO terms regardless of putative role. In a recent autism GWAS described by Wang et al,4 the authors applied a hypothesis-testing candidate gene approach using their own methodology11 to examine whether a group of cadherin and neurexin genes showed enrichment in their association data. The authors conclude that there was association enrichment for both a group of cadherin, and cadherin plus neurexin genes (P=0.02, 0.004, respectively). We applied our approach to these gene lists in our data (data not shown). Using the pedSRT, which differs in statistical method and gene-to-SNP assignment to that of Wang and co-workers, we do not observe significant association signal enrichment in either the discovery or replication dataset for any of the analytic groupings.
To further explore potential overlap of our data and other GWAS, we examined whether previously implicated genes from recent autism GWAS were present in the GO terms identified in this study. None of the genes that overlap with the top-associated SNPs from previous GWAS described by Wang and co-workers4 (CDH22, CTNNA3, DMD, FEZF2, LOC100132914, LRRC1 and SYT17) and Weiss and co-workers6 (ACTN2, ADA, CENPC1, CRIM1, CTNNA3, CUGBP2, GAS2, IQGAP2, JARID2, SGCD and XG) appeared in the 88 unique GO terms showing overlap in these analyses. Moreover, we do not observe overlap with those genes highlighted by the authors as residing close to their maximal association peaks, namely SEMA5A, TAS2R1 and CDH9, CDH10.
The GO database is continuously updated as evidence is gathered on gene biology. The build of the database used in these analyses contains information on 17703 genes, compared with less than 5000 for databases such as KEGG. However, not all genes are tagged to GO terms. This is exemplified by the MACROD2 gene, which contained SNPs showing the strongest association signal from our previous GWAS analyses.7 Over time more information will be gathered on the biological role and interactions between these genes to further annotate these terms.
In addition to single gene effects such as MACROD2, data presented in this analysis may offer some additional insight into biological processes, within which genetic risk for autism may lie. This can include hypothesis-free gene lists such as those in the GO dataset, or more hypothesis-driven candidate gene lists highlighting previous linkage, association or biology. The application of pedSRT to our GWAS data has highlighted biological processes previously implicated in autism and offers impetus to re-examine these processes based on evidence from genome-wide investigation. Association enrichment analysis provides additional evidence from GWAS data to identify genetic risk variants and genes and prioritize biological processes for further research into areas such as biomarker discovery, gene–gene interaction analyses and identification of putative drug targets.
We gratefully acknowledge the families participating in the study and the main funders of the AGP: Autism Speaks (USA), the Health Research Board (HRB, Ireland; AUT/2006/1, AUT/2006/2, PD/2006/48), The Medical Research Council (MRC, UK), Genome Canada/Ontario Genomics Institute and the Hilibrand Foundation (USA). Additional support for individual groups was provided by the US National Institutes of Health (NIH Grants: HD055751, HD055782, HD055784, MH52708, MH55284, MH061009, MH06359, MH066673, MH080647, MH081754, MH66766, NS026630, NS042165, NS049261), the Canadian Institutes for Health Research (CIHR), Assistance Publique – Hôpitaux de Paris (France), Autism Speaks UK, Canada Foundation for Innovation/Ontario Innovation Trust, Deutsche Forschungsgemeinschaft (Grant: Po 255/17-4) (Germany), EC Sixth FP AUTISM MOLGEN, Fundação Calouste Gulbenkian (Portugal), Fondation de France, Fondation FondaMental (France), Fondation Orange (France), Fondation pour la Recherche Médicale (France), Fundação para a Ciência e Tecnologia (Portugal), the Hospital for Sick Children Foundation and University of Toronto (Canada), INSERM (France), Institut Pasteur (France), the Italian Ministry of Health (convention 181 of 19 October 2001), the John P Hussman Foundation (USA), McLaughlin Centre (Canada), Ontario Ministry of Research and Innovation (Canada), the Seaver Foundation (USA), the Swedish Science Council, The Centre for Applied Genomics (Canada), the Utah Autism Foundation (USA) and the Wellcome Trust core award 075491/Z/04 (UK). DP is supported by fellowships from the Royal Netherlands Academy of Arts and Sciences (TMF/DA/5801) and the Netherlands Organization for Scientific Research (Rubicon 825.06.031). SWS holds the GlaxoSmithKline-CIHR Pathfinder Chair in Genetics and Genomics at the University of Toronto and the Hospital for Sick Children (Canada).
EMK and COD developed the principle of the SRT experiments. RJLA designed the study, developed and implemented the pedSRT methodology and wrote the manuscript. EMK, COD, BLY, MG and LG aided in manuscript preparation. RJLA, EMK, COD, BY, EP, JDB and JSS discussed research strategies and data through the 'pathway-based analysis working group‘. Additional intellectual support and guidance was provided through the AGP including BD, ADP, EHC, PS, JTG, CK, KW, HH and EM.
The authors declare no conflict of interest.
Supplementary Information accompanies the paper on European Journal of Human Genetics website (http://www.nature.com/ejhg)