BACKGROUND
killer immunoglobulin-like receptors  are transmembrane glycoproteins expressed by natural killer cells and subsets of t cells.
the kir region lies in a  <dig> kb gene cluster located within the  <dig> mb leukocyte receptor complex on chr19q <dig> .
the region exhibits great haplotype and copy number diversity, which has prevented complete assessment of the kir genes in standard genome-wide association studies , despite their strong candidacy for immune-related traits.
targeted quantitative polymerase chain reaction  assays have been used to detect presence or absence of individual kir genes and more recently, determine copy numbers  <cit> .
nevertheless these remain expensive and labour intensive compared to snp arrays.
we show that snps often discarded in gwas, because they exhibit non-typical number of genotype clusters, can be informative of kir gene copy numbers.
by applying supervised classification, we are able to use qpcr results in a modest number of samples to impute copy numbers into a larger sample for which snp array signals are available.
we illustrate this method by applying it to two genes in the kir complex, kir3dl <dig> and kir3ds <dig>  which are suitable candidates for t1d association due to their interaction with hla class i molecules.
specifically, the kir3dl <dig> protein is known to interact with the hla class i allotypes that contain the hla-bw <dig> serological epitope  <cit> , whereas the protein encoded by kir3ds <dig>  which shares 97% sequence similarity to kir3dl <dig>  is thought to bind the more restrictive hla-bw4-80i epitope subset  <cit> .
the grouping of hla-a and hla-b alleles according to hla-bw <dig> serological epitope  <cit>  is given in additional file 1: table s <dig> and includes several hla class i alleles which are associated with t1d risk after conditioning on the major hla class ii effects  <cit> .
to date, kir3dl1/3ds <dig> association with t1d has only been studied using qpcr assays in limited sample sizes, which assess presence or absence of each kir gene  <cit> .
we used qpcr copy number calls in  <dig> samples as a training set, and imputed copy number in a further  <dig> samples from raw genotyping signals in snp array probes targeting the kir region.
we thus tested association of kir3dl1/3ds <dig> copy number with t1d, either directly, or through interaction with hla-bw <dig>  to the best of our knowledge, the sample size of our study is twenty-fold larger than any previous study of kir3dl1/3ds <dig> in t1d, and the first to test copy number variation rather than simply presence or absence  <cit> .
the hybrid method we advocate, leverages the information available from targeted qpcr assays in modest samples to the level of sample coverage required for modern, well-powered genetic studies.
it has the potential to be applied to other genes in the kir region or, indeed, to other chromosome regions that exhibit similar copy number variation and sequence complexity.
methods
subjects
dna was available from  <dig> individuals:  <dig> cases  from the genetic resource investigating diabetes  cohort, and  <dig> controls from the british  <dig> birth cohort .
all subjects were of white european ancestry with written informed consent and ethics committee/institutional review board approval.
the grid cohort was approved by the cambridgeshire  <dig> research ethics committee, study title “developing targets for diabetes prevention by the study of the genetics of type  <dig> diabetes” .
the 1958bc cohort was approved by the north west ethics committee, study title “ <dig> birth cohort tissue bank” .
the use of these samples was approved by the cambridgeshire  <dig> research ethics committee, study title “investigating genes and phenotypes associated with type  <dig> diabetes” .
ancestry was confirmed by pca analysis of earlier gwas data in these samples  <cit> .
the dna for the cases and controls was prepared using the same protocols in cambridge and in bristol respectively.
all samples were cell-line derived.
hla genotypes were available on a subset of  <dig> individuals,  <dig> cases and  <dig> controls.
hla-a and hla-b genes were typed at four-digit allele resolution using dynal reli sso assays  .
the epitope classification of hla-a and hla-b alleles is given in additional file 1: table s <dig>  all  <dig> samples have been genotyped using immunochip, a custom illumina 200k infinium high-density snp array  <cit> , according to the manufacturer’s protocol, processed at the university of virginia in charlottesville, usa.
a random subset of  <dig> samples,  <dig> cases and  <dig> controls, for which hla genotype was available, were selected for qpcr.
these samples were arrayed on bespoke 96-well plates, randomised half cases, half controls.
the source plates were chosen as those containing samples with the most hla typing available at the time .
design of multiplexed qpcr kir3dl1/3ds <dig> copy number assay
the qpcr platform used was the lightcycler  <dig> real-time pcr instrument on which we ran eighteen 384-well plates and four repeated plates.
the  <dig> samples,  <dig> cases and  <dig> controls, selected for qpcr were arrayed evenly, half-cases, half-controls, across the plates.
on the four repeated plates,  <dig> samples were arrayed .
on each plate, every sample was replicated across four neighbouring wells, resulting in a maximum of  <dig> samples per plate.
all plates, except for one, contained four repeated calibrator samples of known kir3dl1/3ds <dig> copy number that included, two samples with kir3ds1-kir3dl <dig> copy number 1- <dig> and two samples with kir3ds1-kir3dl <dig> copy number 2- <dig>  these are represented in additional file 1: figure s <dig> as black points.
to detect contamination, each plate also included one water well.
four plates were analysed in duplicate in order to assess reproducibility .
the qpcr probes and forward/reverse primers were carefully designed, in collaboration with jiang et al.
<cit> , to target and amplify most known kir3dl <dig> and kir3ds <dig> alleles, as well as the reference gene stat <dig>  known to always be present in two copies.
the probe and primer sequences are summarised in additional file 1: table s <dig>
each qpcr well was multiplexed, so that the copy numbers of kir3dl <dig>  kir3ds <dig> and stat <dig> were simultaneously assayed as part of the same qpcr reaction.
to allow for this, the probes were conjugated with three distinct dyes: fam for kir3ds <dig>  cy <dig> for kir3dl <dig> and dfo for stat <dig>
each qpcr well reaction was prepared with  <dig> μl of dna at  <dig> ng μl- <dig> and  <dig> μl of quantifast multiplex pcr mastermix .
the qpcr conditions were 95°c for  <dig> min, followed by  <dig> cycles at 95°c for  <dig> s and 66°c for  <dig> s. data was collected at 66°c.
quality control and normalisation of the qpcr data
the experiment files exported from the lightcycler gave us three crossingpoint  values per well, one for each of the dye-dna conjugates.
the ct value is representative of the number of qpcr cycles required for the dye-dna conjugate to be sufficiently amplified for the fluorescence to cross the detection threshold.
hence a larger ct value usually implies a smaller underlying copy number.
for each well, by subtracting the ct of fam-kir3dl <dig> and cy5-kir3ds <dig> from the ct value of the dfo-stat <dig>  the reference dye-dna conjugate, we obtained the baseline relative Δct value for kir3ds <dig> and kir3dl <dig>  then for each sample, we took the median over the four replicate wells to obtain per sample, kir3ds <dig> and kir3dl1Δct values .
however, certain wells did not yield a stat <dig> ct value since the detection threshold was not crossed within the  <dig> pcr cycles.
we found  <dig> samples that did not yield a dfo-stat <dig> ct reading in all four well replicates and these were excluded in the first step of our quality control .
visual inspection of the kir3dl <dig> and kir3ds1Δct distributions by plate led us to drop plate  <dig>  as it appeared to contain numerous outliers.
this resulted in a further  <dig> samples being dropped as part of our qc.
following qc, we were left with  <dig> unique samples,  <dig> cases and  <dig> controls, over  <dig> plates.
individual distributions of kir3ds <dig> and kir3dl1Δct were not aligned between plates , this prevented pooling of all plates for copy number calling.
to align the Δct distributions across the  <dig> plates, we first applied the k-medoids algorithm within each plate to the Δct kir3ds <dig> and kir3dl <dig> separately to identify the location of the most distinguishable copy number groups, one and two copies.
we then normalised across plates by a linear transformation so that the median Δct of the two groups mapped to  <dig> and  <dig> across all  <dig> plates .
after normalisation, negative Δct values were assigned to zero to reflect their expected copy number state.
following qc and normalisation, samples which were repeated across different plates showed good reproducibility .
these were summarised by the median of their repeated value.
copy number calling and multiple imputation in the subset of samples with qpcr data
samples with three or more missing fam-kir3dl <dig> or cy5-kir3ds <dig> ct values out of the four well replicates, were assigned to zero copies of kir3dl <dig> or kir3ds <dig> respectively.
for the remainder of the samples, copy number calling was done jointly on kir3dl <dig> and kir3ds <dig> using unsupervised clustering with a finite mixture model.
we called copy number groups by fitting a mixture of eight bivariate gaussian distributions to exploit the notable correlation between the normalised kir3ds <dig> and kir3dl1Δct values .
we allowed for eight kir3ds1-kir3dl <dig> copy number groups: three common groups of two copy numbers  and five rarer groups of lower or higher copy numbers  .
the mixture was fitted using an em algorithm  <cit>  with initial parameters calculated from the clusters returned by k-means with centers set to the eight expected locations of the copy number groups.
after fitting the mixture model, each sample was assigned eight posterior probabilities of belonging to each of the eight copy number groups, allowing for uncertainty in copy number calling.
these posterior probabilities were used to simulate ten plausible kir3ds1-kir3dl <dig> copy number datasets.
these ten multiply imputed datasets allowed for statistical analysis to be conducted in parallel and inference to be combined across datasets using the methods described by little and rubin   <cit>  and implemented in the mitools and mice r packages  <cit> .
copy number imputation into the extended samples
we extended our sample size by using the subset of samples common between the qpcr and snp datasets,  <dig> cases and  <dig> controls, to train a k-nearest neighbour  classifier to predict kir3dl1/3ds <dig> copy number using the r and θ signals from immunochip snps.
illumina arrays, such as the immunochip, have two fluorescent probes which differ on one base and allow discrimination of biallelic snps.
the fluorescent intensities are x and y, from which are derived the sum, r = x + y, and the ratio, tan=xy.
each of  <dig> snps lying within the kir3dl1/3ds <dig> region, were assessed for association with either kir3dl <dig> or kir3ds <dig> copy number in individual linear regression of copy number against r and θ .
nineteen snps out of  <dig> were associated , nine of which would have failed immunochip qc , with rs <dig> the most strongly predictive .
we compared running knn with all predictive snps or on various subsets, and found rs <dig> alone, with k= <dig>  minimised the mean leave-one-out cross-validation  error rate over ten multiply imputed qpcr datasets .
in each multiply imputed qpcr dataset, all samples were assigned a single imputed copy number group.
we also explored the effect of varying the size of the training data set by setting kir gene copy numbers to missing for a randomly chosen subset of samples and imputing them in the remaining samples .
testing for association of kir3dl1/3ds <dig> copy number with t1d
we tested for association of t1d with the predicted copy numbers from the qpcr and snp datasets using logistic regression.
we allowed for uncertainty in the copy number call when estimating individual odds ratios by using the ten multiply imputed datasets generated from the qpcr posterior probabilities  <cit> , and averaging the estimates over those with the r mitools package  <cit> .
we allowed for statistical interaction with hla-bw <dig> by repeating the association test in the subsets of carriers of the target ligand hla-bw <dig> epitopes, hla-bw <dig> for kir3dl <dig> and the putative ligand hla-bw4-80i for kir3ds <dig>  we directly tested for interaction with a more powerful case-only χ <dig> test  <cit> .
RESULTS
bivariate clustering enables accurate copy number calling in qpcr data
before normalisation, Δct distributions varied across plates preventing simple visual copy number assignment .
after normalisation, samples repeated across different plates showed good reproducibility .
furthermore, bivariate clustering, on both the kir3ds <dig> and kir3dl1Δct, enabled  <dig> samples to be confidently assigned to a single copy number group.
over 99% of all qpcr samples were classified with a posterior probability of copy number group membership greater than 99%.
further, samples of known copy number, the black points in additional file 1: figure s <dig> corresponding to kir3ds1-kir3dl <dig> copy number 2- <dig> and 1- <dig>  were assigned to the correct copy number group.
we allowed for the limited uncertainty in copy number calling, which mostly stems from distinguishing 2- <dig> from 3- <dig> , by means of multiple imputation of ten datasets as described in the methods.
imputation into extended samples by integration of snp and qpcr data
snp signals, r and θ, showed strong association with individual copy numbers of kir3dl1/3sd <dig> for  <dig> of  <dig> snps in the kir3dl <dig> region .
the strongest example is shown for snp rs <dig> in figure  <dig>  in which clusters can be discerned that correspond closely with qpcr derived kir3dl1/3sd <dig> copy numbers.
this figure also illustrates a number of important points regarding using snp signals for imputation.
firstly, θ corresponds to the ratio of copies of kir3dl <dig> to kir3ds <dig>  while r corresponds to the total copy number.
secondly, some clusters overlap, particularly along the r axis, making them hard to identify without the qpcr data.
consequently, a large proportion of samples were misclassified when we attempted unsupervised clustering using bivariate finite mixture model approaches, first with platinumcnv  <cit> , then with our own, mixture of beta-gaussian distributions, approach.
finally, the clusters are in slightly different positions in cases and controls, reflecting the known sensitivity of genotyping chips to subtle differences in dna preparation and storage conditions since they were prepared and processed in two different centers  <cit> .
instead, we used the qpcr copy numbers as training data to perform supervised classification with knn on the snp signals, which does not explicitly rely on the identification of distinct clusters.
we first explored the validity of our imputation approach by means of loocv in the samples with qpcr data.
we examined using all  <dig> predictive snps, or various subsets, and found optimal knn imputation was achieved with the single most predictive snp, rs <dig> .
by varying the size of the training data, we suggest that only  <dig> samples are required to achieve loocv error rates <5% and  <dig> for error rates < <dig> % .
no evidence of association of kir3dl1/3ds <dig> copy number with t1d
finally, we tested for association of kir3dl1/3ds <dig> copy number with t1d status.
we found no significant evidence of association, in the qpcr data , nor in the extended snp data , either overall or with any single copy number group .
kir3ds1
-
kir3dl1
, and in the marginal
kir3dl1
and
kir3ds1
copy number groups
kir3ds1-kir3dl1
kir3dl1
kir3ds1
no evidence of a significant, joint or marginal, effect was detected in the qpcr dataset,  <dig> cases and  <dig> controls, nor in the snp dataset,  <dig> cases and  <dig> controls.
case-control counts shown are derived from the most likely copy number assignment across the ten multiply imputed qpcr and snp datasets.
statistical inference for association is derived from the multiply imputed datasets using the r mitools package  <cit> .
the last row of each table contains the pooled p-value for each association test using the r mice package  <cit> .
by expanding to these large samples, which would be infeasible to genotype directly with qpcr, we are able to exclude odds ratios outside of the range  for the common copy number groups with 95% certainty.
we also repeated the association tests in the subset of individuals, carriers of the hla-bw <dig> epitope, and again detected no significant association .
a disadvantage of subsetting by hla-bw <dig> is that we lose power by greatly reducing the sample size.
a more powerful test for interaction between unlinked genes is a case-only test  <cit> .
if there were an interaction between kir3dl1/3ds <dig> and hla-bw <dig> then this should be detectable as a difference in kir3dl1/3ds <dig> copy number frequencies across hla-bw <dig> strata in the cases.
however, we found no significant evidence for association in either the qpcr or snp data sets, before or after summarising the kir copy number by presence/absence to increase power by reducing the degrees of freedom .
kir3ds1
-
kir3dl1
, and in the marginal
kir3dl1
and
kir3ds1
copy number groups
kir3ds1-kir3dl1
kir3dl1
kir3ds1
association is tested in the subset of individuals carriers of an hla-bw <dig> epitope for the joint kir3ds1-kir3dl <dig> and marginal kir3dl <dig> copy number groups and, also tested in the subset of individuals carriers of the hla-bw4-80i epitope for the marginal kir3ds <dig> copy number group.
case-control counts shown are derived from the most likely copy number assignment across the ten multiply imputed qpcr and snp datasets.
statistical inference for association is derived from the multiply imputed datasets using the r mitools package  <cit> .
the last row of each table contains the pooled p-value for each association test using the r mice package  <cit> .
χ
2
test for interaction between
kir3ds1
-
kir3dl1
and hla-bw <dig>  across the ten multiply imputed qpcr and snp datasets
counts in each contingency table are derived from the most likely copy number assignment across the multiply imputed datasets.
to reduce the degrees of freedom and improve power, we summarise copy numbers higher or equal to one by presence  and zero by absence .
the pooled p-value of each χ <dig> test, across the multiply imputed datasets, is given in the last row of each contingency table.
we find no significant association with hla-bw <dig>  within cases, in either the joint  or the marginal kir3ds1-kir3dl <dig> distributions.
CONCLUSIONS
regions with great allelic and copy number variation are difficult to properly assess using gwas.
while genome-wide snp arrays are typically cost effective ways to assay common genetic variation, very polymorphic regions can make the design of snp probes that bind uniquely to their target region, difficult or impossible.
this has resulted in low snp coverage in the kir region for the common snp arrays.
the snps that do exist on arrays are often discarded during the qc phase of any gwas because they do not exhibit the expected three clusters.
on the other hand, assaying individual genes can prove expensive.
for example, the qpcr assays used here to target kir3dl <dig> and kir3ds <dig> cost £ <dig> per sample.
further, qpcr derived data, despite careful design and multiplexing, remain subject to noise .
we ameliorated this through qc and normalisation across plates, and then jointly clustering kir3dl <dig> and kir3ds <dig>  to exploit the correlation between the Δct values.
the advantage of joint clustering is demonstrated in qpcr plate  <dig>  where noisy cases  are difficult to assign as one or two copies based solely on their kir3dl1Δct, but are much more clearly distinguishable when we also consider their kir3ds1Δct value .
as receptors for hla class i molecules, kir genes are important candidates for t1d and other diseases that associate with hla variation.
however, researchers have been unable to fully assess their candidacy due to lack of coverage in gwas and the complexity and expense of kir gene-specific assays.
so far, kir association studies for t1d have involved small samples sizes and have probed for presence/absence of multiple kir genes whilst ignoring the respective copy numbers of these genes, with mixed results and no consistent pattern of association emerging.
however, copy number variation in kir could be important, as it is a mechanism which gives rise to a large diversity of haplotypes  <cit> .
our hybrid approach, as outlined in additional file 1: figure s <dig>  allowed us to perform the  largest study of kir3dl1/3ds <dig> copy number in t1d to date, and to test for association in eight of the most frequent copy number groups .
in  <dig> samples, we found no association of kir3dl1/3ds <dig> copy number with t1d, alone or conditional on presence of the hla-bw <dig> epitope.
our results suggest that, despite the association of certain hla-a and hla-b alleles with t1d and the established biological interaction between hla-bw <dig> and kir3dl <dig>  copy number variation in kir3dl1/3ds <dig> is unlikely to have a significant effect on the risk of developing t1d.
other kir genes that are in high linkage disequilibrium with kir3dl <dig> and kir3ds <dig> are also unlikely to be associated.
according to the allele frequency net database  <cit> , these include kir2ds <dig>  and kir2dl <dig> , for kir3dl <dig> and, kir2dl <dig> , kir2ds <dig>  and kir2ds <dig> , for kir3ds <dig> .
thus, copy number variation in kir3dl1/3ds <dig> or neighbouring genes is unlikely to be an important risk factor in t1d.
in order to better understand why rs <dig> is the best snp for predicting copy number variation in kir3dl1/3ds <dig>  we used blat  <cit>  to match the probe sequences of rs <dig> on immunochip against the allelic sequences of all kir genes available from the immuno polymorphism database  <cit> .
interestingly, we found that the snp probes do not target kir3dl1/3ds <dig> but instead bind uniquely to the fifth intron of kir2dl <dig>  a neighbouring framework gene.
examining the kir2dl <dig> alleles matched by the rs <dig> probes, we discovered that the snp probes are in fact picking up copy number variation of kir2dl4* <dig>  an allele of kir2dl <dig> that undergoes copy number variation along with kir3dl1/3ds <dig> <cit> .
this explains the small but persistent misclassification error rate of 2% since our imputation is based on linkage disequilibrium between rs <dig> and kir3dl1/3ds <dig> rather than on perfect discrimination between our target genes.
we have identified  <dig> samples, which we believe, are consistently misclassified due to imperfect linkage disequilibrium .
we have observed other snps with more than three clusters that may correlate with copy number of other kir genes and, given the availability of qpcr results, could be imputable in a similar manner.
abbreviations
cnv: copy number variation; hla: human leukocyte antigen; kir: killer immunoglobulin-like receptor; knn: k-nearest neighbour; loocv: leave-one-out cross-validation; qpcr: quantitative polymerase chain reaction; snp: single nucleotide polymorphism; t1d: type  <dig> diabetes.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
conceived and designed the experiments: ja todd.
designed the qpcr assays: jj, wj, ja traherne and jt.
performed the qpcr experiments for kir3dl1/3ds1: djs.
performed immunochip genotyping: so-g, w-mc, pc and ssr.
prepared data: hs, nmw, hg and osb.
statistical analysis: np, cw, jmmh and hs.
wrote the manuscript: np and cw.
reviewed and revised the manuscript: jmmh, ds, ja traherne, nmw and ja todd.
all authors have read and approved the manuscript.
supplementary material
additional file 1
supplementary figures and tables.
click here for file
acknowledgements
this work was funded by the jdrf, the wellcome trust  and the national institute for health research  cambridge biomedical research centre.
the research leading to these results has received funding from the european union’s 7th framework programme  under grant agreement no.
<dig> .
we thank david dunger, barry widmer, and the british society for paediatric endocrinology and diabetes for the tid case collection.
we would like to thank the uk medical research council and wellcome trust for funding the collection of dna for the british  <dig> birth cohort .
dna control samples were prepared and provided by s. ring, r. jones, m. pembrey, w. mcardle, d. strachan and p. burton.
np is supported by the medical research council doctoral training grant .
cw and hg are funded by the wellcome trust .
jt and jat are supported by the medical research council  and the wellcome trust .
the cambridge institute for medical research  is in receipt of a wellcome trust strategic award .

 killer immunoglobulin-like receptors are transmembrane glycoproteins expressed by natural killer cells and subsets of t cells. the kir region lies in a kb gene cluster located within the mb leukocyte receptor complex on chr19q. the region exhibits great haplotype and copy number diversity, which has prevented complete assessment of the kir genes in standard genome-wide association studies, despite their strong candidacy for immune-related traits. targeted quantitative polymerase chain reaction assays have been used to detect presence or absence of individual kir genes and more recently, determine copy numbers. nevertheless these remain expensive and labour intensive compared to snp arrays. we show that snps often discarded in gwas, because they exhibit non-typical number of genotype clusters, can be informative of kir gene copy numbers. by applying supervised classification, we are able to use qpcr results in a modest number of samples to impute copy numbers into a larger sample for which snp array signals are available. we illustrate this method by applying it to two genes in the kir complex, kir3dl and kir3ds which are suitable candidates for t1d association due to their interaction with hla class i molecules. specifically, the kir3dl protein is known to interact with the hla class i allotypes that contain the hla-bw serological epitope, whereas the protein encoded by kir3ds which shares 97% sequence similarity to kir3dl is thought to bind the more restrictive hla-bw4-80i epitope subset. the grouping of hla-a and hla-b alleles according to hla-bw serological epitope is given in additional file 1: table s and includes several hla class i alleles which are associated with t1d risk after conditioning on the major hla class ii effects. to date, kir3dl1/3ds association with t1d has only been studied using qpcr assays in limited sample sizes, which assess presence or absence of each kir gene. we used qpcr copy number calls in samples as a training set, and imputed copy number in a further samples from raw genotyping signals in snp array probes targeting the kir region. we thus tested association of kir3dl1/3ds copy number with t1d, either directly, or through interaction with hla-bw to the best of our knowledge, the sample size of our study is twenty-fold larger than any previous study of kir3dl1/3ds in t1d, and the first to test copy number variation rather than simply presence or absence. the hybrid method we advocate, leverages the information available from targeted qpcr assays in modest samples to the level of sample coverage required for modern, well-powered genetic studies. it has the potential to be applied to other genes in the kir region or, indeed, to other chromosome regions that exhibit similar copy number variation and sequence complexity. methods subjects dna was available from individuals: cases from the genetic resource investigating diabetes cohort, and controls from the british birth cohort. all subjects were of white european ancestry with written informed consent and ethics committee/institutional review board approval. the grid cohort was approved by the cambridgeshire research ethics committee, study title “developing targets for diabetes prevention by the study of the genetics of type diabetes”. the 1958bc cohort was approved by the north west ethics committee, study title “ birth cohort tissue bank”. the use of these samples was approved by the cambridgeshire research ethics committee, study title “investigating genes and phenotypes associated with type diabetes”. ancestry was confirmed by pca analysis of earlier gwas data in these samples. the dna for the cases and controls was prepared using the same protocols in cambridge and in bristol respectively. all samples were cell-line derived. hla genotypes were available on a subset of individuals, cases and controls. hla-a and hla-b genes were typed at four-digit allele resolution using dynal reli sso assays. the epitope classification of hla-a and hla-b alleles is given in additional file 1: table s all samples have been genotyped using immunochip, a custom illumina 200k infinium high-density snp array, according to the manufacturer’s protocol, processed at the university of virginia in charlottesville, usa. a random subset of samples, cases and controls, for which hla genotype was available, were selected for qpcr. these samples were arrayed on bespoke 96-well plates, randomised half cases, half controls. the source plates were chosen as those containing samples with the most hla typing available at the time. design of multiplexed qpcr kir3dl1/3ds copy number assay the qpcr platform used was the lightcycler real-time pcr instrument on which we ran eighteen 384-well plates and four repeated plates. the samples, cases and controls, selected for qpcr were arrayed evenly, half-cases, half-controls, across the plates. on the four repeated plates, samples were arrayed. on each plate, every sample was replicated across four neighbouring wells, resulting in a maximum of samples per plate. all plates, except for one, contained four repeated calibrator samples of known kir3dl1/3ds copy number that included, two samples with kir3ds1-kir3dl copy number 1- and two samples with kir3ds1-kir3dl copy number 2- these are represented in additional file 1: figure s as black points. to detect contamination, each plate also included one water well. four plates were analysed in duplicate in order to assess reproducibility. the qpcr probes and forward/reverse primers were carefully designed, in collaboration with jiang et al., to target and amplify most known kir3dl and kir3ds alleles, as well as the reference gene stat known to always be present in two copies. the probe and primer sequences are summarised in additional file 1: table s each qpcr well was multiplexed, so that the copy numbers of kir3dl kir3ds and stat were simultaneously assayed as part of the same qpcr reaction. to allow for this, the probes were conjugated with three distinct dyes: fam for kir3ds cy for kir3dl and dfo for stat each qpcr well reaction was prepared with μl of dna at ng μl- and μl of quantifast multiplex pcr mastermix. the qpcr conditions were 95°c for min, followed by cycles at 95°c for s and 66°c for s. data was collected at 66°c. quality control and normalisation of the qpcr data the experiment files exported from the lightcycler gave us three crossingpoint values per well, one for each of the dye-dna conjugates. the ct value is representative of the number of qpcr cycles required for the dye-dna conjugate to be sufficiently amplified for the fluorescence to cross the detection threshold. hence a larger ct value usually implies a smaller underlying copy number. for each well, by subtracting the ct of fam-kir3dl and cy5-kir3ds from the ct value of the dfo-stat the reference dye-dna conjugate, we obtained the baseline relative Δct value for kir3ds and kir3dl then for each sample, we took the median over the four replicate wells to obtain per sample, kir3ds and kir3dl1Δct values. however, certain wells did not yield a stat ct value since the detection threshold was not crossed within the pcr cycles. we found samples that did not yield a dfo-stat ct reading in all four well replicates and these were excluded in the first step of our quality control. visual inspection of the kir3dl and kir3ds1Δct distributions by plate led us to drop plate as it appeared to contain numerous outliers. this resulted in a further samples being dropped as part of our qc. following qc, we were left with unique samples, cases and controls, over plates. individual distributions of kir3ds and kir3dl1Δct were not aligned between plates, this prevented pooling of all plates for copy number calling. to align the Δct distributions across the plates, we first applied the k-medoids algorithm within each plate to the Δct kir3ds and kir3dl separately to identify the location of the most distinguishable copy number groups, one and two copies. we then normalised across plates by a linear transformation so that the median Δct of the two groups mapped to and across all plates. after normalisation, negative Δct values were assigned to zero to reflect their expected copy number state. following qc and normalisation, samples which were repeated across different plates showed good reproducibility. these were summarised by the median of their repeated value. copy number calling and multiple imputation in the subset of samples with qpcr data samples with three or more missing fam-kir3dl or cy5-kir3ds ct values out of the four well replicates, were assigned to zero copies of kir3dl or kir3ds respectively. for the remainder of the samples, copy number calling was done jointly on kir3dl and kir3ds using unsupervised clustering with a finite mixture model. we called copy number groups by fitting a mixture of eight bivariate gaussian distributions to exploit the notable correlation between the normalised kir3ds and kir3dl1Δct values. we allowed for eight kir3ds1-kir3dl copy number groups: three common groups of two copy numbers and five rarer groups of lower or higher copy numbers. the mixture was fitted using an em algorithm with initial parameters calculated from the clusters returned by k-means with centers set to the eight expected locations of the copy number groups. after fitting the mixture model, each sample was assigned eight posterior probabilities of belonging to each of the eight copy number groups, allowing for uncertainty in copy number calling. these posterior probabilities were used to simulate ten plausible kir3ds1-kir3dl copy number datasets. these ten multiply imputed datasets allowed for statistical analysis to be conducted in parallel and inference to be combined across datasets using the methods described by little and rubin and implemented in the mitools and mice r packages. copy number imputation into the extended samples we extended our sample size by using the subset of samples common between the qpcr and snp datasets, cases and controls, to train a k-nearest neighbour classifier to predict kir3dl1/3ds copy number using the r and θ signals from immunochip snps. illumina arrays, such as the immunochip, have two fluorescent probes which differ on one base and allow discrimination of biallelic snps. the fluorescent intensities are x and y, from which are derived the sum, r = x + y, and the ratio, tan=xy. each of snps lying within the kir3dl1/3ds region, were assessed for association with either kir3dl or kir3ds copy number in individual linear regression of copy number against r and θ. nineteen snps out of were associated, nine of which would have failed immunochip qc, with rs the most strongly predictive. we compared running knn with all predictive snps or on various subsets, and found rs alone, with k= minimised the mean leave-one-out cross-validation error rate over ten multiply imputed qpcr datasets. in each multiply imputed qpcr dataset, all samples were assigned a single imputed copy number group. we also explored the effect of varying the size of the training data set by setting kir gene copy numbers to missing for a randomly chosen subset of samples and imputing them in the remaining samples. testing for association of kir3dl1/3ds copy number with t1d we tested for association of t1d with the predicted copy numbers from the qpcr and snp datasets using logistic regression. we allowed for uncertainty in the copy number call when estimating individual odds ratios by using the ten multiply imputed datasets generated from the qpcr posterior probabilities, and averaging the estimates over those with the r mitools package. we allowed for statistical interaction with hla-bw by repeating the association test in the subsets of carriers of the target ligand hla-bw epitopes, hla-bw for kir3dl and the putative ligand hla-bw4-80i for kir3ds we directly tested for interaction with a more powerful case-only χ test.  bivariate clustering enables accurate copy number calling in qpcr data before normalisation, Δct distributions varied across plates preventing simple visual copy number assignment. after normalisation, samples repeated across different plates showed good reproducibility. furthermore, bivariate clustering, on both the kir3ds and kir3dl1Δct, enabled samples to be confidently assigned to a single copy number group. over 99% of all qpcr samples were classified with a posterior probability of copy number group membership greater than 99%. further, samples of known copy number, the black points in additional file 1: figure s corresponding to kir3ds1-kir3dl copy number 2- and 1- were assigned to the correct copy number group. we allowed for the limited uncertainty in copy number calling, which mostly stems from distinguishing 2- from 3-, by means of multiple imputation of ten datasets as described in the methods. imputation into extended samples by integration of snp and qpcr data snp signals, r and θ, showed strong association with individual copy numbers of kir3dl1/3sd for of snps in the kir3dl region. the strongest example is shown for snp rs in figure in which clusters can be discerned that correspond closely with qpcr derived kir3dl1/3sd copy numbers. this figure also illustrates a number of important points regarding using snp signals for imputation. firstly, θ corresponds to the ratio of copies of kir3dl to kir3ds while r corresponds to the total copy number. secondly, some clusters overlap, particularly along the r axis, making them hard to identify without the qpcr data. consequently, a large proportion of samples were misclassified when we attempted unsupervised clustering using bivariate finite mixture model approaches, first with platinumcnv, then with our own, mixture of beta-gaussian distributions, approach. finally, the clusters are in slightly different positions in cases and controls, reflecting the known sensitivity of genotyping chips to subtle differences in dna preparation and storage conditions since they were prepared and processed in two different centers. instead, we used the qpcr copy numbers as training data to perform supervised classification with knn on the snp signals, which does not explicitly rely on the identification of distinct clusters. we first explored the validity of our imputation approach by means of loocv in the samples with qpcr data. we examined using all predictive snps, or various subsets, and found optimal knn imputation was achieved with the single most predictive snp, rs. by varying the size of the training data, we suggest that only samples are required to achieve loocv error rates <5% and for error rates < %. no evidence of association of kir3dl1/3ds copy number with t1d finally, we tested for association of kir3dl1/3ds copy number with t1d status. we found no significant evidence of association, in the qpcr data, nor in the extended snp data, either overall or with any single copy number group. kir3ds1 - kir3dl1, and in the marginal kir3dl1 and kir3ds1 copy number groups kir3ds1-kir3dl1 kir3dl1 kir3ds1 no evidence of a significant, joint or marginal, effect was detected in the qpcr dataset, cases and controls, nor in the snp dataset, cases and controls. case-control counts shown are derived from the most likely copy number assignment across the ten multiply imputed qpcr and snp datasets. statistical inference for association is derived from the multiply imputed datasets using the r mitools package. the last row of each table contains the pooled p-value for each association test using the r mice package. by expanding to these large samples, which would be infeasible to genotype directly with qpcr, we are able to exclude odds ratios outside of the range for the common copy number groups with 95% certainty. we also repeated the association tests in the subset of individuals, carriers of the hla-bw epitope, and again detected no significant association. a disadvantage of subsetting by hla-bw is that we lose power by greatly reducing the sample size. a more powerful test for interaction between unlinked genes is a case-only test. if there were an interaction between kir3dl1/3ds and hla-bw then this should be detectable as a difference in kir3dl1/3ds copy number frequencies across hla-bw strata in the cases. however, we found no significant evidence for association in either the qpcr or snp data sets, before or after summarising the kir copy number by presence/absence to increase power by reducing the degrees of freedom. kir3ds1 - kir3dl1, and in the marginal kir3dl1 and kir3ds1 copy number groups kir3ds1-kir3dl1 kir3dl1 kir3ds1 association is tested in the subset of individuals carriers of an hla-bw epitope for the joint kir3ds1-kir3dl and marginal kir3dl copy number groups and, also tested in the subset of individuals carriers of the hla-bw4-80i epitope for the marginal kir3ds copy number group. case-control counts shown are derived from the most likely copy number assignment across the ten multiply imputed qpcr and snp datasets. statistical inference for association is derived from the multiply imputed datasets using the r mitools package. the last row of each table contains the pooled p-value for each association test using the r mice package. χ 2 test for interaction between kir3ds1 - kir3dl1 and hla-bw across the ten multiply imputed qpcr and snp datasets counts in each contingency table are derived from the most likely copy number assignment across the multiply imputed datasets. to reduce the degrees of freedom and improve power, we summarise copy numbers higher or equal to one by presence and zero by absence. the pooled p-value of each χ test, across the multiply imputed datasets, is given in the last row of each contingency table. we find no significant association with hla-bw within cases, in either the joint or the marginal kir3ds1-kir3dl distributions.  regions with great allelic and copy number variation are difficult to properly assess using gwas. while genome-wide snp arrays are typically cost effective ways to assay common genetic variation, very polymorphic regions can make the design of snp probes that bind uniquely to their target region, difficult or impossible. this has resulted in low snp coverage in the kir region for the common snp arrays. the snps that do exist on arrays are often discarded during the qc phase of any gwas because they do not exhibit the expected three clusters. on the other hand, assaying individual genes can prove expensive. for example, the qpcr assays used here to target kir3dl and kir3ds cost £ per sample. further, qpcr derived data, despite careful design and multiplexing, remain subject to noise. we ameliorated this through qc and normalisation across plates, and then jointly clustering kir3dl and kir3ds to exploit the correlation between the Δct values. the advantage of joint clustering is demonstrated in qpcr plate where noisy cases are difficult to assign as one or two copies based solely on their kir3dl1Δct, but are much more clearly distinguishable when we also consider their kir3ds1Δct value. as receptors for hla class i molecules, kir genes are important candidates for t1d and other diseases that associate with hla variation. however, researchers have been unable to fully assess their candidacy due to lack of coverage in gwas and the complexity and expense of kir gene-specific assays. so far, kir association studies for t1d have involved small samples sizes and have probed for presence/absence of multiple kir genes whilst ignoring the respective copy numbers of these genes, with mixed results and no consistent pattern of association emerging. however, copy number variation in kir could be important, as it is a mechanism which gives rise to a large diversity of haplotypes. our hybrid approach, as outlined in additional file 1: figure s allowed us to perform the largest study of kir3dl1/3ds copy number in t1d to date, and to test for association in eight of the most frequent copy number groups. in samples, we found no association of kir3dl1/3ds copy number with t1d, alone or conditional on presence of the hla-bw epitope. our results suggest that, despite the association of certain hla-a and hla-b alleles with t1d and the established biological interaction between hla-bw and kir3dl copy number variation in kir3dl1/3ds is unlikely to have a significant effect on the risk of developing t1d. other kir genes that are in high linkage disequilibrium with kir3dl and kir3ds are also unlikely to be associated. according to the allele frequency net database, these include kir2ds and kir2dl, for kir3dl and, kir2dl, kir2ds and kir2ds, for kir3ds. thus, copy number variation in kir3dl1/3ds or neighbouring genes is unlikely to be an important risk factor in t1d. in order to better understand why rs is the best snp for predicting copy number variation in kir3dl1/3ds we used blat to match the probe sequences of rs on immunochip against the allelic sequences of all kir genes available from the immuno polymorphism database. interestingly, we found that the snp probes do not target kir3dl1/3ds but instead bind uniquely to the fifth intron of kir2dl a neighbouring framework gene. examining the kir2dl alleles matched by the rs probes, we discovered that the snp probes are in fact picking up copy number variation of kir2dl4* an allele of kir2dl that undergoes copy number variation along with kir3dl1/3ds. this explains the small but persistent misclassification error rate of 2% since our imputation is based on linkage disequilibrium between rs and kir3dl1/3ds rather than on perfect discrimination between our target genes. we have identified samples, which we believe, are consistently misclassified due to imperfect linkage disequilibrium. we have observed other snps with more than three clusters that may correlate with copy number of other kir genes and, given the availability of qpcr results, could be imputable in a similar manner. abbreviations cnv: copy number variation; hla: human leukocyte antigen; kir: killer immunoglobulin-like receptor; knn: k-nearest neighbour; loocv: leave-one-out cross-validation; qpcr: quantitative polymerase chain reaction; snp: single nucleotide polymorphism; t1d: type diabetes. competing interests the authors declare that they have no competing interests. authors’ contributions conceived and designed the experiments: ja todd. designed the qpcr assays: jj, wj, ja traherne and jt. performed the qpcr experiments for kir3dl1/3ds1: djs. performed immunochip genotyping: so-g, w-mc, pc and ssr. prepared data: hs, nmw, hg and osb. statistical analysis: np, cw, jmmh and hs. wrote the manuscript: np and cw. reviewed and revised the manuscript: jmmh, ds, ja traherne, nmw and ja todd. all authors have read and approved the manuscript. supplementary material additional file 1 supplementary figures and tables. click here for file acknowledgements this work was funded by the jdrf, the wellcome trust and the national institute for health research cambridge biomedical research centre. the research leading to these results has received funding from the european union’s 7th framework programme under grant agreement no.. we thank david dunger, barry widmer, and the british society for paediatric endocrinology and diabetes for the tid case collection. we would like to thank the uk medical research council and wellcome trust for funding the collection of dna for the british birth cohort. dna control samples were prepared and provided by s. ring, r. jones, m. pembrey, w. mcardle, d. strachan and p. burton. np is supported by the medical research council doctoral training grant. cw and hg are funded by the wellcome trust. jt and jat are supported by the medical research council and the wellcome trust. the cambridge institute for medical research is in receipt of a wellcome trust strategic award.
