BACKGROUND
more than  <dig>  genome-wide association studies  of complex human diseases have been performed since 2008
<cit> .
the majority of gwas were conducted using common single nucleotide polymorphism  arrays targeting markers that were identified from the international hapmap project
<cit> .
these studies are based on the assumption that common traits are driven by common low-penetrance polymorphisms with a frequency of more than one or five percent in the population
<cit> .
a vast proportion of the heritability of complex traits remains unexplained
<cit> .
however, advances in genomic technologies now allow for the search of rare variants of modest to intermediate penetrance
<cit> .
snp arrays offer the possibility of rapid genotyping of thousands of samples with highly reliable results at low cost.
several commercial arrays now include a large fraction of rare single nucleotide variants  discovered by high-throughput sequencing technologies.
the latter, while still expensive compared to snp arrays, allows for the discovery of all variants, rare and common, located in the genome of sequenced individuals.
the illumina humanexome beadchip provides a compromise between genotyping snp arrays and next generation sequencing by enabling the genotyping of rare snvs in thousands of samples at relatively low cost.
the humanexome beadchip is enriched for rare and low frequency coding variations previously identified from the sequenced exomes of approximately  <dig>  individuals of diverse populations for variations seen in more than two individuals and in more than two sequencing projects
<cit> .
compared to other genotyping platforms targeting millions of markers, the proportion of rare variants  included in the humanexome beadchip is considerably larger.
a recent review of clustering tools for widely used illumina beadchip arrays was presented by ritchie et al.
<cit> .
they reported that some tools performed marginally better than others for common and rare variants .
the authors noted that methods borrowing information from other snps  to genotype rare variants could theoretically outperform reference-based methods, as these would suffer from the limited information available in the training sets for the homozygous and heterozygous clusters of rare alleles.
such within-sample methods are implemented in genosnp <cit> .
some tools, such as m <dig> <cit>  and opticall <cit> , use a mixture of between- and within-sample approaches.
other tools, such as gencall
<cit> , rely on a reference cluster file to cluster marker genotypes one at a time.
the zcall tool exclusively genotypes markers that have been previously “missed” by a another tool, and was also recently described
<cit> .
for this project, the performance of gencall, genosnp, opticall and zcall for clustering markers from the humanexome beadchip have been analysed and compared.
with the growing interest of the community for studies with rare variants
<cit> , this “head to head” comparison will provide guidance for study design, tool selection and results interpretation.
methods
clustering tools
four clustering tools were compared: gencall, part of the genomestudio software version v <dig>
<cit> , genosnp version  <dig>
<cit> , opticall version  <dig> .3
<cit>  and zcall version  <dig>
<cit> .
all four tools differ with respect to their genotype calling method.
gencall
gencall is illumina’s proprietary tool and is available through the genomestudio software.
for a complete description of this tool, refer to ritchie et al.
<cit> .
in brief, this tool uses the normalized microarray intensities for both alleles  to compute the associated polar coordinates  for each marker/sample pair.
next, using a between-sample model, meaning that it calls one marker by looking at the population of samples, it assigns genotypes by determining the nearest cluster using a reference containing the expected position of each genotype cluster for every marker as determined from the hapmap data.
if required, the user may modify the position of each of the expected clusters to be more representative of the data at hand.
by pre-assigning the expected position of each cluster, this method can readily provide a genotyping assignation of rare variants for studies having only a small number of samples.
however, due to experimental variabilities and genomic variations in different populations, the position of the observed cluster’s centroid might shift when compared with the expected one.
a considerable amount of manual cluster adjustments might be needed to achieve good genotype calls.
for this project, a custom cluster file was created by modifying the expected position of the genotype cluster’s centroid for a subset of markers by using all samples from the dataset.
the markers selected for manual inspection were:  markers with a high heterozygous frequency,  markers with a low mean intensity,  markers with a low call frequency,  markers with a low minor allele frequency with no heterozygous calls,  markers showing an excess of heterozygous calls,  markers with low aa t means or low bb t meansa,  mitochondrial markers,  markers on sex chromosomes,  markers that fail reproducibility tests,  markers with a small cluster separation or  markers with low quality score.
to compare the efficiency of this modified cluster file, the results from gencall with the original cluster file were also included in the analysis.
genosnp
genosnp uses a within-sample model, meaning that it assigns genotypes to all markers of a single sample at once
<cit> .
it uses raw x and y allele intensities extracted using the genomestudio software and calls genotypes by fitting a four-components mixture of student’s t-distributions on different subsets of markers  by the mean of a variational bayes expectation maximization algorithm .
the use of this method improves robustness by allowing uncertainty in the statistical model, in contrast to standard expectation maximisation methods.
genosnp computes the posterior probability of the marker genotype calls.
as this tool calls one sample at a time, it offers the flexibility of providing final sample genotyping results before the whole study dataset is ready to be processed.
it can also be parallelized by running the tool on multiple samples at a time, and it does not require a reference panel.
it is generally expected that this tool would perform well with rare variants, as their genotypes will be clustered with higher frequency variants according to measured x and y allele intensities .
to speed up the genotyping process, samples were called using genosnp as soon as they were released from the genotyping center.
a posterior probability cutoff of  <dig>  was used to achieve higher quality calls.
to ascertain the quality of the results once all samples were genotyped, the mean and the median intensities of all calls for each sample were plotted.
opticall
the opticall tool uses a mixture of between- and within-sample models.
it uses a subset of  intensities from random samples at a random marker to find a prior distribution to the statistical model used to call genotypes across markers.
this distribution is inferred by using an em algorithm to fit a four-class mixture of student’s t-distributions.
the initial values used by the em algorithm are obtained by using the kmeans++ algorithm
<cit>  and the individual genotype’s a priori probabilities are assumed to be uniform .
then, a second mixture of t-distributions is used for the between-sample clustering, where the previously estimated priors are used in a maximum-a-posteriori  estimate of its parameters.
to measure the quality of genotype calls, opticall relies on deviation from the hardy-weinberg equilibrium .
the tool will try to improve the genotype calls when the hwe test fails  by fitting the previously described model without a prior.
zcall
this tool functions as a post-processing tool
<cit> .
the zcall tool separates the clusters for rare variants by partitioning the  intensity space using horizontal and vertical thresholds.
their positions are derived from the mean and variance of the homozygote clusters for common variants that were previously called and are scaled according to a z-score factor to optimize concordance with the default tool.
genotypes are then assigned with respect to their position relative to the z-score scaled coordinates.
accordingly, rare variants are called by defining a distance threshold.
the homozygote threshold for the major allele is estimated from the first calling tool’s genotypes, and the rare allele’s threshold is estimated by linear regression from the means and standard deviations of x and y intensities of common markers.
as recommended by the authors, zcall was used as a post-processing step after gencall .
version  <dig>  was used, where all z thresholds were derived from genomestudio’s final report, from which samples were filtered out based on call rate and global heterozygosity.
a z threshold of  <dig> was used after comparing the concordance with the original calls .
only missing genotypes from the original genomestudio report were recalled by zcall.
dataset
the four tools were applied to a dataset consisting of  <dig>  unique samples from the montreal heart institute  cohort.
we also included  <dig> experimental replicates of na <dig>   <dig> replicates of na <dig> and  <dig> replicates of na <dig> from the coriell institute .
finally,  <dig> and  <dig> mhi cohort samples were replicated  <dig> and  <dig> times, respectively.
all  <dig>  samples  were genotyped using the illumina humanexome beadchip, assessing  <dig>  markers including  <dig> insertions/deletions .
some of these markers  are present in nhlbi go exome sequencing project  database
<cit> , and 93% of these are rare variants with a minor allele frequency  below 5% in the european american population according to the esp database
<cit> .
the research protocol was approved by the montreal heart institute research ethics review board and all participants signed an informed consent.
agreement between tools
we used cohen’s kappa  and fleiss’ pi , two widely used statistics, to compute the extent of agreement between raters
<cit> , or in this case, genotype calling tools.
cohen’s κ computes the extent of agreement between two tools by first computing the overall agreement probability , using a two-way contingency table , for the distribution of n samples by tools  and genotype category, where n
kl
indicates the number of samples that tool  <dig> and  <dig> classified into genotypes k and l, respectively.
pa=∑k=1qnkkn
in a comparable contingency table, cohen’s κ is estimated using equation  <dig>
γκ=pa-pe|κ1-pe|κ,wherepe|κ=∑k=1qnaknbkn <dig>
where p
a
is the observed proportion of agreement  and p
e|κ
is the proportion of agreement expected by chance.
the data can be summarized in a frequency table , where, for a given sample i and genotype k, r
ik
represents the total number of tools that called genotype k for sample i. fleiss’ π is then defined by equation 3
<cit> .
γ^π=pa-pe|π1-pe|πpa=1n∑i=1n∑k=1qrikrpe|π=∑k=1qπ^k <dig> andπ^k=1n∑i=1nrikr
the possible set of genotypes included the no call genotype, as all tools might agree that a marker is impossible to be categorized in either of the three genotype clusters  due to quality issues .
error rates
several methods for error rate estimation  with pedigree data have been proposed and reviewed by liu et al.
for their use with unrelated samples
<cit> .
the genotypic model  provides a proper estimation of the error rate and proposes different constraints on the parameter space in order to make the model parameters identifiable.
ϵ=2±2+1226wherec1≃naanc2≃nabnc3≃nbbn
the genotypic model was tested by liu et al.
for common variants.
however, we found that the possible values of ϵ were out of bounds  with the humanexome beadchip data where a majority of markers are rare.
this can be explained by the proportion of the minor allele in the population, p <dig>  which is almost null.
for these cases, ϵ was approximated using ϵ≃/ <dig> .
RESULTS
clustering quality
genosnp
the genosnp tool returns the probability of belonging to one of the three genotype clusters  for each evaluated genotype.
the maximal probability is used to define the genotype to be called.
we observed a majority of samples with a high proportion of low quality calls .
cluster plots were created for some of the poorly performing samples which raised concern about the calling quality .
suspecting a lack of convergence, modifications to the tool were made to increase the number of iterations of the vb-em procedure.
this improved the quality of the clustering for some of the samples .
further modifications to the initial x and y intensities’ variance parameter from  <dig>  to  <dig>  greatly improved the quality of the clustering tool .
to efficiently ascertain the effect of the optimized parameters on the calling quality, plots of the mean and the median of the maximal probability of each sample have been created , showing a net increase of the probabilities.
to further improve the comparison between the original genosnp and the modified version , the two tools  were used for comparison with the other tools in this study.
zcall
according to goldstein et al.
<cit> , the optimal value of the z threshold should be determined by trying different values of z to find the one with the most concordance to the original gencall calls.
here, the optimal value of z was determined to be  <dig>  having a concordance of  <dig> % with the original data.
missing rates
one important property of a calling tool is its capacity to assign a genotype to the majority of samples and markers .
the sample and marker missing rates of the six tools were compared .
the original version of genosnp had the highest missing rate .
the optimized version of genosnp increased both the mean sample and marker calling rates by  <dig> % .
these rates were inferior to those of the other tools:  <dig> %,  <dig> %,  <dig> % and  <dig> % for gencall , opticall and zcall, respectively.
it is important to note that the missing rate of zcall is bound to be less than or equal to that of gencall, as zcall will only call missing genotypes from the results produced by gencall.
also, the missing rates of gencall were slightly better when using the optimized cluster file when compared to the original one.
a non-parametric friedman rank sum test comparing tools showed a statistically significant difference in both sample and marker missing rates , even though the rates seemed similar .
pairwise dependent-samples non-parametric sign tests comparing pairwise tools were also significant .
more specifically,  <dig> - and  <dig> -fold decreases, on average, were observed for sample and marker missing rates, respectively, when comparing the original and the optimized cluster file for gencall.
precision estimates
the dataset contained a high number of technical replicates.
the concordance of genotype calls between replicates was computed for the five tools .
the optimized version of genosnp increased concordance between sample duplicates from a mean of  <dig> % to  <dig> %, but did not exceed the performance of gencall .
even though zcall reduced the missing rate from a mean of  <dig> % to  <dig> %, it added variability in the newly called genotypes, slightly decreasing the concordance between replicated samples when compared with the original results .
overall, the call concordance between replicates was similar .
gencall had the highest concordance rate with a mean of  <dig> % between replicated sample pairs.
optimizing gencall’s cluster file had only a minor effect on the mean concordance, with a difference of  <dig> ×10-3%.
accuracy estimates
one sample sequenced by the  <dig> genomes project was included three times to the original dataset to assess the concordance with next-generation sequencing.
the genotypes called from the different tools were compared to the third release of the  <dig> genomes project, comprised of low coverage whole genome and high coverage exome sequencing data.
the comparison with this gold standard was performed using the pygenclean tool
<cit> .
since the majority of markers in the humanexome beadchip are located in exons, the results were mostly compared to the high coverage sequencing data .
table
<dig> shows the concordance for the three replicates of na <dig>  apart from the original genosnp, all the tools had concordance rates greater than 99%.
gencall  had the highest concordance rate , closely followed by zcall .
gencall  also performed well, with a mean concordance of  <dig> %.
then, calls were categorized according to their allele content: either  homozygous calls for the common allele or  involving the rare allele .
the frequencies were computed using the corresponding dataset .
table
<dig> shows the concordance for the three replicates of na <dig> for genotypes called as homozygous for the common allele.
apart from the original genosnp on one replicate, all had concordance rates greater than 99%.
gencall  had the highest concordance rate , closely followed by gencall  and zcall .
table
<dig> summarizes the concordance rates with the gold standard for genotypes involving the rare allele.
gencall  had the highest concordance rate , closely followed by zcall .
the other tool had a lower concordance rate .
call concordance and number of compared markers for the three control replicates when compared to the  <dig> genomes project.
the following six tools were compared: gencall , genosnp , opticall  and zcall.
call concordance and number of compared markers for the three control replicates when compared to the  <dig> genomes project.
for each dataset , only genotypes called as homozygous of the common allele  were kept for analysis.
the following six tools were compared: gencall , genosnp , opticall  and zcall.
call concordance and number of compared markers for the three control replicates when compared to the  <dig> genomes project.
for each dataset , only genotypes called as heterozygous or homozygous of the rare allele  were kept for analysis.
the following six tools were compared: gencall , genosnp , opticall  and zcall.
inter-tool agreement
to estimate tool agreement, three coefficients were computed: cohen’s κ and percent agreement , and fleiss’ π .
since gencall  provided the best result set, its agreement with the other three tools  was considered .
the best agreement was between gencall and zcall, as expected due to the dependence of zcall on the results of gencall.
both metrics showed a high number of outliers according to the interquartile range dispersionb.
to assess the overall agreement between the four tools , fleiss’ κ was computed .
the agreement was fairly good between the tools when comparing all markers .
there was a total of  <dig>  outliers according to the interquartile range dispersion, a majority of which  were rare variants  and 5% were markers that were zeroed out while optimizing gencall’s cluster file.
the proportion of rare variants  was significantly higher in the outlier group .
out of the approximately three thousands outlier markers that were comparable with the  <dig> genomes project data, gencall had the highest concordance for the three replicates of na <dig> , followed by zcall  .
error rate estimates
error rates were estimated by using the genotypic model , where the six tools were evaluated .
except for the original version of genosnp, all tools showed a comparable estimated error rate mean .
the zcall tool provided the lowest estimated error rate with a mean of  <dig> %, followed by gencall  with a mean of  <dig> %.
a non-parametric friedman rank sum test comparing tools showed a significant difference in the error rate distributions .
pairwise dependent-samples non-parametric sign tests comparing pairwise tools were also significant .
<cit> .
CONCLUSIONS
this study compares the performance of widely used clustering tools when applied to genotyping data from illumina’s humanexome beadchip.
this genotyping array includes a large proportion of rare variants that were previously identified by sequencing technologies.
the dataset used here included  <dig>  unique samples along with a high number of technical replicates for quality assessment.
contrary to our original expectations, genosnp, which relies on a within-sample model, did not perform well when used straight “out of the box”.
this might be explained by the high density of markers in each beadpool, which is higher than for illumina’s human 1m beadchip .
the latter was successfully tested by genosnp’s authors.
this problem was mostly obvious with cluster plots  and with graphs showing the summarized quality of the calls per sample .
the concordance of results from the original genosnp tool with the  <dig> genomes project, however, remained high .
the optimized tool  increased the concordance to a mean of  <dig> %.
the quality threshold of  <dig>  provided a better separation of the three clusters, but increased the missing rate of both sample and marker .
gencall relies on a between-sample model that requires reference parameters to perform its clustering.
as such, it is common practice to manually modify the reference cluster parameters to ensure the best quality of results when the population analyzed is different from the one used to generate the reference parameters.
this task requires a significant amount of manual labor which increases with the number of samples and markers.
loading the raw data, normalizing and modifying the original cluster file took one person  <dig> work days, compared to only a few days to generate the intensity files for the other tools.
when we compared the genotype results from gencall using the original cluster file with the ones generated with the optimized cluster file, we saw only limited improvement in the overall concordance of genotype calls with those of the gold standard.
when partitioning calls according to their allele content, we saw a limited decrease in the concordance of homozygous calls involving the common allele when compared with the gold standard, while the concordance of calls involving the rare allele improved by approximately  <dig> %.
a limited improvement in the concordance between technical replicates was observed.
however, the optimization of the cluster file had a greater impact on the missing rate per sample and per loci .
according to normal quality control procedures, a filter to remove samples and markers with a missing rate greater than 2% is typically imposed prior to genetic analysis
<cit> .
by optimizing the cluster file, an average of  <dig>  markers  per sample could be rescued by lowering the marker missing rate below the 2% quality threshold.
it should be noted that this study was limited by the lack of an independent gold standard for concordance analysis.
gencall’s  concordance with the  <dig> genomes project might be overestimated due to the nature of its reference parameters.
indeed, the reference cluster file was created by using estimated cluster position for each of the markers in the hapmap population.
hence, the sample na <dig> is expected to have a higher concordance value than samples from the montreal heart institute cohort owing to the ceu hapmap population.
one possible way to overcome this limitation would be to sequence a few of our own cohort samples, however, it could be argued that the sequencing technology itself would not provide an adequate gold standard in this situation.
many study designs will plan to include experimental replicates chosen from the genotyped cohort to assess the reproducibility of the genotyping pipeline and the precision of the results.
including a reference sample to the study design offers the additional possibility of assessing the accuracy of the results.
high precision is particularly important as it provides optimal power for statistical analysis
<cit>  and can prevent type i errors due to plate biases or subgroup effects
<cit> .
gencall  provided the highest concordance rate between experimental replicates  when compared to other tools.
the other tools  performed relatively well, all providing a mean concordance greater than 99%.
all clustering tools had a high accuracy  when calling common markers .
other metrics have shown that all the tools  performed almost at the same efficiency level on the humanexome beadchip.
the major difference arose when in the presence of rare markers, where the accuracy of all the tools decreased below  <dig> %  for genotypes involving the rare allele.
within-sample tools like genosnp can process samples in an efficient manner by running single samples and smaller batches of samples in parallel instead of having to wait for a large amount of samples to be genotyped and normalized.
unfortunately, even with the proper optimization of initial parameters, genosnp could not outperform the other tools.
however, gencall, illumina’s proprietary tool, performed better than the other tested tools with respect to concordance with the gold standard for genotypes involving the rare allele  and slightly better for the concordance in between technical replicates .
the fact that zcall has been run as a post-process of gencall without the use of the reference cluster file optimization  might explain why it’s accuracy was not as high as the optimized gencall for calls involving the rare allele.
since the third version of zcall derives its thresholds from genomestudio’s report, the call rate will increase and better accuracy might be possible if the original cluster file is optimized beforehand.
recommending a single clustering tool according to the metrics shown in this report is not straightforward.
in general, gencall  outperformed the other tools in terms of precision and accuracy .
its accuracy was also higher for the markers with low inter-tool agreement.
however, when using the optimized cluster file, gencall’s accuracy for the homozygous calls  was lower than when using the default cluster file.
when considering missing and estimated error rates, zcall outperformed the other tools, closely followed by gencall .
it is important to mention that the task of optimizing the cluster file is time consuming.
furthermore, all the other tools presented here require intensity data provided by the genomestudio software and possible file conversion, which increase the total execution time.
the parallel use of multiple clustering tools offers the possibility of identifying discordant markers which can be further investigated.
but notably, the manual optimization of gencall’s cluster file at those loci and the visual inspection of the cluster plots should provide high quality datasets for downstream analysis.
endnotes
aθ values of the center of the aa and bb clusters in normalized polar coordinates, respectively.
b outliers are observations that fall below q1- <dig>  or above q3+ <dig> , where q <dig> and q <dig> are respectively the first and third quartiles and iqr=q3-q <dig> is the interquartile range.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
lplp worked on the experimental design, optimized genosnp’s initial parameters, performed analyses, compared the tools, implemented the statistical tests and drafted the manuscript.
mal contributed to the analysis and to the manuscript.
ab contributed to the statistical analyses and the manuscript.
sp contributed to the experimental design, analysis and the manuscript.
vn contributed to the manual optimization of gencall and the manuscript.
jct leads the montreal heart institute cohort and directs the genomic project with the humanexome beadchip.
mpd supervised the project, participated in the design and coordination of the study, supervised the laboratory and revised the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
additional materials.
supplemental equation s1: the genotypic model for error rate estimation was tested by liu et al.
for common variants only.
however, we found that the possible values of ϵ were out of bound  for a majority of rare markers.
for those cases, ϵ was approximated using ϵ≃/ <dig>  supplemental table s1: distribution of n samples by calling tool in q categories.
the set of possible categories are all possible genotypes .
this table is computed for each marker and for each pair of calling tools.
the overall agreement probability and cohen’s κ are shown in equation  <dig> and  <dig> of the main text, respectively.
supplemental table s2: distribution of r calling tools by n samples and q response categories.
the set of possible categories are all possible genotypes .
this table is computed for each marker and for each calling tool.
fleiss’ π is explained in equation  <dig> of the main text.
supplemental table s3: call concordance and number of compared markers for the three control replicates when compared to the  <dig> genomes project for the markers that were outliers for their fleiss’ π values.
the following four tools were compared: gencall , genosnp , opticall  and zcall.
click here for file
acknowledgements
this project was funded by the montreal heart institute foundation and by the centre of excellence in personalized medicine.
we would like to acknowledge the team at the beaulieu-saucier université de montréal pharmacogenomics center who performed the genotyping of more than  <dig>  samples of the montreal heart institute cohort on the humanexome beadchip together with all the participants of the cohort and the ones responsible for their recruitment.

 more than genome-wide association studies of complex human diseases have been performed since 2008. the majority of gwas were conducted using common single nucleotide polymorphism arrays targeting markers that were identified from the international hapmap project. these studies are based on the assumption that common traits are driven by common low-penetrance polymorphisms with a frequency of more than one or five percent in the population. a vast proportion of the heritability of complex traits remains unexplained. however, advances in genomic technologies now allow for the search of rare variants of modest to intermediate penetrance. snp arrays offer the possibility of rapid genotyping of thousands of samples with highly reliable results at low cost. several commercial arrays now include a large fraction of rare single nucleotide variants discovered by high-throughput sequencing technologies. the latter, while still expensive compared to snp arrays, allows for the discovery of all variants, rare and common, located in the genome of sequenced individuals. the illumina humanexome beadchip provides a compromise between genotyping snp arrays and next generation sequencing by enabling the genotyping of rare snvs in thousands of samples at relatively low cost. the humanexome beadchip is enriched for rare and low frequency coding variations previously identified from the sequenced exomes of approximately individuals of diverse populations for variations seen in more than two individuals and in more than two sequencing projects. compared to other genotyping platforms targeting millions of markers, the proportion of rare variants included in the humanexome beadchip is considerably larger. a recent review of clustering tools for widely used illumina beadchip arrays was presented by ritchie et al.. they reported that some tools performed marginally better than others for common and rare variants. the authors noted that methods borrowing information from other snps to genotype rare variants could theoretically outperform reference-based methods, as these would suffer from the limited information available in the training sets for the homozygous and heterozygous clusters of rare alleles. such within-sample methods are implemented in genosnp. some tools, such as m and opticall, use a mixture of between- and within-sample approaches. other tools, such as gencall, rely on a reference cluster file to cluster marker genotypes one at a time. the zcall tool exclusively genotypes markers that have been previously “missed” by a another tool, and was also recently described. for this project, the performance of gencall, genosnp, opticall and zcall for clustering markers from the humanexome beadchip have been analysed and compared. with the growing interest of the community for studies with rare variants, this “head to head” comparison will provide guidance for study design, tool selection and results interpretation. methods clustering tools four clustering tools were compared: gencall, part of the genomestudio software version v, genosnp version, opticall version.3 and zcall version. all four tools differ with respect to their genotype calling method. gencall gencall is illumina’s proprietary tool and is available through the genomestudio software. for a complete description of this tool, refer to ritchie et al.. in brief, this tool uses the normalized microarray intensities for both alleles to compute the associated polar coordinates for each marker/sample pair. next, using a between-sample model, meaning that it calls one marker by looking at the population of samples, it assigns genotypes by determining the nearest cluster using a reference containing the expected position of each genotype cluster for every marker as determined from the hapmap data. if required, the user may modify the position of each of the expected clusters to be more representative of the data at hand. by pre-assigning the expected position of each cluster, this method can readily provide a genotyping assignation of rare variants for studies having only a small number of samples. however, due to experimental variabilities and genomic variations in different populations, the position of the observed cluster’s centroid might shift when compared with the expected one. a considerable amount of manual cluster adjustments might be needed to achieve good genotype calls. for this project, a custom cluster file was created by modifying the expected position of the genotype cluster’s centroid for a subset of markers by using all samples from the dataset. the markers selected for manual inspection were: markers with a high heterozygous frequency, markers with a low mean intensity, markers with a low call frequency, markers with a low minor allele frequency with no heterozygous calls, markers showing an excess of heterozygous calls, markers with low aa t means or low bb t meansa, mitochondrial markers, markers on sex chromosomes, markers that fail reproducibility tests, markers with a small cluster separation or markers with low quality score. to compare the efficiency of this modified cluster file, the results from gencall with the original cluster file were also included in the analysis. genosnp genosnp uses a within-sample model, meaning that it assigns genotypes to all markers of a single sample at once. it uses raw x and y allele intensities extracted using the genomestudio software and calls genotypes by fitting a four-components mixture of student’s t-distributions on different subsets of markers by the mean of a variational bayes expectation maximization algorithm. the use of this method improves robustness by allowing uncertainty in the statistical model, in contrast to standard expectation maximisation methods. genosnp computes the posterior probability of the marker genotype calls. as this tool calls one sample at a time, it offers the flexibility of providing final sample genotyping results before the whole study dataset is ready to be processed. it can also be parallelized by running the tool on multiple samples at a time, and it does not require a reference panel. it is generally expected that this tool would perform well with rare variants, as their genotypes will be clustered with higher frequency variants according to measured x and y allele intensities. to speed up the genotyping process, samples were called using genosnp as soon as they were released from the genotyping center. a posterior probability cutoff of was used to achieve higher quality calls. to ascertain the quality of the results once all samples were genotyped, the mean and the median intensities of all calls for each sample were plotted. opticall the opticall tool uses a mixture of between- and within-sample models. it uses a subset of intensities from random samples at a random marker to find a prior distribution to the statistical model used to call genotypes across markers. this distribution is inferred by using an em algorithm to fit a four-class mixture of student’s t-distributions. the initial values used by the em algorithm are obtained by using the kmeans++ algorithm and the individual genotype’s a priori probabilities are assumed to be uniform. then, a second mixture of t-distributions is used for the between-sample clustering, where the previously estimated priors are used in a maximum-a-posteriori estimate of its parameters. to measure the quality of genotype calls, opticall relies on deviation from the hardy-weinberg equilibrium. the tool will try to improve the genotype calls when the hwe test fails by fitting the previously described model without a prior. zcall this tool functions as a post-processing tool. the zcall tool separates the clusters for rare variants by partitioning the intensity space using horizontal and vertical thresholds. their positions are derived from the mean and variance of the homozygote clusters for common variants that were previously called and are scaled according to a z-score factor to optimize concordance with the default tool. genotypes are then assigned with respect to their position relative to the z-score scaled coordinates. accordingly, rare variants are called by defining a distance threshold. the homozygote threshold for the major allele is estimated from the first calling tool’s genotypes, and the rare allele’s threshold is estimated by linear regression from the means and standard deviations of x and y intensities of common markers. as recommended by the authors, zcall was used as a post-processing step after gencall. version was used, where all z thresholds were derived from genomestudio’s final report, from which samples were filtered out based on call rate and global heterozygosity. a z threshold of was used after comparing the concordance with the original calls. only missing genotypes from the original genomestudio report were recalled by zcall. dataset the four tools were applied to a dataset consisting of unique samples from the montreal heart institute cohort. we also included experimental replicates of na replicates of na and replicates of na from the coriell institute. finally, and mhi cohort samples were replicated and times, respectively. all samples were genotyped using the illumina humanexome beadchip, assessing markers including insertions/deletions. some of these markers are present in nhlbi go exome sequencing project database, and 93% of these are rare variants with a minor allele frequency below 5% in the european american population according to the esp database. the research protocol was approved by the montreal heart institute research ethics review board and all participants signed an informed consent. agreement between tools we used cohen’s kappa and fleiss’ pi, two widely used statistics, to compute the extent of agreement between raters, or in this case, genotype calling tools. cohen’s κ computes the extent of agreement between two tools by first computing the overall agreement probability, using a two-way contingency table, for the distribution of n samples by tools and genotype category, where n kl indicates the number of samples that tool and classified into genotypes k and l, respectively. pa=∑k=1qnkkn in a comparable contingency table, cohen’s κ is estimated using equation γκ=pa-pe|κ1-pe|κ,wherepe|κ=∑k=1qnaknbkn where p a is the observed proportion of agreement and p e|κ is the proportion of agreement expected by chance. the data can be summarized in a frequency table, where, for a given sample i and genotype k, r ik represents the total number of tools that called genotype k for sample i. fleiss’ π is then defined by equation 3. γ^π=pa-pe|π1-pe|πpa=1n∑i=1n∑k=1qrikrpe|π=∑k=1qπ^k andπ^k=1n∑i=1nrikr the possible set of genotypes included the no call genotype, as all tools might agree that a marker is impossible to be categorized in either of the three genotype clusters due to quality issues. error rates several methods for error rate estimation with pedigree data have been proposed and reviewed by liu et al. for their use with unrelated samples. the genotypic model provides a proper estimation of the error rate and proposes different constraints on the parameter space in order to make the model parameters identifiable. ϵ=2±2+1226wherec1≃naanc2≃nabnc3≃nbbn the genotypic model was tested by liu et al. for common variants. however, we found that the possible values of ϵ were out of bounds with the humanexome beadchip data where a majority of markers are rare. this can be explained by the proportion of the minor allele in the population, p which is almost null. for these cases, ϵ was approximated using ϵ≃/.  clustering quality genosnp the genosnp tool returns the probability of belonging to one of the three genotype clusters for each evaluated genotype. the maximal probability is used to define the genotype to be called. we observed a majority of samples with a high proportion of low quality calls. cluster plots were created for some of the poorly performing samples which raised concern about the calling quality. suspecting a lack of convergence, modifications to the tool were made to increase the number of iterations of the vb-em procedure. this improved the quality of the clustering for some of the samples. further modifications to the initial x and y intensities’ variance parameter from to greatly improved the quality of the clustering tool. to efficiently ascertain the effect of the optimized parameters on the calling quality, plots of the mean and the median of the maximal probability of each sample have been created, showing a net increase of the probabilities. to further improve the comparison between the original genosnp and the modified version, the two tools were used for comparison with the other tools in this study. zcall according to goldstein et al., the optimal value of the z threshold should be determined by trying different values of z to find the one with the most concordance to the original gencall calls. here, the optimal value of z was determined to be having a concordance of % with the original data. missing rates one important property of a calling tool is its capacity to assign a genotype to the majority of samples and markers. the sample and marker missing rates of the six tools were compared. the original version of genosnp had the highest missing rate. the optimized version of genosnp increased both the mean sample and marker calling rates by %. these rates were inferior to those of the other tools: %, %, % and % for gencall, opticall and zcall, respectively. it is important to note that the missing rate of zcall is bound to be less than or equal to that of gencall, as zcall will only call missing genotypes from the results produced by gencall. also, the missing rates of gencall were slightly better when using the optimized cluster file when compared to the original one. a non-parametric friedman rank sum test comparing tools showed a statistically significant difference in both sample and marker missing rates, even though the rates seemed similar. pairwise dependent-samples non-parametric sign tests comparing pairwise tools were also significant. more specifically, - and -fold decreases, on average, were observed for sample and marker missing rates, respectively, when comparing the original and the optimized cluster file for gencall. precision estimates the dataset contained a high number of technical replicates. the concordance of genotype calls between replicates was computed for the five tools. the optimized version of genosnp increased concordance between sample duplicates from a mean of % to %, but did not exceed the performance of gencall. even though zcall reduced the missing rate from a mean of % to %, it added variability in the newly called genotypes, slightly decreasing the concordance between replicated samples when compared with the original results. overall, the call concordance between replicates was similar. gencall had the highest concordance rate with a mean of % between replicated sample pairs. optimizing gencall’s cluster file had only a minor effect on the mean concordance, with a difference of ×10-3%. accuracy estimates one sample sequenced by the genomes project was included three times to the original dataset to assess the concordance with next-generation sequencing. the genotypes called from the different tools were compared to the third release of the genomes project, comprised of low coverage whole genome and high coverage exome sequencing data. the comparison with this gold standard was performed using the pygenclean tool. since the majority of markers in the humanexome beadchip are located in exons, the results were mostly compared to the high coverage sequencing data. table shows the concordance for the three replicates of na apart from the original genosnp, all the tools had concordance rates greater than 99%. gencall had the highest concordance rate, closely followed by zcall. gencall also performed well, with a mean concordance of %. then, calls were categorized according to their allele content: either homozygous calls for the common allele or involving the rare allele. the frequencies were computed using the corresponding dataset. table shows the concordance for the three replicates of na for genotypes called as homozygous for the common allele. apart from the original genosnp on one replicate, all had concordance rates greater than 99%. gencall had the highest concordance rate, closely followed by gencall and zcall. table summarizes the concordance rates with the gold standard for genotypes involving the rare allele. gencall had the highest concordance rate, closely followed by zcall. the other tool had a lower concordance rate. call concordance and number of compared markers for the three control replicates when compared to the genomes project. the following six tools were compared: gencall, genosnp, opticall and zcall. call concordance and number of compared markers for the three control replicates when compared to the genomes project. for each dataset, only genotypes called as homozygous of the common allele were kept for analysis. the following six tools were compared: gencall, genosnp, opticall and zcall. call concordance and number of compared markers for the three control replicates when compared to the genomes project. for each dataset, only genotypes called as heterozygous or homozygous of the rare allele were kept for analysis. the following six tools were compared: gencall, genosnp, opticall and zcall. inter-tool agreement to estimate tool agreement, three coefficients were computed: cohen’s κ and percent agreement, and fleiss’ π. since gencall provided the best result set, its agreement with the other three tools was considered. the best agreement was between gencall and zcall, as expected due to the dependence of zcall on the results of gencall. both metrics showed a high number of outliers according to the interquartile range dispersionb. to assess the overall agreement between the four tools, fleiss’ κ was computed. the agreement was fairly good between the tools when comparing all markers. there was a total of outliers according to the interquartile range dispersion, a majority of which were rare variants and 5% were markers that were zeroed out while optimizing gencall’s cluster file. the proportion of rare variants was significantly higher in the outlier group. out of the approximately three thousands outlier markers that were comparable with the genomes project data, gencall had the highest concordance for the three replicates of na, followed by zcall. error rate estimates error rates were estimated by using the genotypic model, where the six tools were evaluated. except for the original version of genosnp, all tools showed a comparable estimated error rate mean. the zcall tool provided the lowest estimated error rate with a mean of %, followed by gencall with a mean of %. a non-parametric friedman rank sum test comparing tools showed a significant difference in the error rate distributions. pairwise dependent-samples non-parametric sign tests comparing pairwise tools were also significant..  this study compares the performance of widely used clustering tools when applied to genotyping data from illumina’s humanexome beadchip. this genotyping array includes a large proportion of rare variants that were previously identified by sequencing technologies. the dataset used here included unique samples along with a high number of technical replicates for quality assessment. contrary to our original expectations, genosnp, which relies on a within-sample model, did not perform well when used straight “out of the box”. this might be explained by the high density of markers in each beadpool, which is higher than for illumina’s human 1m beadchip. the latter was successfully tested by genosnp’s authors. this problem was mostly obvious with cluster plots and with graphs showing the summarized quality of the calls per sample. the concordance of results from the original genosnp tool with the genomes project, however, remained high. the optimized tool increased the concordance to a mean of %. the quality threshold of provided a better separation of the three clusters, but increased the missing rate of both sample and marker. gencall relies on a between-sample model that requires reference parameters to perform its clustering. as such, it is common practice to manually modify the reference cluster parameters to ensure the best quality of results when the population analyzed is different from the one used to generate the reference parameters. this task requires a significant amount of manual labor which increases with the number of samples and markers. loading the raw data, normalizing and modifying the original cluster file took one person work days, compared to only a few days to generate the intensity files for the other tools. when we compared the genotype results from gencall using the original cluster file with the ones generated with the optimized cluster file, we saw only limited improvement in the overall concordance of genotype calls with those of the gold standard. when partitioning calls according to their allele content, we saw a limited decrease in the concordance of homozygous calls involving the common allele when compared with the gold standard, while the concordance of calls involving the rare allele improved by approximately %. a limited improvement in the concordance between technical replicates was observed. however, the optimization of the cluster file had a greater impact on the missing rate per sample and per loci. according to normal quality control procedures, a filter to remove samples and markers with a missing rate greater than 2% is typically imposed prior to genetic analysis. by optimizing the cluster file, an average of markers per sample could be rescued by lowering the marker missing rate below the 2% quality threshold. it should be noted that this study was limited by the lack of an independent gold standard for concordance analysis. gencall’s concordance with the genomes project might be overestimated due to the nature of its reference parameters. indeed, the reference cluster file was created by using estimated cluster position for each of the markers in the hapmap population. hence, the sample na is expected to have a higher concordance value than samples from the montreal heart institute cohort owing to the ceu hapmap population. one possible way to overcome this limitation would be to sequence a few of our own cohort samples, however, it could be argued that the sequencing technology itself would not provide an adequate gold standard in this situation. many study designs will plan to include experimental replicates chosen from the genotyped cohort to assess the reproducibility of the genotyping pipeline and the precision of the results. including a reference sample to the study design offers the additional possibility of assessing the accuracy of the results. high precision is particularly important as it provides optimal power for statistical analysis and can prevent type i errors due to plate biases or subgroup effects. gencall provided the highest concordance rate between experimental replicates when compared to other tools. the other tools performed relatively well, all providing a mean concordance greater than 99%. all clustering tools had a high accuracy when calling common markers. other metrics have shown that all the tools performed almost at the same efficiency level on the humanexome beadchip. the major difference arose when in the presence of rare markers, where the accuracy of all the tools decreased below % for genotypes involving the rare allele. within-sample tools like genosnp can process samples in an efficient manner by running single samples and smaller batches of samples in parallel instead of having to wait for a large amount of samples to be genotyped and normalized. unfortunately, even with the proper optimization of initial parameters, genosnp could not outperform the other tools. however, gencall, illumina’s proprietary tool, performed better than the other tested tools with respect to concordance with the gold standard for genotypes involving the rare allele and slightly better for the concordance in between technical replicates. the fact that zcall has been run as a post-process of gencall without the use of the reference cluster file optimization might explain why it’s accuracy was not as high as the optimized gencall for calls involving the rare allele. since the third version of zcall derives its thresholds from genomestudio’s report, the call rate will increase and better accuracy might be possible if the original cluster file is optimized beforehand. recommending a single clustering tool according to the metrics shown in this report is not straightforward. in general, gencall outperformed the other tools in terms of precision and accuracy. its accuracy was also higher for the markers with low inter-tool agreement. however, when using the optimized cluster file, gencall’s accuracy for the homozygous calls was lower than when using the default cluster file. when considering missing and estimated error rates, zcall outperformed the other tools, closely followed by gencall. it is important to mention that the task of optimizing the cluster file is time consuming. furthermore, all the other tools presented here require intensity data provided by the genomestudio software and possible file conversion, which increase the total execution time. the parallel use of multiple clustering tools offers the possibility of identifying discordant markers which can be further investigated. but notably, the manual optimization of gencall’s cluster file at those loci and the visual inspection of the cluster plots should provide high quality datasets for downstream analysis. endnotes aθ values of the center of the aa and bb clusters in normalized polar coordinates, respectively. b outliers are observations that fall below q1- or above q3+, where q and q are respectively the first and third quartiles and iqr=q3-q is the interquartile range. competing interests the authors declare that they have no competing interests. authors’ contributions lplp worked on the experimental design, optimized genosnp’s initial parameters, performed analyses, compared the tools, implemented the statistical tests and drafted the manuscript. mal contributed to the analysis and to the manuscript. ab contributed to the statistical analyses and the manuscript. sp contributed to the experimental design, analysis and the manuscript. vn contributed to the manual optimization of gencall and the manuscript. jct leads the montreal heart institute cohort and directs the genomic project with the humanexome beadchip. mpd supervised the project, participated in the design and coordination of the study, supervised the laboratory and revised the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 additional materials. supplemental equation s1: the genotypic model for error rate estimation was tested by liu et al. for common variants only. however, we found that the possible values of ϵ were out of bound for a majority of rare markers. for those cases, ϵ was approximated using ϵ≃/ supplemental table s1: distribution of n samples by calling tool in q categories. the set of possible categories are all possible genotypes. this table is computed for each marker and for each pair of calling tools. the overall agreement probability and cohen’s κ are shown in equation and of the main text, respectively. supplemental table s2: distribution of r calling tools by n samples and q response categories. the set of possible categories are all possible genotypes. this table is computed for each marker and for each calling tool. fleiss’ π is explained in equation of the main text. supplemental table s3: call concordance and number of compared markers for the three control replicates when compared to the genomes project for the markers that were outliers for their fleiss’ π values. the following four tools were compared: gencall, genosnp, opticall and zcall. click here for file acknowledgements this project was funded by the montreal heart institute foundation and by the centre of excellence in personalized medicine. we would like to acknowledge the team at the beaulieu-saucier université de montréal pharmacogenomics center who performed the genotyping of more than samples of the montreal heart institute cohort on the humanexome beadchip together with all the participants of the cohort and the ones responsible for their recruitment.
