the genome wide association study was conducted in a discovery sample of 750 patients of predominantly british descent with clinical and radiological evidence of pdb in whom mutations of the sqstm1 gene had been excluded by dna sequencing . these comprised subjects who had participated in the prism study ( n=597 ) a randomised trial of two different treatment strategies for pdb38 ; clinic - based patients from the uk with sporadic pdb ( n=55 ) ; and patients with a family history of pdb derived from the uk ( n=20 ) , australia ( n=66 ) , new zealand ( n=8 ) and italy ( n=4 ) . details of the 1002 control subjects have previously been described 9 , but in brief , they comprised healthy subjects of scottish descent with no clinical evidence of pdb . for the replication study we conducted genotyping in an additional 500 pdb patients without sqstm1 mutations diagnosed according to standard techniques . these comprised subjects with sporadic pdb who had been recruited from hospital clinics in the uk ( n=226 ) , italy ( n=20 ) and spain ( n=200 ) ; patients with sporadic pdb who had participated in the prism study ( n=43 ) and subjects with a positive family history of pdb who had been recruited from hospital clinics in australia ( n=10 ) , and the uk ( n=1 ) . the 535 replication controls comprised subjects from the uk who had been referred for investigation of osteoporosis but who had been found to have normal bone density on examination by dual energy x - ray absorptiometry ( n=248 ) ; spouses of participants of the prism study who were not known to be affected by pdb ( n=252 ) and clinic based controls from spain ( n=35 ) . the discovery sample had 96 % power to detect disease associated allele with maf = 0.2 and genotype relative risk of 1.6 assuming a multiplicative model and a disease with population prevalence of 2% . genotyping of pdb cases was performed at the genetics core of the wellcome trust clinical research facility ( edinburgh , uk ) using illumina humanhap300 beadchip version 2 duo . genotypes for cases and controls were called using beadstudio v3.2 ( illumina , inc . ) by following the manufacturer s recommended protocol . genotype data for control subjects were provided after applying the quality control ( qc ) measures described previously9 . for the cases , we used a no - call threshold of 0.15 in beadstudio and quality control metrics such as cluster separation , ab t mean and ab r mean to exclude badly performing snps . samples with a call rate of less than 90% were excluded ( n = 30 ) . the data were then subjected to further qc measures using plink39 to exclude snps with a call rate of less than 95% , those with hardy - weinberg equilibrium ( hwe ) p - values of less than 1.0  10 in controls , and those with a minor allele frequency of less than 1% . this left a total of 294,663 snps common to cases and controls with at least 95% call rates in each set . samples with excess heterozygosity ( 1 case ) , non - european ancestry ( 21 cases and 1 control ) , and related subjects ( 6 cases ) were excluded before analysis leaving a final total of 692 cases and 1001 controls with an average (  sd ) genotype call rate of 99.63  1.0 . the genotype cluster plots for all snps showing association with pdb at p < 1.0  10 were visually inspected in beadstudio . population ancestry was determined using multidimensional scaling analysis of identity - by - state ( ibs ) distances matrix of all individuals after combining genotype data from the hapmap project ( release 22 ) samples of european ( ceu ) , asian ( chb and jpt ) , and african ( yri ) ancestry . for this analysis we first removed snps in areas of extended linkage disequilibrium ( ld ) ( chr2 : 134.0 -138.0 , chr6 : 25.0 - 34.0 , chr8 : 8.0  12.0 , chr11 : 45.0  57.0)40 , and those with r > 0.2 within a 150 snp window . snps with call rate < 99% , maf < 5% , and hwe p < 1.0  10 in cases or controls were also excluded leaving a total of 63,528 snps . the genome wide average ibs distances matrix for all pair of individuals was then calculated based on the 63,528 snps using plink and then used for multidimensional scaling analysis . supplementary fig 1 is a plot of the first two component of the multidimensional scaling analysis showing three clusters corresponding to ceu , chb+jpt , and yri samples with the majority of cases and controls located within the european ceu cluster . we identified 21 cases and 1 control as outliers from the ceu cluster and these were excluded from further analysis . based on genome wide ibs distance , we identified 5 identical pairs ( ibs distance > 99% ) and 1 related pair ( ibs distance > 85% ) of samples from the cases cohort ; the sample with the lowest call rate was excluded from each pair before further analysis . genotyping of replication samples was performed by sequenom ( hamburg , germany ) using the massarray iplex platform . dna from cases and controls were distributed into 384 well plates so that each plate had the same number of cases and controls to minimize genotyping bias due variations between runs . replication samples with call rate < 95% were excluded ( 19 cases and 15 controls ) leaving a total of 481 cases and 520 controls with an average genotype call rate of 99.61% . genotypes were imputed for untyped variants located within 2.0 mb of snps identified in stage 1 with genome wide significant association with pdb using mach41 . the hapmap european ( ceu ) genotype data from release 22 were used as a reference . to minimize end effects , only data from the middle 2 mb of each 4 mb imputed segment was used for further analysis . we used 200 rounds of markov chain iterations to estimate allele dosage and best guess genotypes in stage 1 data . imputation quality was assessed by estimating the correlation ( r ) between imputed and true genotypes . analysis of imputed data was performed using logistic regression implemented in mach2dat42 in which the imputed allelic dosage was used to account for uncertainty in imputed genotypes . statistical analyses were performed using plink ( version 1.07 ) 39 . in stage 1 , genotyped snps were tested for association with pdb using a stratified cochran - mantel - haenszel ( cmh ) test . samples were stratified based on their genome - wide ibs similarity so that individuals assigned to one cluster were not genetically different ( p > 0.001 from pair wise population concordance test ) . the quantile - quantile plot and genomic control  were used to assess overdispersion of the test statistics and were calculated using the statistical package r version 2.7.2 based on the 90% least significant snps as described previously10 . stepwise logistic regression was used to test for independent effects of snp where the allelic dosage of the conditioning snp was entered as a covariate in the regression model along with the population clusters identified by ibs sharing analysis described above to adjust for population substructure . haplotype analysis was performed by logistic regression looking at the presence or absence of the test haplotype and including the population clusters as a covariate in the model . haplotypes were phased using the e - m algorithm implemented in plink and only haplotypes with a frequency of  1% were analyzed . the cut off point for genome wide significance was set as p < 1.7  10 ( 0.05/294,663 ) for stage 1 , and p < 3  10 ( 0.05/16 ) for the replication stage . for the combined analysis we set the threshold for significance as p < 5  10 as recently proposed43 . the replication and combined data sets were analysed as described above except that the replication dataset was considered as a separate cluster when population clusters were used in a stratified cmh test or as a covariate in logistic regression models . the population attributable risk ( par ) for markers showing association with pdb was calculated according to the following formula : par = p(or-1)/[p(or-1)+1 ] ; where p is the frequency of the risk allele in controls and or is the risk allele odds ratio . the cumulative par was calculated as follows : cumulative par = 1- ( 1n ( 1-pari ) ) ; where n is the number of variants and pari is the individual par for the ith snp . the genome wide association study was conducted in a discovery sample of 750 patients of predominantly british descent with clinical and radiological evidence of pdb in whom mutations of the sqstm1 gene had been excluded by dna sequencing . these comprised subjects who had participated in the prism study ( n=597 ) a randomised trial of two different treatment strategies for pdb38 ; clinic - based patients from the uk with sporadic pdb ( n=55 ) ; and patients with a family history of pdb derived from the uk ( n=20 ) , australia ( n=66 ) , new zealand ( n=8 ) and italy ( n=4 ) . details of the 1002 control subjects have previously been described 9 , but in brief , they comprised healthy subjects of scottish descent with no clinical evidence of pdb . for the replication study we conducted genotyping in an additional 500 pdb patients without sqstm1 mutations diagnosed according to standard techniques . these comprised subjects with sporadic pdb who had been recruited from hospital clinics in the uk ( n=226 ) , italy ( n=20 ) and spain ( n=200 ) ; patients with sporadic pdb who had participated in the prism study ( n=43 ) and subjects with a positive family history of pdb who had been recruited from hospital clinics in australia ( n=10 ) , and the uk ( n=1 ) . the 535 replication controls comprised subjects from the uk who had been referred for investigation of osteoporosis but who had been found to have normal bone density on examination by dual energy x - ray absorptiometry ( n=248 ) ; spouses of participants of the prism study who were not known to be affected by pdb ( n=252 ) and clinic based controls from spain ( n=35 ) . the discovery sample had 96 % power to detect disease associated allele with maf = 0.2 and genotype relative risk of 1.6 assuming a multiplicative model and a disease with population prevalence of 2% . genotyping of pdb cases was performed at the genetics core of the wellcome trust clinical research facility ( edinburgh , uk ) using illumina humanhap300 beadchip version 2 duo . genotypes for cases and controls were called using beadstudio v3.2 ( illumina , inc . ) by following the manufacturer s recommended protocol . genotype data for control subjects were provided after applying the quality control ( qc ) measures described previously9 . for the cases , we used a no - call threshold of 0.15 in beadstudio and quality control metrics such as cluster separation , ab t mean and ab r mean to exclude badly performing snps . samples with a call rate of less than 90% were excluded ( n = 30 ) . the data were then subjected to further qc measures using plink39 to exclude snps with a call rate of less than 95% , those with hardy - weinberg equilibrium ( hwe ) p - values of less than 1.0  10 in controls , and those with a minor allele frequency of less than 1% . this left a total of 294,663 snps common to cases and controls with at least 95% call rates in each set . samples with excess heterozygosity ( 1 case ) , non - european ancestry ( 21 cases and 1 control ) , and related subjects ( 6 cases ) were excluded before analysis leaving a final total of 692 cases and 1001 controls with an average (  sd ) genotype call rate of 99.63  1.0 . the genotype cluster plots for all snps showing association with pdb at p < 1.0  10 were visually inspected in beadstudio . population ancestry was determined using multidimensional scaling analysis of identity - by - state ( ibs ) distances matrix of all individuals after combining genotype data from the hapmap project ( release 22 ) samples of european ( ceu ) , asian ( chb and jpt ) , and african ( yri ) ancestry . for this analysis we first removed snps in areas of extended linkage disequilibrium ( ld ) ( chr2 : 134.0 -138.0 , chr6 : 25.0 - 34.0 , chr8 : 8.0  12.0 , chr11 : 45.0  57.0)40 , and those with r > 0.2 within a 150 snp window . snps with call rate < 99% , maf < 5% , and hwe p < 1.0  10 in cases or controls were also excluded leaving a total of 63,528 snps . the genome wide average ibs distances matrix for all pair of individuals was then calculated based on the 63,528 snps using plink and then used for multidimensional scaling analysis . supplementary fig 1 is a plot of the first two component of the multidimensional scaling analysis showing three clusters corresponding to ceu , chb+jpt , and yri samples with the majority of cases and controls located within the european ceu cluster . we identified 21 cases and 1 control as outliers from the ceu cluster and these were excluded from further analysis . based on genome wide ibs distance , we identified 5 identical pairs ( ibs distance > 99% ) and 1 related pair ( ibs distance > 85% ) of samples from the cases cohort ; the sample with the lowest call rate was excluded from each pair before further analysis . genotyping of replication samples was performed by sequenom ( hamburg , germany ) using the massarray iplex platform . dna from cases and controls were distributed into 384 well plates so that each plate had the same number of cases and controls to minimize genotyping bias due variations between runs . replication samples with call rate < 95% were excluded ( 19 cases and 15 controls ) leaving a total of 481 cases and 520 controls with an average genotype call rate of 99.61% . genotypes were imputed for untyped variants located within 2.0 mb of snps identified in stage 1 with genome wide significant association with pdb using mach41 . the hapmap european ( ceu ) genotype data from release 22 were used as a reference . to minimize end effects , only data from the middle 2 mb of each 4 mb imputed segment was used for further analysis . we used 200 rounds of markov chain iterations to estimate allele dosage and best guess genotypes in stage 1 data . imputation quality was assessed by estimating the correlation ( r ) between imputed and true genotypes . analysis of imputed data was performed using logistic regression implemented in mach2dat42 in which the imputed allelic dosage was used to account for uncertainty in imputed genotypes . statistical analyses were performed using plink ( version 1.07 ) 39 . in stage 1 , genotyped snps were tested for association with pdb using a stratified cochran - mantel - haenszel ( cmh ) test . samples were stratified based on their genome - wide ibs similarity so that individuals assigned to one cluster were not genetically different ( p > 0.001 from pair wise population concordance test ) . the quantile - quantile plot and genomic control  were used to assess overdispersion of the test statistics and were calculated using the statistical package r version 2.7.2 based on the 90% least significant snps as described previously10 . stepwise logistic regression was used to test for independent effects of snp where the allelic dosage of the conditioning snp was entered as a covariate in the regression model along with the population clusters identified by ibs sharing analysis described above to adjust for population substructure . haplotype analysis was performed by logistic regression looking at the presence or absence of the test haplotype and including the population clusters as a covariate in the model . haplotypes were phased using the e - m algorithm implemented in plink and only haplotypes with a frequency of  1% were analyzed . the cut off point for genome wide significance was set as p < 1.7  10 ( 0.05/294,663 ) for stage 1 , and p < 3  10 ( 0.05/16 ) for the replication stage . for the combined analysis we set the threshold for significance as p < 5  10 as recently proposed43 . the replication and combined data sets were analysed as described above except that the replication dataset was considered as a separate cluster when population clusters were used in a stratified cmh test or as a covariate in logistic regression models . the population attributable risk ( par ) for markers showing association with pdb was calculated according to the following formula : par = p(or-1)/[p(or-1)+1 ] ; where p is the frequency of the risk allele in controls and or is the risk allele odds ratio . the cumulative par was calculated as follows : cumulative par = 1- ( 1n ( 1-pari ) ) ; where n is the number of variants and pari is the individual par for the ith snp .
paget s disease of bone ( pdb ) is a common disorder with a strong genetic component characterised by focal increases in bone turnover which in some cases is caused by sqstm1 mutations . to identify additional susceptibility genes we performed a genome wide association study in 750 pdb cases without sqstm1 mutations and 1002 controls and identified three candidate loci for the disease which were replicated in an independent set of 500 cases and 535 controls . 
 the strongest signal was with rs484959 on 1p13 close to the csf1 gene ( p = 5.38  1024 ) and significant associations were also observed with rs1561570 on 10p13 within the optn gene ( p = 6.09  1013 ) and with rs3018362 on 18q21 close to the tnfrsf11a gene ( p = 5.27  1013 ) . 
 these studies provide new insights into the pathogenesis of pdb and identify optn , csf1 and tnfrsf11a as novel candidate genes for disease susceptibility .