our approach could potentially be applied to other kir genes to allow cost effective assaying of gene copy number in large samples. however, their cost is prohibitive at the sample sizes required for detecting effects typically observed in complex genetic diseases. we find no association between kir3dl1/3ds copy number and t1d in cases and controls; a sample size twenty-fold larger than in any previous kir association study. due to our sample size, we can exclude odds ratios larger than for the common kir3dl1/3ds copy number groups at the 5% significance level. killer immunoglobulin-like receptors are surface receptors of natural killer cells that bind to their corresponding human leukocyte antigen class i ligands, making them interesting candidate genes for hla-associated autoimmune diseases, including type diabetes. however, allelic and copy number variation in the kir region effectively mask it from standard genome-wide association studies: single nucleotide polymorphism probes targeting the region are often discarded by standard genotype callers since they exhibit variable cluster numbers. we applied this method to assess copy number variation in two kir genes, kir3dl and kir3ds which are suitable candidates for t1d susceptibility since they encode the only kir molecules known to bind with hla-bw epitopes. five other kir genes, kir2ds kir2dl kir2dl kir2ds and kir2ds in high linkage disequilibrium with kir3dl and kir3ds are also unlikely to be significantly associated. first, we showed that noise and batch effects in multiplexed qpcr assays are addressed through normalisation and simultaneous copy number calling of multiple genes. then, we used supervised classification to impute copy numbers of specific kir genes from snp signals.
 killer immunoglobulin-like receptors are surface receptors of natural killer cells that bind to their corresponding human leukocyte antigen class i ligands, making them interesting candidate genes for hla-associated autoimmune diseases, including type diabetes. however, allelic and copy number variation in the kir region effectively mask it from standard genome-wide association studies: single nucleotide polymorphism probes targeting the region are often discarded by standard genotype callers since they exhibit variable cluster numbers. quantitative polymerase chain reaction assays address this issue. however, their cost is prohibitive at the sample sizes required for detecting effects typically observed in complex genetic diseases.  we propose a more powerful and cost-effective alternative, which combines signals from snps with more than three clusters found in existing datasets, with qpcr on a subset of samples. first, we showed that noise and batch effects in multiplexed qpcr assays are addressed through normalisation and simultaneous copy number calling of multiple genes. then, we used supervised classification to impute copy numbers of specific kir genes from snp signals. we applied this method to assess copy number variation in two kir genes, kir3dl and kir3ds which are suitable candidates for t1d susceptibility since they encode the only kir molecules known to bind with hla-bw epitopes. we find no association between kir3dl1/3ds copy number and t1d in cases and controls; a sample size twenty-fold larger than in any previous kir association study. due to our sample size, we can exclude odds ratios larger than for the common kir3dl1/3ds copy number groups at the 5% significance level.  we found no evidence of association of kir3dl1/3ds copy number with t1d, either overall or dependent on hla-bw epitope. five other kir genes, kir2ds kir2dl kir2dl kir2ds and kir2ds in high linkage disequilibrium with kir3dl and kir3ds are also unlikely to be significantly associated. our approach could potentially be applied to other kir genes to allow cost effective assaying of gene copy number in large samples. kir3dl1kir3ds1kir2ds4kir2dl3kir2dl5kir2ds5kir2ds1hla-bw4cnvqpcrimmunochipkirimputationt1d 