the reproducibility of genotype calls for each tool was assessed using experimental duplicates which provided concordance rates above 99%. while the four tools provided concordance with the gold standard above 99% for common alleles, some of them performed poorly for rare alleles. the inter-tool agreement of genotype calls was high for approximately 95% of variants. the use of experimental replicates provided a valuable quality control tool for genotyping projects with rare variants. we have compared the performance of four clustering tools for the genotyping of over samples using the illumina’s humanexome beadchip, which includes variants, 90% of which have a minor allele frequency below 5% in a population of european ancestry. the performance of clustering tools for common polymorphisms is well established, while their performance when conducted with a large proportion of rare variants is less known. along with the improvement of high throughput sequencing technologies, the genetics community is showing marked interest for the rare variants/common diseases hypothesis. while sequencing can still be prohibitive for large studies, commercially available genotyping arrays targeting rare variants prove to be a reasonable alternative. different reference parameters for gencall and different initial parameters for genosnp were tested. concordance of genosnp’s calls with the gold standard was below expectations and was increased by changing the tool’s initial parameters. the genosnp clustering tool could not be run straight “out of the box” with the humanexome beadchip, as modification of hard coded parameters was necessary to achieve optimal performance. overall, gencall marginally outperformed the other tools for the humanexome beadchip.
 along with the improvement of high throughput sequencing technologies, the genetics community is showing marked interest for the rare variants/common diseases hypothesis. while sequencing can still be prohibitive for large studies, commercially available genotyping arrays targeting rare variants prove to be a reasonable alternative. a technical challenge of array based methods is the task of deriving genotype classes by clustering intensity data points. the performance of clustering tools for common polymorphisms is well established, while their performance when conducted with a large proportion of rare variants is less known. we have compared the performance of four clustering tools for the genotyping of over samples using the illumina’s humanexome beadchip, which includes variants, 90% of which have a minor allele frequency below 5% in a population of european ancestry. different reference parameters for gencall and different initial parameters for genosnp were tested. genotyping accuracy was assessed using data from the genomes project as a gold standard, and agreement between tools was measured.  concordance of genosnp’s calls with the gold standard was below expectations and was increased by changing the tool’s initial parameters. while the four tools provided concordance with the gold standard above 99% for common alleles, some of them performed poorly for rare alleles. the reproducibility of genotype calls for each tool was assessed using experimental duplicates which provided concordance rates above 99%. the inter-tool agreement of genotype calls was high for approximately 95% of variants. most tools yielded similar error rates, except for zcall which performed better with a mean error rate.  the genosnp clustering tool could not be run straight “out of the box” with the humanexome beadchip, as modification of hard coded parameters was necessary to achieve optimal performance. overall, gencall marginally outperformed the other tools for the humanexome beadchip. the use of experimental replicates provided a valuable quality control tool for genotyping projects with rare variants. 