|Home | About | Journals | Submit | Contact Us | Français|
DAvH and CW designed, co-ordinated and led the study. Experiments were performed in the labs of CW, DAvH, CAM, PD and PMG. Major contributions were: 1) DNA sample preparation - PCAD, GT, KAH, JR, AZ, PS 2) genotyping - PCAD, GT, KAH, AC, JR, RG; 3) expression data generation - HJMG, LHvdB, RAO, RKW, LF; 4) case/control association analyses - PCAD, GT, LF, JCB, DAvH; 5) expression analyses - LF, GAH, RSNF; 6) manuscript preparation - PCAD, GT, LF, RSNF, GAH, JCB, CW, DAvH. Other authors contributed variously to sample collection and all other aspects of the study. All authors reviewed the final manuscript.
We performed a second-generation genome wide association study of 4,533 celiac disease cases and 10,750 controls. We genotyped 113 selected SNPs with PGWAS<10−4, and 18 SNPs from 14 known loci, in a further 4,918 cases and 5,684 controls. Variants from 13 new regions reached genome wide significance (Pcombined<5×10−8), most contain immune function genes (BACH2, CCR4, CD80, CIITA/SOCS1/CLEC16A, ICOSLG, ZMIZ1) with ETS1, RUNX3, THEMIS and TNFRSF14 playing key roles in thymic T cell selection. A further 13 regions had suggestive association evidence. In an expression quantitative trait meta-analysis of 1,469 whole blood samples, 20 of 38 (52.6%) tested loci had celiac risk variants correlated (P<0.0028, FDR 5%) with cis gene expression.
Celiac disease is a common heritable chronic inflammatory condition of the small intestine induced by dietary wheat, rye and barley, as well as other unidentified environmental factors, in susceptible individuals. Specific HLA-DQA1 and HLA-DQB1 risk alleles are necessary, but not sufficient, for disease development1,2. The well defined role of HLA-DQ heterodimers, encoded by these alleles, is to present cereal peptides to CD4+ T cells, activating an inflammatory immune response in the intestine. A single genome wide association study (GWAS) has been performed in celiac disease, which identified the IL2/IL21 risk locus1. Subsequent studies probing the GWAS information in greater depth have identified a further 12 risk regions. Most of these regions contain a candidate gene functional in the immune system, although only in the case of HLA-DQA1 and HLA-DQB1 have the causal variants been established3-5. Many of the known celiac loci overlap with other immune-related diseases6. In order to identify additional risk variants, particularly of smaller effect size, we performed a second-generation GWAS using over six times as many samples as the previous GWAS and a denser genome-wide SNP set. We followed up promising findings in a large collection of independent samples.
The GWAS included five European celiac disease case and control sample collections including the previously reported celiac disease dataset1. We performed stringent data quality control (Online Methods), including calling genotypes using a custom algorithm on both large sample sets, and where possible cases and controls together (Online Methods). We tested 292,387 non-HLA SNPs from the Illumina Hap300 marker set for association in 4,533 celiac disease cases and 10,750 controls of European descent (Table 1). A further 231,362 additional non-HLA markers from the Illumina Hap550 marker set were tested for association in a subset of 3,796 celiac disease cases and 8,154 controls. All markers were from autosomes or the X chromosome. Genotype call rates were >99.9% in both datasets. The overdispersion factor of association test statistics, λGC=1.12, was similar to that observed in other GWAS of this sample size7,8. Findings were not substantially altered by imputation of missing genotypes for 737 celiac disease cases genotyped on the Hap300 BeadChip and corresponding controls (Table 1, collection 1). Here we present results for directly genotyped samples, as around half the additional Hap550 markers cannot be accurately imputed from Hap300 data9 (including the new ETS1 locus finding in this study). Results for the top 1000 markers are available in Supplementary Data 1, but because of concerns regarding identity detection of individuals10, results for all markers are available only on request to the corresponding author.
For follow-up, we first inspected genotype clouds for the 417 non-HLA SNPs meeting PGWAS<10−4, being aware that top GWAS association signals may be enriched for genotyping artefact, and excluded 22 SNPs from further analysis using a low threshold for possible bias. We selected SNPs from 113 loci for replication. Markers that passed design and genotyping quality control included: a) 18 SNPs from all 14 previously identified celiac disease risk loci (including a tag SNP for the major celiac disease associated HLA-DQ2.5cis haplotype1); b) 13 SNPs from all 7 novel regions with PGWAS<5×10−7; c) 86 SNPs from 59 of 68 novel regions with 5×10−7> PGWAS >5×10−5 in stage 1; d) 14 SNPs from 14 of 30 novel regions with 5×10−5> PGWAS >10−4 in stage 1 (for this last category, we mostly chose regions with immune system genes). Two SNPs were selected per region for: regions with stronger association; regions with possible multiple independent associations; and/or containing genes of obvious biological interest. We successfully genotyped 131 SNPs in 7 independent follow-up cohorts comprising 4,918 celiac disease cases and 5,684 controls of European descent. Genotype call rates were >99.9% in each collection. Primary association analyses of the combined GWAS and follow-up data were performed with a two-sided 2×2×12 Cochran-Mantel-Haenszel test.
The HLA locus and all 13 other previously reported celiac disease risk loci showed evidence for association at a genome wide significance threshold (Pcombined<5×10−8) (Table 2). We note that some loci were previously reported using less stringent criteria (e.g. the P<5×10−7 recommended by the 2007 WTCCC study11), but that in the current, much larger sample set, all known loci meet recently proposed P<5×10−8 thresholds12,13.
We identified 13 novel risk regions with genome-wide significant evidence (Pcombined <5×10−8) of association, including regions containing the BACH2, CCR4, CD80, CIITA/SOCS1/CLEC16A, ETS1, ICOSLG, RUNX3, THEMIS, TNFRSF14, and ZMIZ1 genes which are of obvious immunological function (Table 2). A further 13 regions met ‘suggestive’ criteria for association (either 10−6> Pcombined >5×10−8 and/or PGWAS<10−4 and Pfollowup<0.01). These regions also contain multiple genes of obvious immunological function, including CD247, FASLG/TNFSF18/TNFSF4, IRF4, TLR7/TLR8, TNFRSF9 and YDJC. Six of the 39 non-HLA regions show evidence for the presence of multiple independently associated variants in a conditional logistic regression analysis (Supplementary Table 2).
We tested the 40 SNPs with the strongest association (Table 2) from each of the known genome-wide significant, new genome-wide significant, and new suggestive loci for evidence of heterogeneity across the 12 collections studied. Only the HLA region was significant (Breslow-Day test P<0.05 / 40 tests, rs2187668 P=4.8×10−8) which is consistent with the well described North-South gradient in HLA allele frequency in European populations, and more specifically for HLA-DQ in celiac disease14.
We observed no evidence for interaction between each of the 26 genome-wide significant non-HLA loci, which is consistent with what has been reported for complex diseases so far. However, we did observe weak evidence for lower effect sizes at non-HLA loci in high risk HLA-DQ2.5cis homozygotes, similar to what has been previously observed in type 1 diabetes7.
To obtain more insight into the functional relatedness of the celiac loci, we applied GRAIL, a statistical tool that utilizes text mining of PubMed abstracts to annotate candidate genes from loci associated with common disease risk15,16. In order to assess the sensitivity of this tool (using known loci as a positive control), we first performed a ‘leave-one-out’ analysis of the 27 genome-wide significant celiac disease loci (including HLA-DQ). GRAIL scores of Ptext<0.01 were obtained for 12 of the 27 loci (44% sensitivity, Table 2). Factors that limit the sensitivity of GRAIL include biological pathways being both known (a 2006 dataset is used to avoid GWAS era studies), and published in the literature. We then applied GRAIL analysis, using the 27 known regions as a seed, to all 49 regions (49 SNPs) with 10−3 >Pcombined >5×10−8 and obtained GRAIL Ptext <0.01 for 9 regions (18.4%). As a control, only 5.5% (279 of 5033) of randomly selected Hap550 SNPs reached this threshold. In addition to the five ‘suggestive’ loci shown in Table 2, GRAIL annotated four further interesting gene regions of lower significance in the combined association results: rs944141/PDCD1LG2 (Pcombined=4.4×10−6), rs976881/TNFRSF8(Pcombined =2.1×10−4), rs4682103/CD200/BTLA(Pcombined=6.8×10−6) and rs4919611/NFKB2 (Pcombined=6.1×10−5). There appeared to be further enrichment for genes of immunological interest which are not GRAIL annotated in the 10−3>Pcombined>5×10−8 significance window, including rs3828599/TNIP1 (Pcombined=1.55×10−4), rs8027604/PTPN9 (Pcombined=1.4×10−6), rs944141/CD274 (Pcombined=4.4×10−6). Some of these findings, for which neither genome-wide significant nor suggestive association is achieved, are likely to comprise part of a longer tail of disease predisposing common variants, of weaker effect sizes. Definitive assessment of these biologically plausible regions would require genotyping and association studies using much larger sample collections than the present study.
We previously showed considerable overlap between celiac disease and type 1 diabetes risk loci17, as well as celiac disease and rheumatoid arthritis18, and more generally, there is now substantial evidence for shared risk loci between the common chronic immune mediated diseases6. To update these observations, we searched ‘A Catalog of Published Genome Wide Association Studies’ (18 Nov 2009)19 and the HuGE database20. We found some evidence (requiring a published association report of P<1×10−5) of shared loci with at least one other inflammatory or immune mediated disease for 18 of the current 27 genome-wide significant celiac risk regions. We defined shared regions as the broad LD block, however different SNPs are often reported in different diseases, and at only three of the 18 shared regions are associations across all diseases with the same SNP or a proxy SNP in r2>0.8 in HapMap CEU. Currently 9 regions appear celiac disease specific and may reflect distinctive disease biology including the regions containing rs296547 and rs9792269, and the regions around CCR4, CD80, ITGA4, LPP, PLEK, RUNX3 and THEMIS. In fact, locus sharing between diseases is probably greater, due to both stochastic variation in results from sample size limitations, and regions with a genuinely stronger effect size in one disease and weaker effect size in another.
Genetic variation in ETS1 has recently been reported to be associated with systemic lupus erythematosus (SLE) in the Chinese population, although is not associated with SLE in European populations21. The most strongly associated celiac (European population) SNP rs11221332 and the most strongly associated SLE (Chinese population) SNP rs6590330 map 70kb apart. Inspection of the HapMap phase II data shows broadly similar linkage disequilibrium patterns between Chinese (CHB) and European (CEU) populations in this region, with the two associated SNPs in separate non-adjacent linkage disequilibrium blocks. Thus distinct common variants within the same gene are predisposing to different autoimmune diseases across different ethnic groups.
Celiac risk variants in the HLA alter protein structure and function4. However we identified only four non-synonymous SNPs with evidence for celiac disease association (PGWAS<10−4) from the other 26 genome-wide significant associated regions (rs3748816/MMEL1, rs3816281/PLEK, rs196432/RUNX3, rs3184504/SH2B3). Although comprehensive regional resequencing is required to test the possibility that coding variants contribute to the observed association signals, more subtle effects of genetic variation on gene expression are the more likely major functional mechanism for complex disease genes. With this in mind, we performed a meta-analysis of new and published genome-wide expression quantitative trait loci (eQTL) datasets comprising 1,469 human whole blood (PAXgene) samples reflecting primary leucocyte gene expression. We applied a new method, Transcriptional Components, to remove a substantial proportion of inter-individual non-genetic expression variation and performed eQTL meta-analysis on the residual expression variation (Online Methods).
We assessed 38 of the 39 genome-wide significant and suggestive celiac disease associated non-HLA loci (Table 2) for cis expression - genotype correlations. We tested the SNP with the strongest association from each region. However for five regions the most associated SNP was not genotyped in the eQTL samples (Hap300 data), instead for four of these we tested a proxy SNP (r2>0.5 in HapMap CEU). In addition, for six loci showing evidence of multiple independent associations in conditional regression analyses, we tested a second SNP showing independent celiac disease association for eQTL analysis. In total we assessed 44 independent non-HLA SNP associations in peripheral whole blood samples genotyped on the Illumina Hap300 BeadChip and either Illumina Ref8 or HT12 expression arrays, correlating each SNP with data from gene-probes mapping within a 1Mb window.
We identified significant (Spearman P<0.0028, corresponding to 5% false discovery rate) eQTLs at 20 of 38 (52.6%) non-HLA celiac loci tested (Table 3, Supplementary Figures 2 & 3). Some loci had evidence of eQTLs with multiple probes, genes or SNPs (Table 3). We assessed whether the number of SNPs with cis-eQTL effects out of the 44 SNPs that we tested, was significantly higher than expected. We observed that eQTL SNPs on average have a substantially higher MAF than non-eQTL SNPs in the 294,767 SNPs tested. In order to correct for this we selected 44 random SNPs that had an equal MAF distribution, and determined for how many of these MAF-matched SNPs eQTLs were observed. We observed a significantly higher number of eQTL SNPs (P=9.3 × 10−5, 106 permutations) amongst the celiac associated SNPs than expected by chance (22 observed eQTL SNPs, vs. 7.8 expected eQTL SNPs). Therefore the celiac disease associated regions are greatly enriched for eQTLs. These data suggest some risk variants may influence celiac disease susceptibility through a mechanism of altered gene expression. Candidate genes with a significant eQTL, where the peak eQTL signal and peak case/control association signal are similar (Supplementary Figure 3), include MMEL1, NSF, PARK7, PLEK, TAGAP, RRP1, UBE2L3 and ZMIZ1.
We also assessed co-expression of genes mapping within 500kb of SNPs showing strongest case/control association from the 40 genome-wide significant and suggestive celiac disease loci in an analysis of the 33,109 human Affymetrix Gene Expression Omnibus dataset. This analysis loses power to detect tissue specific correlations by use of numerous tissue types, but greatly gains power by the large sample size. We detected several distinct co-expression clusters (Pearson correlation coefficient between genes >0.5), including four clusters of immune-related genes which contain at least one gene from 37 of the 40 genome-wide significant and suggestive loci (Fig. 1). These data further demonstrate that genes from celiac disease risk loci map to multiple distinct immunological pathways involved in disease pathogenesis.
We previously reported that most celiac genetic risk variants mapped near genes that are functional in the immune system22, and this remains true for the 13 new genome-wide significant, and 13 new suggestive, risk variants from the current study. We can now refine these observations and highlight specific immunological pathways relevant to celiac disease pathogenesis:
The rs802734 LD block contains the recently identified gene THEMIS ‘THymus-Expressed Molecule Involved in Selection’. THEMIS plays a key regulatory role in both positive and negative T-cell selection during late thymocyte development23. Furthermore, the rs10903122 LD block contains RUNX3, a master regulator of CD8+ T lymphocyte development in the thymus24,25. TNFRSF14 (LIGHTR, rs3748816 LD block) has widespread peripheral leucocyte functions as well as a critical role in promoting thymocyte apoptosis26. The ETS1 transcription factor (rs11221332 LD block) is also active in peripheral leucocytes, however it is also a key player in thymic CD8+ lineage differentiation, acting in part by promoting RUNX3 expression27.
The importance of the thymus in autoimmune disease pathogenesis has been previously emphasised by the established role of thymectomy in the treatment of myasthenia gravis. In type 1 diabetes, it was shown that disease associated genetic variation in the insulin gene INS causes altered thymic insulin expression and subsequent T cell tolerance for insulin as a self-protein28. However, the importance of thymic T cell regulation has not been previously recognised in the aetiology of celiac disease. It is conceivable that the associated variants may alter biological processes prior to thymic MHC-ligand interactions. Alternatively it is now clear that exogenous antigen presentation and selection occurs in the thymus via migratory dendritic cells - this has been demonstrated for skin, and has been hypothesised for food antigens29,30. These findings suggest research into immuno-/pharmaco-logical modifiers of T cell tolerance more generally in autoimmune diseases.
Although the association signal at rs5979785 (Pcombined=6.36×10−8) in the TLR7/TLR8 region is just outside our genome wide significance threshold, we observe a strong effect of rs5979785 on TLR8 expression in whole blood. Both TLRs recognise viral RNA. Taken together with the recent observation of rare loss of function mutations in the enteroviral response gene IFIH1 protective against type 1 diabetes31, these findings suggest viral infection (and the nature of the host response to infection) as a putative environmental trigger common to these autoimmune diseases.
This class of molecules controls the strength and nature of the response to T or B (immunoglobulin) cell receptor activation by antigens. We observe multiple regions with genes (CTLA4/ICOS/CD28, TNFRSF14, CD80, ICOSLG, TNFRSF9, TNFSF4) from this class of ligand-receptor pairs suggesting fine control of the adaptive immune response might be altered in at-risk individuals.
Our previous report discussed the function of the 2q11-12 interleukin receptor cluster (IL18RAP, etc), the 3p21 chemokine receptor cluster (CCR5, etc.) and the loci containing IL2/IL21 and IL12A22. We now report additional loci containing TNFSF18 and CCR4.
We estimate that the current celiac disease variants, including the major celiac disease associated HLA variant, HLA-DQ2.5cis, less common celiac disease associated haplotypes in the HLA (HLA-DQ8; HLA-DQ2.5trans; HLADQ2.2), and the additional 26 definitively implicated loci explain about 20% of total celiac disease variance, which would represent 40% of genetic variance, assuming a heritability of 0.5. A long tail of low effect size common variants, along with highly penetrant rare variants (both at the established loci and elsewhere in the genome), may substantially contribute to the remaining heritability.
We observed different haplotypes within the ETS1 region associated with coeliac disease in Europeans, and SLE in the Chinese population. We further note for some autoimmune diseases studied in European origin populations, that although the same LD block has been associated, the association is with a different haplotype. In some cases, the same variants are associated, but the direction of association is opposite (e.g. rs917997/IL18RAP in celiac disease versus type 1 diabetes). We believe further exploration of these signals may reveal critical differences in the nature of the immune system perturbation between these diseases.
Previous investigators have observed that only a small proportion of GWAS associations are coding variants, and have suggested that these may instead influence regulation of gene expression. Here, we show that over half the celiac disease associated variants are correlated with expression changes in nearby genes. This mechanism is likely to explain the function of some risk variants for other common, complex diseases. Further research, however, is needed to definitively determine at each locus both the celiac disease causal variants and their functional mechanisms.
We thank Celiac UK for assistance with direct recruitment of celiac disease individuals, and UK clinicians (L.C. Dinesen, G.K.T. Holmes, P.D. Howdle, J.R.F. Walters, D.S. Sanders, J. Swift, R. Crimmins, P. Kumar, D.P. Jewell, S.P.L. Travis, K. Moriarty) who recruited celiac disease blood samples described in our previous studies1,22. We thank the genotyping facility of the UMCG (J. Smolonska, P. van der Vlies) for generating part of the GWAS and replication data and the gene expression data; R. Booij and M. Weenstra are thanked for preparation of Italian samples. H. Ahola, A. Heimonen, L. Koskinen, E. Einarsdottir and K. Löytynoja are thanked for their work on Finnish sample collection, preparation and data handling, and E. Szathmári, J.B.Kovács, M. Lörincz and A. Nagy for their work with the Hungarian families. Health2000 organization, Finrisk consortium, K. Mustalahti, M. Perola, K. Kristiansson and J. Koskinen are thanked for providing the Finnish control genotypes. We thank D.G. Clayton and N. Walker for providing T1DGC data in the required format. We thank the Irish Transfusion Service and Trinity College Dublin Biobank for control samples; V. Trimble, E. Close, G. Lawlor, A. Ryan, M. Abuzakouk, C. O’Morain, G. Horgan, for celiac sample collection and preparation We acknowledge DNA provided by Mayo Clinic Rochester, and Prof. M. Bonamico and Prof. M. Barbato (Department of Paediatrics, Sapienza University of Rome, Italy) for recruiting individuals. We thank Polish clinicians for recruitment of celiac disease individuals (Z. Domagala, A. Szaflarska-Poplawska, B. Oralewska, W. Cichy, B. Korczowski, K. Fryderek, E. Hapyn, K. Karczewska, A. Zalewska, I. Sakowska-Maliszewska, R. Mozrzymas, A. Zabka, M. Kolasa, B. Iwanczak). We thank M. Szperl for isolating DNA from blood samples provided by Children’s Memorial Health Institute (Warsaw, Poland).
Dutch and UK genotyping for the second celiac disease GWAS was funded by the Wellcome Trust (084743 to D.A.vH.). Italian genotyping for the second celiac disease GWAS was funded by the Coeliac Disease Consortium, an Innovative Cluster approved by the Netherlands Genomics Initiative and partially funded by the Dutch Government (BSIK03009 to C.W.) and by the Netherlands Organisation for Scientific Research (NWO, VICI grant 918.66.620 to C.W.). E.G. is funded by the Italian Ministry of Healthy (grant RC2009). L.H.v.d.B. acknowledges funding from the Prinses Beatrix Fonds, the Adessium foundation, and the Amyotrophic Lateral Sclerosis Association. L.F. received a Horizon Breakthrough grant from the Netherlands Genomics Initiative (93519031) and a VENI grant from NWO (ZonMW grant 916.10.135). P.C.D. is a MRC Clinical Training Fellow (G0700545). G.T. received a Ter Meulen Fund grant from the Royal Netherlands Academy of Arts and Sciences (KNAW). The gene expression study was funded in part by COPACETIC (EU grant 201379). This study makes use of data generated by the Wellcome Trust Case-Control Consortium 2 (WTCCC2). A full list of the WTCCC2 investigators who contributed to the generation of the data is available from www.wtccc.org.uk. Funding for the WTCCC2 project was provided by the Wellcome Trust under award 085475. This research utilizes resources provided by the Type 1 Diabetes Genetics Consortium, a collaborative clinical study sponsored by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK), National Institute of Allergy and Infectious Diseases (NIAID), National Human Genome Research Institute (NHGRI), National Institute of Child Health and Human Development (NICHD), and Juvenile Diabetes Research Foundation International (JDRF) and supported by U01 DK062418. We acknowledge the use of BRC Core Facilities provided by the financial support from the Department of Health via the National Institute for Health Research (NIHR) comprehensive Biomedical Research Centre award to Guy’s & St Thomas’ NHS Foundation Trust in partnership with King’s College London and King’s College Hospital NHS Foundation Trust. We acknowledge funding from NIH: DK050678 and DK081645 (to S.L.N.), NS058980 (to R.A.O.); DK57892 and DK071003 (to J.A.M.). Collection of Finnish and Hungarian patients was funded by EU Commission (MEXT-CT-2005-025270), the Academy of Finland, Hungarian Scientific Research Fund (contract OTKA 61868), the University of Helsinki Funds, the Competitive Research Funding of the Tampere University Hospital, the Foundation of Pediatric Research, Sigrid Juselius Foundation, and the Hungarian Academy of Sciences (2006TKI247 to RA). Funding for the Polish samples collection and genotyping was provided by UMC Cooperation Project (6/06/2006/NDON). R.McM is funded by Science Foundation Ireland. C. Núñez has a FIS contract (CP08/0213). The Dublin Centre for Clinical Research contributed to patient sample collection and is funded by the Irish Health Research Board and the Wellcome Trust
Finally we thank all individuals with celiac disease and control individuals for participating in this study.
Written informed consent was obtained from all subjects, with Ethics Committee / Institutional Review Board approval. All individuals are of European ancestry. Affected celiac individuals were diagnosed according to standard clinical, serological and histopathological criteria including small intestinal biopsy. DNA samples were from blood, lymphoblastoid cell lines or saliva. A more detailed description of subjects is provided in Supplementary Information.
See Table 1. UK(1) case and control genotyping was previously described1,7. Illumina 670-Quad and 1.2M-Duo (custom chips designed for the WTCCC2 and comprising Hap550 / 1M and common CNV content) and 610-Quad genotyping was performed in London, Hinxton and Groningen. Bead intensity data was normalized for each sample in BeadStudio, R and theta values exported and genotype calling performed using a custom algorithm1,35. A detailed description of genotype calling steps is provided in Supplementary Information.
Quality control steps were performed in the following order: Very low call rate samples and SNPs were first excluded. SNPs were excluded from all sample collections if any collection showed call rates were <95% or deviation from Hardy-Weinberg equilibrium (P<0.0001) in controls. Samples were excluded for call rate <98%, incompatible recorded gender and genotype inferred gender, ethnic outliers (identified by multi-dimensional scaling plots of samples merged with HapMap Phase II data), duplicates and first degree relatives. 22 of 417 SNPs showing apparent association (PGWAS <10−4) were excluded after visual inspection of R theta plots suggested possible bias.
The over-dispersion factor of association test statistics (genomic control inflation factor), λGC, was calculated using observed versus expected values for all SNPs in PLINK.
See Table 1. Finnish controls (12) were genotyped on the 610-Quad BeadChip, other samples were genotyped using Illumina GoldenGate BeadXpress assays in London, Hinxton and Groningen. Genotyping calling was performed in BeadStudio for combined cases and controls in each separate collection, with the exception of the Finnish collection, and whole genome amplified samples (89 Irish cases and 106 Spanish controls). Quality control steps were performed as for the GWAS. 131 of 144 SNPs passed quality control and visual inspection of genotype clouds.
Analyses were performed using PLINK v1.0736, mostly using the Cochran-Mantel-Haenzel test. Logistic regression analyses were used to define the independence of association signals within the same linkage disequilibrium block, with group membership included as a factorized covariate
Genotype imputation was performed for samples genotyped on the Hap300 using BEAGLE and CEU, TSI, MEX and GIH reference samples from HapMap3. Association analysis was performed using logistic regression on posterior genotype probabilities, with group membership included as a factorized covariate.
Structured association tests were performed using PLINK, as described using genetically matched cases and controls within collections identified by identity by state similarity across autosomal non-HLA SNPs as described34 (settings --ppc 0.001 --cc, clusters constrained by the 5 collections). Principal components analysis was performed using EIGENSTRAT and a set of 12,810 autosomal non-HLA SNPs chosen for low LD and ancestry information37,38, association tests were corrected for the top 10 principal components and combined using weighted Z scores.
The fraction of additive variance was calculated using a liability threshold model39 assuming a population prevalence of 1%. Effect sizes and control allele frequencies were estimated from the combined replication panel. Genetic variance was calculated assuming 50% heritability.
We performed GRAIL analysis (http://www.broadinstitute.org/mpg/grail/grail.php) using HG18 and Dec2006 PubMed datasets, default settings for SNP rs number submission, and the 27 genome-wide significant celiac disease risk loci (most associated SNP) as seeds. As a query we used either associated SNPs, or 101 × 50 randomly chosen Hap550 SNP datasets (5050 SNPs, of which 5033 mapped to the GRAIL database).
We noted that the power of eQTL studies in humans is limited by substantial observed inter-individual variation in expression measurements due to non-genetic factors, and therefore developed a method, Transcriptional Components, to remove a large component of this variation (manuscript in preparation). Expression data from 42,349 heterogeneous human samples hybridized to Affymetrix HG-U133A (GEO accession number: GPL96) or HG-U133 Plus 2.0 (GEO accession number: GPL570) Genechips were downloaded 40 (Fig. 1, step 1). Samples missing data for >150 probes were excluded, and only probes available on both platforms were analysed, resulting in expression data for 22,106 probes and 41,408 samples. We performed quantile normalization using the median rank distribution41 and log2 transformed the data - ensuring an identical distribution of expression signals for every sample, discarding previous normalization and transformation steps.
Initial quality control (QC) was performed by applying principal component analysis (PCA) on the sample correlation matrix (pair-wise Pearson correlation coefficients between all samples). The first principal component (PC), explaining ~80-90% of the total variance42,43, describes probe-specific variance. 6,375 samples with correlation R < 0.75 of the sample array with this PC were considered outliers of lesser quality and excluded from analysis. We excluded entire GEO datasets where >25% of the samples were outliers (probably expression ratios versus a reference, not absolute data). The final dataset comprised 33,109 samples (17,568 GPL96 and 15,541 GPL570 samples), and we repeated the normalization and transformation on the originally deposited expression values of these post-quality control samples.
We next applied PCA on the pairwise 22,106 × 22,106 probe Pearson correlation coefficient matrix assayed on the 33,109 sample dataset (our fast C++ tool, MATool, is available upon request), attempting to simplify the structure of the data. Here, PCA represents a transformation of a set of correlated probes into sets of uncorrelated linear additions of probe expression signals (eigenvectors) that we name Transcriptional Components (TCs). Each TC is a weighted sum of probe expression signals and eigenvector probe coefficients. These TC-scores can be calculated for each observed expression array sample (reflecting the TC activity per sample).
We obtained peripheral blood DNA and RNA (PAXgene) from Dutch and UK individuals who were disease cases or controls for GWAS studies (Supplementary Table 1). All samples had been genotyped for a common SNP set on Illumina platforms. Analysis was confined to 294,767 SNPs that had a MAF ≥ 5%, call-rate ≥ 95% and exact HWE P>0.001. RNA from the samples was either hybridized to Illumina HumanRef-8 v2 arrays (229 samples, Ref-8v2) or Illumina HumanHT-12 arrays (1,240 samples, HT-12), and raw probe intensity extracted using BeadStudio. The Ref-8v2 samples were jointly quantile normalized and log2 transformed, and similarly for the HT-12 samples. Subsequent analyses were also conducted separately for these datasets, up to the eventual eQTL mapping, that uses a meta-analysis framework, combining eQTL results from both arrays. HT-12 and Ref-8v2 arrays are different, but share many probes with identical probe sequences. Illumina sometimes use different probe identifiers for the same probe sequences - in meta-analysis and Table 3, the label HT-12 was used if both HT-12 and Ref-8v2 had the same sequence.
If probes mapped incorrectly, or cross-hybridized to multiple genomic loci, it might be that an eQTL would be detected that would be deemed a trans-eQTLs. To prevent this, we used a mapping approach versus a known reference that we developed for high-throughput short sequence RNAseq data44. We took the DNA sequence as synthesised for each cDNA probe, and aligned it versus a transcript masked gDNA genome combined with cDNA sequences. A more detailed description of probe re-mapping is provided in Supplementary Information. Probes that did not map, or mapped to multiple different locations were removed.
TC-scores can be inferred in new (non-Affymetrix) datasets for every new individual sample. For the Illumina samples (used for the cis-eQTL mapping), only Illumina probes that could be mapped to any of our 22,106 Affymetrix probes were used (www.switchtoi.com/probemapping.ilmn). The TC-score of sample i for the jth TC is defined as: , where vtj is defined as the tth Affymetrix probe coefficient for the jth TC; ati is the Illumina expression measurement for the tth mapped probe for sample i. We inferred the Illumina TC-scores for the top 1,000 TCs.
Because our Illumina eQTL dataset (n = 1,469) is much less heterogeneous than the Affymetrix dataset (n = 33,109), we expect that some TCs will hardly vary. We therefore performed a PCA on the covariance matrix of the top 1,000 inferred TC-scores for the Illumina dataset to effectively compress the TC data into a small set of ‘aggregate TCs’ (aTCs). As aTCs are orthogonal, we used linear regression to eliminate the effect of the top 50 aTCs. We correlated the TC-scores for each peripheral blood sample with probe expression levels. We then used the resulting residual gene expression data for subsequent cis-eQTL mapping.
We used the residual gene expression data (Fig. 1) in a meta-analysis framework, as described45,46. In brief, analyses were confined to those probe-SNP pairs for which the distance from probe transcript midpoint to SNP genomic location was less than 500 kb. To prevent spurious associations due to outliers, a non-parametric Spearman’s rank correlation analysis was performed. When a particular probe-SNP pair was present in both the HT12 and H8v2 datasets, an overall, joint p-value was calculated using a weighted Z-method (square root of the dataset sample number). To correct for multiple testing we controlled the false discovery rate (FDR). The distribution of observed P values was used to calculate the FDR, by permuting expression phenotypes relative to genotypes 1000 times within the HT12 and H8v2 dataset. Finally, we removed any probes from analysis which contain a known SNP (1000Genomes CEU SNP data, April 2009 release).
The authors declare no competing financial interests.
DATABASE ACCESSION NUMBERS
Expression data is available in GEO (http://www.ncbi.nlm.nih.gov/geo/) as GSE11501 and GSE20142.