SUPPLEMENTARY METHODS
Stage 1 samples.
Cases.Within the arcOGEN consortium, five locations (London, Nottingham, Oxford, Sheffield and Southampton) had existing DNA collectionsand contributed OA DNA samples for stage 1 of the study. All cases were collected in the UK, were unrelated and of European origin. The phenotype was determined by two criteria; radiographic evidence of disease and clinical evidence of disease to a level requiring joint replacement. Each investigator defined radiographic disease as a Kellgren-Lawrence (KL) grade[1] of ≥2. In addition, each centre had some site-specific variations.
London contributed 348 cases from the TwinsUK cohort[2] and the Chingford Study[3]with radiographic OA, of which 51 also had clinical disease and had undergone a total joint replacement (20 knee, 29 hip and 2 both hip and knee). Radiographic knee OA was defined as a KL grade ≥2 for the tibio-femoral compartment and radiographic hip OA as a Croft grade ≥2.[4]
Nottingham contributed 879 cases from two sources: 321 from the “Nottingham” collection and 558 from the Generalised Osteoarthritis collection.Cases were over age 40 with clinically severe primary hip or knee OA referred to hospital for joint surgery.[5,6] All participants underwent clinical enquiry and examination to exclude other arthropathies or cause of hip disease, and radiographs were examined by the same trained observer to confirm radiographic OA (all had KL grade =>2).
TheOxford collection contributed 1540 OA cases. All patients had undergone joint replacement surgery for primary osteoarthritis of the hip or knee. They had a minimum KL grade of 3.
Sheffield contributed 240 patients, ascertained through joint replacement surgeryfor primary osteoarthritis of the hip (KL grade >2). All subjects were free from known inflammatory disorders at the time of recruitment, and were not taking any drug medications known to affect bone metabolism. The characteristics of this study population are published in detail elsewhere.[7]
The Southampton collection contributed 170 patients. All cases were subjects with symptomatic, radiographic knee OA defined by a KL score of 2 or more that were recruited from the community into a randomized, placebo controlled trial of vitamin D replacement (VIDEO Study). All participants underwent clinical enquiry and examination to exclude other arthropathies.
All stage 1 cases listed in Supplementary Table 1 have primary OA with exclusion of inflammatory arthritis (rheumatoid, polyarthritic or autoimmune disease), post-traumatic arthritis, post-septic arthritis, skeletal dysplasia and developmental dysplasia.
Controls. The main study used 4,894 population-based UK controls from an early release of the Wellcome Trust Case Control Consortium 2 data.[8] Control samples came from 2 distinct sources (the 1958 Birth Cohort and the UK Blood Donor Service) and were unrelated. For sensitivity analysis purposes, we used data from the TwinsUK cohort, which consists of a group of twins ascertained to study the heritability and genetics of age-related diseases.[9]These unselected twins were recruited from the general population through national media campaigns in the UK and shown to be comparable to age-matched population singletons in terms of disease-related and lifestyle characteristics.[10]
Stage 1 genome-wide genotyping.
DNA sample preparation and quality control.Genomic DNA from extant cohorts and blood samples from new cases were shipped by the collection centres to The Centre for Integrated Genomic Medical research (CIGMR), University of Manchester, for processing. DNA was extracted from blood samples using a PSS Magtration XL DNA Extraction System. All DNA samples and associated patient information were logged into a Laboratory Information Management System (LIMS) database. 2D bar coded matrix plates and patient data manifests for each set of 96 DNAs were shipped to The Wellcome Trust Sanger Institute, Hinxton. The quality of the DNA and subject identity were validated using the Sequenom iPLEX assay designed to genotype 4 gender-specific SNPs and 26 SNPs present on the Illumina Human 610-Quad array. DNA concentrations were quantified using a PicoGreen assay (Invitrogen) and an aliquot assayed by agarose gel electrophoresis. A DNA sample was considered to pass quality control if the original DNA concentration was >50ng/ul, the DNA was not degraded, the gender assignment from the iPLEX assay matched that provided in the patient data manifest and genotypes were obtained for over 65% of the SNPs on the iPLEX.
Genotyping and allele calling. DNA samples passing QC were rearrayed into 96-well plates for genotyping using Illumina Human 610-Quad BeadChips. Approximately 200ng of genomic DNA was used for genotyping. Briefly, each sample was whole-genome amplified, fragmented, precipitated and resuspended in appropriate hybridization buffer. Denatured samples were hybridised on prepared Illumina Human 610-Quad BeadChips overnight. After hybridisation, the BeadChip oligonucleotides were extended by addition of a single fluorescently labelled base which was detected by fluorescence imaging with an Illumina Bead Array Reader. Normalised bead intensity data obtained for each sample were converted into SNP genotypes using a customised genotype calling algorithm.[11]
Stage 1 genotype quality control.
SampleQC.Unless otherwise stated, quality control was performed using plink [12] and was carried out separately for OA cases and each control dataset.Samples were excluded if their call rate was <97%, and if they showed gender discrepancies (estimated from genotypic data against external information). Individuals were also excluded on the basis of excess genome-wide heterozygosity or homozygosity. For each sample group, heterozygosity histograms were inspected to determine exclusion thresholds empirically. For the cases, these were set to >35% and <28% for excess heterozygosity and excess homozygosity respectively. These thresholds were set to >36% and <32% for both of the control groups. We identified samples that were accidentally duplicated or closely-related by calculating genome-wide IBD (given IBS information) for pairs of individuals. Multidimensional scaling (MDS) was performed in conjunction with data from the three HapMap phase II populations in order to identify and exclude individuals of non-European descent (Supplementary Figure 5). In addition, we used GoldSurfer2 to carry out principal component analysis (PCA).[13] Of the 3324 OA cases, 2587 UKBS controls and 2482 1958C controls present in the original dataset, these sample quality control filters resulted in the inclusion of 3177 cases, 2474 UKBS and 2420 1958BC controls in downstream analyses.
SNP QC. SNPs were excluded from further analysis based on the following criteria: Non autosomal, call rate <95% if minor allele frequency (MAF)≥5% or call rate <99% if MAF<5%, HWE exact p values <0.0001 and MAF <1% (MAF exclusion doesn’t apply for the rare variant analysis). GC/AT SNPs were also removed. Of the 582,585 autosomal SNPs genotyped in the OA cases using the Illumina Human610 chip, 537,772 passed these filters. 517,868 of these overlapped with SNPs genotyped on the Illumina Human1M chip in the two control groups. SNP and sample quality control of the control groups was performed after extracting the 517,868 overlapping clean SNPs. Furthermore, 43 SNPs with multiplicative p value <5.7x10-7 from the association test between the two control groups were removed before pooling them together as one control group. This resulted in 514,964 SNPs to be used for case/control analysis. Subsequently, intensity plots were inspected for SNPs with multiplicative p value <10-4 from the case/control analysis and 66 SNPs were further removed, resulting in a final dataset comprising 514,898 SNPs.
Association analysis of stage 1 genotype data. 514,898 autosomal SNPs were analysed under the multiplicative (or log additive) model using plink.[12] The genomic control (GC) value for directly-typed SNPs was estimated to be 1.077; after correcting for population structure, by including the first ancestry informative principal component (PC)[14] as a covariate, the GC value was estimated to be 1.078 (Supplementary Results). We also carried out stratified analysis by joint (hip, knee) and gender.
Prioritisation of SNPs for follow-up. We prioritised 102 SNPs (Supplementary Table 2) for follow-up by in silico replication by the deCODE, Rotterdam and Framingham studies based on a set of heuristic criteria we developed. Briefly, we focused on signals with multiplicative model p values<0.0001 in the OA v. controls and joint-stratified cases v. controls analyses. We defined independent SNPs on the basis of r2<0.4. We investigated QC properties (call rate, exact HWE p value in cases and controls, differential missingness between cases and controls) for these independent SNPs with p<0.0001 and examined cluster plots for the cases and two sets of controls separately. We preferentially selected SNPs with frequencies over 0.10 to boost power of replication studies. We additionally examined the genic location of variants and up-weighted signals within or near strong candidate genes for OA. We used GeneSniffer[15] to assign candidacy scores to genes within 0.1centiMorgans either side of each index SNP. We selected a subset of 52 SNPs for follow-up by in silico replication by the TwinsUK study primarily by ranking the selected 102 signals on the basis of statistical significance.
After obtaining in silico replication data from the deCODE, Rotterdam and TwinsUK studies, we carried out a meta-analysis across all available data and prioritised 36 of the 102 SNPs for de novo replication mainly on the basis of their overall statistical significance.
Stage 1 genotype imputation and analysis of imputed data. We imputed genotypes for autosomal SNPs that were present in HapMap Phase II but were not present in the genome-wide chip or did not pass direct genotyping QC. In each sample, genotypes were imputed using the directly typed data and phased HapMap II genotype data from the 60 CEU HapMap founders.[16] Genotypes were imputed using the program IMPUTE,[17] which determines the probability distribution of missing genotypes conditional on a set of known haplotypes and an estimated fine-scale recombination map. Imputation was based on 514,865 autosomal SNPs with MAF>0.01 (excluding SNPs that demonstrated poor genotype clustering upon manual inspection). We analysed 2,018,811 imputed SNPs (excluding the directly typed SNPs passing QC) that had MAF>0.01 in cases and controls using SNPTEST, appropriately accounting for the probability distributions of imputed genotypes.[18] We included 1,829,948 imputed autosomal SNPs with MAF>0.01 in cases and controls and an imputation information score >0.8 in downstream association analyses of all OA cases v. controls, 1,830,059 SNPs in hip OA v. controls, 1,829,812 SNPs in knee OA v. controls, 1,829,458 SNPs in female OA v. female controls, 1,829,623 SNPs in female hip OA v. female controls, 1,829,251 SNPs in female knee OA v. female controls, 1,829,817 SNPs in male OA v. male controls, 1,829,690 SNPs in male hip OA v. male controls, and 1,829,567 SNPs in male knee OA v. male controls. The genomic control values for imputed SNPs in each of the analyses above were 1.079, 1.056, 1.081, 1.060, 1.037, 1.057, 1.035, 1.021, and 1.036 respectively.
In silico replication samples, genotyping and analysis.Four OA cohorts that had previously been subjected to a GWAS were used for in silico replication of promising signals from the arcOGEN GWAS in the first instance: the Rotterdam study, deCODE, Framingham and TwinsUK (Supplementary Table 1). Summary association statistics (under the log-additive model) were shared for 102 SNPs (Rotterdam, deCODE and Framingham) and 52 SNPs (TwinsUK) respectively. The results of these analyses guided the selection of 36 SNPs for further de novo and in silico (in the Estonian Genome Center, University of Tartu dataset) follow-up.
The Rotterdam Study. The study population comprises men and women aged 55 years and older from the Rotterdam Study, which is a prospective population-based study on determinants of chronic disabling diseases. The rationale and study design have been described previously.[19] The medical ethics committee of Erasmus University Medical School approved the study and written informed consent was obtained from each participant. Subjects were scored for the presence of OA using standardized radiographs of the hip and knee with cases having a Kellgren and Lawrence (KL) grade of ≥2and controls a KL grade of <2. Genotyping of the samples with the HumanHap550v3 Genotyping BeadChip (Illumina, San Diego, USA) was carried out at the Genetic Laboratory of the Department of Internal Medicine of Erasmus Medical Center, Rotterdam, the Netherlands. The BeadStudio GenCall algorithm was used for genotype calling and quality control procedures were as described previously.[20]The following sample exclusion criteria were applied: call rate <97.5%; gender mismatches with typed X-linked markers; autosomal heterozygosity >0.336~FDR>0.1%, duplicates and/or 1st or 2nd degree relatives using IBS probabilities >97% from plink; ethnic outliers using IBS distances >3SD from plink.[12] Analysis was restricted to SNPs with a call rate ≥98%, minor allele frequency (MAF) ≥0.01, andHWE p value ≥1x10-6. MACH software was used for imputation and statistical analysis was carried out using MACH2QTL.[21]
deCODE. The study population comprises patients with OA of the knee and/or hip obtained on the basis of patients’ records at hospitals and health care centres in Iceland.[22]All OA patients had undergone joint replacement surgery of the knee or the hip. A clinician reviewed the patient records to verify the diagnosis.Control individuals were recruited as part of various genetic programs at deCODE. All individuals that had chip genotype data that passed quality control were included in the control group with the exception of those on any OA listever, who were excluded. The study was approved by the Data Protection Authority of Iceland and the National Bioethics Committee of Iceland. Informed consent was obtained from all participants. The samples were assayed with the Infinium HumanHap300 or humanCNV370 SNP chips (Illumina) and called with BeadStudio. All of the genotyped SNPs tested in this report passed quality filtering (call rate >96%, MAF>0.01, HWE p>10-6 on any of the three chip types used [HumanHap300, HumanHap300-duo and HumanCNV370]). SNPs were imputed using the IMPUTE software[17] and phased haplotypes for the HapMap CEU sample set v22. The difference in frequency of each SNP was tested, assuming an additive model, using the SNPTEST program[17] with the –proper option which uses a likelihood ratio test to incorporate the uncertainty in imputed genotypes. Any samples with a yield <98% were excluded from the analysis.
The Framingham osteoarthritis study. This is a longitudinal population-based cohort study established in 1948 in Framingham, Massachusetts to examine risk factors for heart disease.[23] In addition to the original cohort, a study of the offspring and their spouses of this cohort was initiated in 1971. The Framingham OA study, which includes participants of both cohorts, was developed to study the inheritance of OA.[24] Knee OA cases and controls were available for association analyses, with cases having a KL grade of ≥2 in at least one knee and controls having KL grades of 0 or 1 in both knees. The Boston University Medical Center IRB approved the Osteoarthritis Protocol. Written informed consent was obtained from all subjects. Samples were genotyped using the Affymetrix GeneChip® Human Mapping 500K array set and the 50K supplemental array set focused on coding SNPs and SNPs tagging protein-coding genes (Santa Clara, California) as part of the SHARe initiative. Sample exclusion criteria included call rate <97% and a per subject heterozygosity >±5 standard deviations from the mean, while any samples that demonstrated excessive Mendelian errors were excluded. Imputation was performed using MACH (version 1.0.15)[21] to impute all autosomal SNPs using the publicly available phased haplotypes from HapMap (release 22, build 26, CEU population) as a reference population. A sample of 200 known unrelated participants with high call rates and low Mendelian errors were used to determine parameter estimates that were subsequently applied in a model for all subjects. Dosage estimates outputted by MACH were used in analysis.
TwinsUK. The study participants were white monozygotic and dizygotic twin pairs from the TwinsUK adult twin registry, a group used to study the heritability and genetics of age-related diseases.[25] These unselected twins were recruited from the general population through national media campaigns in the UK. In this study, cases had a KL grade of ≥2 at the hip or the knee, whereas controls had a KL grade of <2 at the hip or the knee. Ethics approval was obtained from the Guy’s and St. Thomas’ Hospital Ethics Committee. Written informed consent was obtained from every participant. Samples were genotyped with the Infinium HumanHap300 assay (Illumina) at the Duke University Genotyping Center (NC USA), Helsinki University (Finland) and the Wellcome Trust Sanger Institute. The Illuminus calling algorithm was used for genotype calling. Analysis was restricted to SNPs with a call rate >90%, MAF ≥0.01, and HWE p value ≥1x10-4[20] whilst imputation was performed using IMPUTE[17]. At imputed loci, all genotypes with posterior probabilities <0.9 were discarded and the imputed loci were filtered out using usual QC filters.
Estonian Genome Center, University of Tartu. The Estonian cohort is from the population-based biobank of the Estonian Genome Project of University of Tartu(EGCUT). The whole project is conducted according to the Estonian Gene Research Act and all participants have signed the broad informed consent.[26-27] The current cohort size is over 47,000, from 18 years of age and up, which reflects closely the age distribution in the adult Estonian population. Subjects are recruited by the general practitioners (GP) and physicians in the hospitals were randomly selected from individuals visiting GP offices or hospitals. Each participant filled out a Computer Assisted Personal interview during 1-2 hours at a doctor’s office, including personal data (place of birth, place(s) of living, nationality etc.), genealogical data (family history, three generations), educational and occupational history and lifestyle data (physical activity, dietary habits, smoking, alcohol consumption, women’s health, quality of life).Anthropometric and physiological measurements were also taken. Osteoarthritis was diagnosed by a specialist as a clinical finding and was usually confirmed by a radiograph (KL score>2). The OA cases for the current study had an ICD10 M16 and/or M17 diagnosis. All diseases are defined according to the ICD10 coding.[28] All the samples are genotyped with Illumina HumanCNV370 or HumanOmniExpress according to the Illumina protocol[29] in the Estonian Biocenter Genotyping Core Facility. Data quality control was performed with plink[12] (SNP call rate98%; sample call rate>95%; MAF >0.01; HWE p>10-6; cryptic relatedness). Imputation was performed with IMPUTE v1.0 (CEU HapMap rel22 build 36) and association analyses were carried out with SNPTEST.[17] Inflation factors for directly genotyped and imputed data were 1.02 and 1.01respectively.