Cindy Boer

Cartilage Thickness and Hip Osteoarthritis | 87 2.1 croarrays were hybridized. Sample pairs were randomly dispersed over the microar- rays, however each pair was measured on a single chip. Microarrays were read using an Illumina Beadarray 500GX scanner and after basic quality checks using Beadstudio software data were analyzed in R statistical programming language. Intensity values were normalized using the “rsn” option in the Lumi-package and absence of large scale between-chip effects was confirmed using the Globaltest-package in which the individ- ual chip numbers were tested for association to the raw data. After removal of probes that were not optimally measured (detection p-value>0.05 in more than 50% of the samples) a paired t-test was performed on all sample pairs while adjusting for chip (to adjust for possible batch effects) and using multiple testing correction as implement- ed in the “BH” (Benjamini and Hochberg) option in the Limma-package. Analyses for differential expression between OA and healthy and between preserved and healthy cartilage was performed likewise, adjusting in addition for sex and for age. Exome sequencing Exome sequencing was performed in 2628 individuals from the Rotterdam Study the average mean coverage was 55x, corresponding to approximately 80% of the targeted regions covered by at least 20 reads. The exome sequencing was performed in house (HuGe-F, www.glimDNA.org ). Details of the technical procedure and variant calling are given in S2 Text. We tested the exome variants for association with mJSW and/or hip OA in two ways. Each individual variant was tested for association with mJSW using the single variant option within RV-test, while adjusting for age and sex. In addition, we did a burden test for each of the selected genes by using SNP-set kernel association test (SKAT-O). SKAT aggregates individual score test statistics of SNPs in a SNP set and computes SNP-set level p-values for a gene[23]. Visualization of the regulatory landscape of mJSW associated loci For each of the top mJSW GWAS associated SNPs the LD region was determined using the 1000G Phase 1 population using the Haploreg tool [47]. The LD threshold was set at r 2 ≥0.8. For each of these SNPs it was determined if the variant was located in a potential enhancer region using the Roadmap consortium reference epigenomes data set[27]. Heatmaps were constructed by calculating the percentage of variants in LD with the top mJSW GWAS found SNP located in enhancer regions as defined by the Roadmap epigenome chromatin states. The reference epigenomes were downloaded from the of- ficial data portal accompanying[27]. Reference epigenome data was used from mesen- chymal stem cell derived chondrocyte cultured cells, Osteoblast, Bone marrow derived cultured mesenchymal stem cells, K562, HUVEC, HeLA and NHEK cells. Reference epig- enomes were chromatin state models based on ChIPseq data of 5 core histone marks