BACKGROUND
in population genetics, inference of population structures is important for various purposes such as assessment of genetic diversity, detection of genetic discontinuities in natural wildlife habitats, and correction for stratification in association studies.
to infer population structures without prior knowledge about the population, various statistical approaches using neutral molecular markers have been proposed  <cit> .
bayesian approaches using markov chain monte carlo  methods have been widely used to infer population structures since pritchard et al.
<cit>  proposed the bayesian clustering algorithms implemented in the well-known program structure.
this program can infer the assignment of individuals to populations or the admixture proportions of individuals for a given number of populations .
researchers have extended bayesian algorithms for various purposes such as to take advantage of spatial information  <cit> , estimate inbreeding coefficients  <cit> , allow for allele mutations  <cit> , and infer k values  <cit> .
pella and masuda  <cit>  used a dirichlet process  to infer k values.
dp is a stochastic process that was proposed by ferguson  <cit>  to treat nonparametric problems in bayesian frameworks.
the merit of using dp to infer k is that k can take any value between  <dig> and the number of individuals , and thus, few assumptions about k are needed for inference.
pella and masuda  <cit>  considered k and the assignment of individuals to populations as random variables using dp as a prior distribution for k and allele frequencies unique to populations.
huelsenbeck and andolfatto  <cit>  also used the dp prior for the inference of population structures, and reich and bondell  <cit>  proposed a clustering algorithm using the dp prior, which can incorporate spatial information.
besides the inference of population structures, dp priors have been used to infer the number of ancestral haplotype blocks  <cit> , to model nonsynonymous/synonymous rate ratios  <cit> , and to model the selfing rates of individuals  <cit> .
to date, two clustering programs that implement the dp prior have been provided, hwler  <cit>  and structurama  <cit> .
both programs implement the gibbs sampling procedure to infer the posterior distribution.
these programs differ in their approach to improve the mixing of mcmc algorithms.
hwler implements the sequentially-allocated merge-split  sampler, which moves multiple observations simultaneously  <cit> , and structurama implements the metropolis-coupled mcmc  technique  <cit> , which runs multiple chains, some of which are closer to a uniform distribution than the target distribution, and attempts to swap states among chains.
although hwler and structurama are useful and have been used in some recent studies  <cit> , their application to real data sets has been less common compared with that of structure.
this may be because the properties of these methods have not been investigated in detail.
when results obtained with different methods are inconsistent, it is difficult to interpret the results.
for example, hwler and structurama may provide inconsistent results.
in addition, although hwler detects three populations, the result obtained with structure under the assumption that k is  <dig> may be different from that obtained with hwler.
the purpose of this study was to characterize the bayesian method based on the dp prior and to provide information that will be useful in practice.
first, we evaluated the sams sampler because its effectiveness has not been examined.
second, we investigated the effect of a hyperparameter, named λ, which defines the prior distribution of allele frequencies, on the performance of this method.
as described by pella and masuda  <cit> , hwler set λ to  where jl is the number of alleles at locus l for any loci.
however, this specification resulted in inaccurate inferences in some cases.
third, we compared the dp prior method with other bayesian methods that infer the assignment of individuals for a given k value.
we focused primarily on the effect of unbalanced sample sizes among populations on the behavior of these methods because unbalanced sample sizes were found to affect the behavior of these methods in our preliminary study.
methods
assumption and goal
we assumed that a sampled data set consists of diploid individuals derived from an unknown number of populations, which are defined by unique values of allele frequencies.
although our novel program, named dpart, can also analyze haploid populations, we assumed only diploid individuals in this study for simplicity.
our goal was to infer the number of populations and assign individuals to their populations based on their genotypes.
the estimation of admixture proportions was not of interest.
dp prior
the mathematical details of dp can be found in studies conducted by ferguson  <cit> , antoniak  <cit> , and neal  <cit> .
here we provided a brief description of the dp prior.
approximately and intuitively, dp can be seen as a stochastic process that converts a continuous distribution to a discrete distribution.
we assumed that observations  are generated from the following model .
let g <dig> denote the prior distribution of allele frequencies at the locus, which is continuous.
dp is defined by g <dig> and the concentration parameter α .
dp divides g <dig> into a number of classes, each of which is represented by a single point and yields a discrete probability distribution, g. the number of classes, which is determined by the number of observations and α, can be infinite.
a vector of allele frequencies  for each genotype is drawn from g instead of g <dig>  and genotypes are drawn from the corresponding allele frequency vectors.
this model can be written as follows:
where xi is the genotype of individual i, ϕi is the allele frequency vector for genotype i, and dp  is dp.
this model is known as the dp mixture model  <cit> .
because g is discrete, allele frequency vectors for some genotypes may have values in common, i.e., these genotypes can be seen as members of the same population, which is characterized by the shared allele frequency vectors.
parameters
by integrating out g, we can obtain a simpler representation of the model.
when g is integrated out, the predictive  distribution for the allele frequency vector ϕi conditional on {ϕ <dig>  ..., ϕi-1} can be written as follows:
where nj is the number of allele frequency vectors that share values with ϕj.
this sequence of predictive distributions is known as a polya urn scheme  <cit> .
to represent the clustering property of the model more explicitly, we let {ϕ <dig> ...,ϕk} denote allele frequency vectors unique to populations { <dig> ...,k}.
in addition, we introduced parameters that represent the partition of individuals { <dig> ...,n}, according to the parameterization of dahl  <cit> .
we let η = {s <dig> ...,sk} define a partition for { <dig> ...,n} such that , si ∩ sj = ∅ for all i ≠ j, and si ≠ ∅ for all i.
eq.
includes the following two types of prior information:
where η<i is η before assigning the ith individual and |sj| is the number of individuals included in sj, which is the prior for the number of populations  and assignment of individuals, and
which is the prior for allele frequencies for each population.
now, we can rewrite eq.
as follows:
where x = {x <dig> ...,xn} is a vector of the genotypes of the individuals { <dig> ...,n}, i{i ∈ sj} =  <dig> if i ∈ sj, and otherwise  <dig>  and
where Γ is the gamma function.
note that eq.
results from the products of eq.
.
integration of allele frequencies
although some clustering methods infer allele frequencies , allele frequencies can be analytically integrated out as in hwler and structurama.
thus, our goal is to infer not f  but f .
by bayes' theorem, this can be written as follows:
where
where  denotes the genotypes of individuals included in sj and  is the posterior distribution of ϕj updated from g <dig> on the basis of the genotypes of individuals assigned to sj preceding the ith individual.
g <dig> is assumed to be a flat dirichlet distribution.
the prior probability density of ϕ is given as follows:
where l is the number of loci, jl is the number of alleles at locus l, λl is the hyperparameter for locus l, and ϕjlh is the frequency of allele h at locus l of population j. rannala and mountain  <cit>  provided the following equation:
where yjlh is the number of copies of allele h at locus l in individuals assigned to sj preceding the ith individual and yjl = ∑yjlh.
l in eq.
can be rewritten as follows:
where δ =  <dig>  if the genotype of individual i at locus l is heterozygous and if not is  <dig>  and nilh is the number of copies of allele h at locus l in individual i. now the integral in eq.
can be solved analytically  as follows:
gibbs sampler
neal  <cit>  and maceachern  <cit>  introduced the gibbs sampling procedure, which integrates out model parameters .
let η-i denote η when only the ith individual is removed.
the conditional probability of the assignment of individual i is given as follows:
where b is the normalizing constant.
because individuals are exchangeable, one can treat every individual as the last observation.
the integrals in eq.
can be solved analytically.
the details of gibbs sampling is described in endnote a. in our experiments, one mcmc iteration consisted of one scan of the gibbs sampler.
after burn-in, η is sampled in every predefined number of mcmc iterations.
sams sampler
because the gibbs sampler moves only one observation, a large change hardly ever occurs, and tends to stay at a local mode.
a remedy is to move multiple observations simultaneously.
jain and neal  <cit>  proposed a split-merge sampler, which proposes a new partition by splitting or merging components and accepts it with metropolis-hastings  probability.
when it splits a component, the proposed assignment of observations in the component is generated by repeating gibbs sampling only for these observations.
this method can produce more acceptable proposals than the random assignment of observations.
dahl  <cit>  improved the method in terms of efficiency and proposed an alternative sampler, the sams sampler, which generates proposals by one cycle of allocation of observations instead of repeating gibbs sampling.
the algorithm of the sams sampling is described in endnote b. in this study, the sams sampling was performed once in an mcmc iteration.
four iterations with the sams sampler were followed by one iteration with the gibbs sampler.
this cycle was repeated until the end of the iterations.
inference of the hyperparameter
λl, which is the hyperparameter for the prior distribution of allele frequencies at locus l, can be inferred by the mh method.
the joint posterior probability of η and λ = {λ <dig> ...,λl} can be written as follows:
when the prior distribution for λl is assumed to be uniform, equations in the gibbs and sams samplers can be used without modification.
in this study, the prior for λl was assumed to be uniform, u. the proposal of λl, λl', was generated from the normal distribution n, where δ is an arbitrary value, e.g.,  <dig> .
the proposal was accepted with the probability:
where λ-l is λ when λl is removed.
we can also assume that each locus shares a single value, λa.
in this case the probability becomes the following:
prior distribution of the number of populations
as seen in eq.
or , the prior distribution of k depends on α and the number of individuals.
we can infer this prior distribution by the monte carlo method .
for example, when n =  <dig> and α =  <dig> , the expected k value is approximately  <dig> .
when α increases to  <dig> , k increases to approximately  <dig> .
although it is possible to infer α as well as λ, in this study, we fixed α through the mcmc iterations.
the effect of α has been thoroughly examined by huelsenbeck and andolfatto  <cit> .
the authors reported that the misspecification of α  could affect the results, especially when the number of loci was small.
summarizing the sampled partitions
two methods were used to summarize the partitions that were sampled from the posterior distribution of η. for simulated data sets, we used the mean partition proposed by huelsenbeck and andolfatto  <cit> , which is defined as follows:
where v is the number of sampled partitions and d is the partition distance between the sampled partition ηi and mean partition.
the partition distance between ηi and u is defined as the minimum number of individuals that must be removed from both ηi and u such that the two partitions are the same  <cit> .
the partition distance was calculated using eq.
of konovalov et al.
<cit>  to obtain the cost matrix, and we solved this using the hungarian algorithm.
the algorithm for calculating the mean partition is described in endnote d. this statistic is useful for evaluating the posterior distribution of η automatically.
we divided the partition distances by the number of individuals.
thus, the partition distances ranged from  <dig> to  <dig>
for real data sets, we performed agglomerative hierarchical clustering on the basis of co-assignment probabilities.
this method was introduced by dawson and belkhir  <cit> .
briefly, after all samples were obtained, the co-assignment probabilities for all individual pairs were calculated from the sampled partitions.
next, complete linkage clustering was performed on the basis of the probabilities.
we used the hclust function of the stats package of r for complete linkage clustering  <cit> .
an example r code is provided in additional file  <dig>  note that neither the mean partition nor hierarchical clustering are affected by the label-switching problem that often emerges in analyses using methods implementing the dp prior.
programs and mcmc parameters
we wrote the clustering programs in c. the program implementing the dp prior is referred to as dpart.
in addition to dpart, we also produced two clustering programs that infer only the assignment of individuals under a given k value.
these programs are based on algorithms proposed by pritchard et al.
<cit>  and falush et al.
<cit> .
in these algorithms, the prior probability that an individual belongs to each population is equal among the populations, and individuals are assigned to populations according to likelihood.
these algorithms infer the allele frequencies of populations that are integrated out in dpart in addition to the assignment of individuals.
one program assumes that there is no correlation among the allele frequencies of populations, and the other assumes a correlation.
thus, the former is equivalent to the "no admixture and no f model" setting of structure and the latter is equivalent to the "no admixture and f model" setting.
these programs are referred to as fixed k and uncorrelated model  and fixed k and correlated model , respectively.
structure infers admixture proportions of individuals until burn-in is half completed even when no admixture option is selected .
because this setting helps the program to avoid the generation of an empty cluster, which has no individuals, we followed the setting.
in fcm, the prior distribution of the drift parameter was assumed to be u. the programs provided in this study are summarized in table  <dig>  dpart is provided in additional files  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>  fum and fcm can be reproduced by structure.
structurama was used for simulated data sets.
the expected k value was set equal to the true k value.
four chains were run simultaneously and the temperature was set to  <dig> .
we also used structure ver.
<dig>  with its default setting, "admixture and f model," for the bull data set because it is a widely used method for the inference of population structures, although the estimation of admixture proportions was not of interest.
for comparison between dpart and structurama and examination of the effect of λ, the number of mcmc iterations was  <dig> , and the first half of the iterations was discarded as burn-in.
sampling was performed every  <dig> iterations.
for comparison among dpart, fum, and fcm, the number of iterations was  <dig> , and the first half of the iterations was discarded.
sampling was performed every  <dig> iterations.
for the real data sets, the number of iterations was  <dig>  and the first  <dig>  iterations were discarded.
sampling was performed every  <dig> iterations.
simulation method 1
this simulation method was used to compare among clustering methods.
easypop was used to simulate populations under the wright-fisher model  <cit> .
the island model was assumed.
the number of populations was  <dig> and the number of individuals per population was  <dig> .
the number of generations, number of loci, mutation rate, and possible number of alleles were  <dig> ,  <dig>   <dig> , and  <dig>  respectively, for microsatellites, and  <dig> ,  <dig>   <dig> × 10- <dig>  and  <dig>  respectively, for single nucleotide polymorphisms .
new alleles occurred in every possible allelic state with equal probability.
diploidy, random mating, one sex, and linkage equilibrium were also assumed.
the genotypes of the first generation were generated by randomly assigning alleles from the possible allelic states.
migration rates were  <dig> ,  <dig> , and  <dig>  for microsatellites and  <dig>  for snps.
easypop was run five times for each scenario.
the resulting data sets are referred to as base population sets.
sampling was performed  <dig> times from each base population set to obtain  <dig> data sets per scenario.
when two, three or four populations were sampled from the base population sets, base populations to be sampled were selected randomly.
similarly, when  <dig> or  <dig> loci were sampled, the loci to be sampled were also selected randomly.
the pairwise fst values between the base populations for each migration rate are summarized in table  <dig>  typically, the observed heterozygosities were approximately  <dig>  for microsatellites and  <dig>  for snps.
the pairwise fst and heterozygosity values were calculated using genetix  <cit> .
m indicates the migration rate.
simulation method 2
this method was used to investigate the effect of λ. the variance of a dirichlet distribution, dirichlet , decreases as λ increases.
when allele frequencies are drawn from the dirichlet distribution, it means that the frequencies closely approach uniformity among alleles as λ increases.
thus, it was expected that λ would affect clustering behavior depending on the uniformity of frequencies among alleles.
to examine this hypothesis, we devised a simulation method on the basis of the allele frequency correlated model  <cit> .
in this method, allele frequencies were generated from the dirichlet distribution
where pa is the allele frequency of a virtual ancestral population, f is the drift parameter, which is equivalent to wright's fst, and ja is the number of ancestral alleles.
we can determine any level of uniformity of frequencies among alleles by varying ancestral allele frequencies.
we assumed two marker types, microsatellites  and snps .
for microsatellites, the ancestral allele frequencies were { <dig> , <dig> ,..., <dig> } or { <dig> , <dig> ,..., <dig> } and f was  <dig> , and for snps, the ancestral frequencies were { <dig> , <dig> } or { <dig> , <dig> } and f was  <dig> .
k was assumed to be  <dig> and the number of individuals was  <dig> per population.
genotypes were generated from the allele frequencies assuming random mating.
if only one allele was observed at a locus, then that locus was excluded.
one hundred data sets were generated for each scenario.
comparison among methods for simulated data sets
to compare methods for simulated data sets, we used the partition distance between the true and inferred partitions, which was denoted as .
we calculated average  over the  <dig> simulated data sets and counted the number of data sets in which  was  <dig>  or less.
for dpart and structurama, we used the mean partition as the inferred partition and calculated the average k values in the mean partitions.
when fum and fcm were used, the inferred partition was determined on the basis of the probabilities of assigning individuals to populations, which were calculated from the sampled partitions.
although fum and fcm may be evaluated with the mean partition, we used the partition that was based on the probabilities of assignment because this is computationally more feasible than the mean partition and label-switching occurs rarely in structure-like algorithms as indicated by pritchard et al.
<cit> .
both approaches provided almost the same partitions in our preliminary study.
note that  can not be used for evaluation without modification when unbalanced sample sizes are present among populations.
suppose that two populations are included in a data set.
if the sample sizes are uniform between the populations and an analysis fails to detect any population structures, i.e., all individuals are assigned to one cluster,  is  <dig> .
however, if the sample sizes are  <dig> and  <dig> and an analysis fails similarly,  decreases to  <dig>  .
thus, in such cases, we calculated  as if an individual in the smaller subset was equivalent to  <dig> individuals in the larger subset.
for example, if the sample sizes are  <dig> and  <dig> and only one individual in the smaller subset is incorrectly assigned to the larger subset,  is  <dig>   but not  <dig>  .
chicken data set
this data set represents  <dig> chickens of  <dig> european breeds, which were genotyped for  <dig> microsatellites by rosenberg et al.
<cit> .
this data set was previously analyzed by pella and masuda  <cit>  using hwler.
with only one run of hwler, pella and masuda were able to obtain a result similar to that obtained by rosenberg et al.
with multiple runs of structure.
here we used this data set to demonstrate the extent to which the choice of λ affects the behavior of the dp prior method.
bull data set
this real data set consists of  <dig> bulls maintained in japan.
these animals were born between  <dig> and  <dig> and had been genotyped for parentage testing.
they included some half sibs but excluded full sibs.
the number of microsatellites was  <dig> and the mean observed heterozygosity was  <dig> .
this data set was used to demonstrate how unbalanced sample sizes among populations affect the behavior of clustering methods.
this data set is provided in additional file  <dig>
RESULTS
evaluation of the sams sampler
the sams sampler was implemented in hwler to improve the mixing of mcmc algorithms  <cit> .
however, the effectiveness of this sampler in population structure analysis was unknown.
therefore, we evaluated the sams sampler with our program, dpart, using simulated data sets generated by simulation method  <dig> and compared it with structurama.
the number of populations  was  <dig>   <dig>  and  <dig>  and the number of individuals per population was  <dig>  for dpart, λ was  <dig> for all loci and α was set such that the expected k value was equal to the true k value.
furthermore, we analyzed the data sets using only the gibbs sampler by dpart for comparison.
although both dpart with the gibbs and sams samplers and structurama were superior to gibbs sampling alone, the former provided more accurate results than the latter, and the difference became prominent as k increased .
these results showed that the sams sampler was effective in population structure analysis.
note that the sams sampler has an additional advantage with regard to calculation time because one attempt of this sampler is faster than one scan with the gibbs sampler, and unlike mcmcmc, the sams sampler does not need multiple chains.
average , which is the partition distance between the true and inferred partition, the number of data sets in which  was  <dig>  or less  in  <dig> simulated data sets, and the average k values in the inferred partitions  are shown.
the number of individuals in each population was  <dig>  the number of loci was  <dig> for microsatellites and  <dig> and for snps.
the migration rate was  <dig>  for microsatellites.
mc <dig> indicates mcmcmc.
effect of λ
in the dp prior method, allele frequencies of populations are assumed to be drawn from the dirichlet distribution dirichlet .
because the variance of the distribution decreases as λ increases, the frequencies approach uniformity among alleles as λ increases.
thus, it was expected that λ would affect clustering behavior depending on the uniformity of frequencies among alleles, i.e., the preferable values of λ would vary depending on uniformity.
we examined this hypothesis by analyzing data sets that were generated with simulation method  <dig>  in this simulation method, the level of uniformity of frequencies among alleles could be determined by varying ancestral allele frequencies.
three scenarios were used for both microsatellites and snps.
in one scenario, frequencies among alleles were relatively uniform at all loci .
in another scenario, frequencies were skewed at all loci , and in the last scenario, frequencies were relatively uniform for half of the loci and skewed for the other half.
we analyzed these data sets with dpart using the gibbs and sams samplers under different settings of λ.
as expected, the results showed that the preferable values of λ varied depending on the uniformity of frequencies among alleles .
when the frequencies were closer to uniformity, higher values were preferable.
in addition, when a data set included loci that differed significantly in the uniformity of frequencies among alleles, analysis with a single λ value was less accurate than that with inferring λ for each locus.
these results suggest that to maximize the performance of the dp prior method, the λ value should to be determined properly for each locus according to the allele frequencies.
although inferring λ for each locus is a solution for this problem, the specification of λ in this manner tended to make inferences less accurate than assuming single values for all loci when the uniformities of allele frequencies were relatively equal among loci.
thus, we recommend that analyses be repeated under two different assumptions, a single λ value for all loci and a unique value for each locus.
in each assumption, inferring λ would be useful.
average , the number of data sets in which  was  <dig>  or less , and the average k values  are shown.
the number of populations was  <dig> and the number of individuals per population was  <dig>  vectors in parentheses indicate ancestral allele frequencies.
"mean major allele frequency" indicates the mean values of major allele frequencies in the data sets.
jl is the number of observed alleles.
"inferred " indicates that a unique λ value was inferred for each locus, and "inferred " indicates that a single value was inferred for all loci.
the number of populations was  <dig> and the number of individuals per population was  <dig>
analysis of the chicken data set
this data set, representing  <dig> chickens of  <dig> european breeds, was analyzed previously using structure and hwler  <cit> .
rosenberg et al.
<cit>  examined k values ranging from  <dig> to  <dig> using structure and selected  <dig> as the proposed k value according to likelihood.
then,  <dig> runs of structure were performed assuming that k was  <dig>  the authors reported that most breeds could be distinguished from each other, but breeds  <dig> and  <dig> shared a cluster in all runs.
pella and masuda  <cit>  analyzed the data set with hwler assuming that k was  <dig> and detected  <dig> clusters.
similar to that in structure analysis, breeds  <dig> and  <dig> were not distinguished.
hwler divided breed  <dig> into two clusters that included  <dig> and  <dig> individuals and detected three additional clusters that included  <dig>   <dig>  and  <dig> individuals in breed  <dig>  the three individuals in breed  <dig> were sampled from a flock of zoo animals, which were reported to be frequently assigned incorrectly in structure analyses.
if the λ value is determined according to the number of alleles at each locus, similar to that in hwler, λ ranges from  <dig>  to  <dig>  and the mean value is  <dig>  because the number of alleles at each locus ranges from  <dig> to  <dig>  the average major allele frequency across all loci and breeds was  <dig> , indicating that the allele frequencies were relatively skewed.
we analyzed this data set with dpart using the gibbs and sams samplers, varying the λ value to demonstrate the extent to which λ affects the behavior of the dp prior method.
α was set to  <dig> , resulting in the expected k value of  <dig> .
the results are summarized in table  <dig>  when a unique λ value was inferred for each locus and λ was  <dig> , the result was the same as that obtained with hwler.
the λ value inferred for each locus ranged from  <dig>  to  <dig>  and the average across loci was  <dig> .
on the other hand, when a single value was inferred for all loci, breed  <dig> was not divided into two clusters.
the λ value inferred for all loci was  <dig> .
as λ was increased from  <dig>  to  <dig>  the number of detected clusters decreased.
thus, the finest partition was obtained when λ was , inferred for each locus, and set to  <dig> .
this indicates that for this data set, the specification of λ according to the number of alleles was appropriate.
this is probably because of the fact that the allele frequencies at any loci were relatively skewed.
these results also indicate that inferring λ was actually useful in empirical data sets.
comparison among dpart, fum, and fcm
we compared dpart, fum, and fcm, focusing on cases with unbalanced sample sizes among populations.
data sets were generated by simulation method  <dig>  first, we assumed situations in which k =  <dig>  the size of the smaller subset was fixed at  <dig>  the number of microsatellites was  <dig>  and the migration rate was  <dig> .
the size of the larger subset was  <dig>   <dig>   <dig>  and  <dig>  hereafter, the sample size is denoted as n , n , and so on.
dpart was used with the gibbs and sams samplers.
the true k value was used for fum and fcm.
the results showed that dpart was insensitive to unbalanced sample sizes .
in contrast, fcm was the most sensitive of the three programs to unbalance.
fum was less sensitive than fcm, but was inferior to dpart when the sizes were n .
the difference among methods was most prominent at n , but it decreased or disappeared when the migration rate decreased or the number of loci increased.
average , the number of data sets in which  was  <dig>  or less , and the average k values  are shown.
nl and m indicate the number of loci and the migration rate, respectively.
n  denotes sample sizes.
next, we increased the size of the smaller subset to  <dig> to create a moderate unbalance.
the number of loci was decreased to  <dig> and the migration rate was increased to  <dig>  in order to compare the differences more clearly.
again, fcm was found to be most sensitive to unbalance .
although the performance of dpart decreased slightly as the sample size became increasingly unbalanced, dpart provided the highest number of data sets in which  was  <dig>  or less at n .
the difference among methods decreased or disappeared when the migration rate decreased or the number of loci increased.
furthermore, we examined whether the number of minor subsets affected performance of these methods.
we compared the methods in situations in which the sample sizes were n , n , n , and n .
the migration rate was  <dig>  in each situation.
when multiple minor subsets were included in the data sets, i.e., at n  and n , dpart outperformed fum and fcm, suggesting that fum and fcm were severely affected by multiple minor subsets .
in contrast, when only one minor subset was included in the data set, i.e., at n  and n , the effect of the minor subset was relatively small.
the migration rate was  <dig>
in these analyses, fum and fcm often assigned individuals to clusters such that the sizes of the clusters were uniform, resulting in the failure of analysis.
for example, at k =  <dig> or k =  <dig> with single minor subsets, the smaller subsets tended to absorb members in the larger subsets.
at k =  <dig> with multiple minor subsets, fum and fcm often failed to distinguish the two smaller subsets and divided the larger subset into three clusters such that their sizes were uniform.
when the sizes were n  and the number of loci was  <dig>  fum frequently generated an empty cluster.
in such cases, fum detected only two clusters, consisting of two smaller subsets and the larger subset.
empty clusters were not observed when the number of loci was  <dig>  because the larger subset had no population structure, fum probably detected the larger subset correctly as one cluster because of the increase in the number of loci.
however, since fum failed to divide the smaller subsets, only two clusters were detected by fum.
we did not observe this phenomenon in analyses with fcm; the difference in performance between fum and fcm may be relevant in this situation.
analysis of the bull data set
the bull data set, representing  <dig> bulls genotyped with  <dig> microsatellites, was analyzed by dpart, fum, fcm, and structure to demonstrate how unbalanced sample sizes among populations affect the results of these methods.
in dpart, the α value was  <dig> , resulting in the expected k value of  <dig>  a single value of λ was inferred for all loci.
as a result, we obtained two partitions, a partition with five clusters denoted as clusters a to e and a partition with four clusters.
when a unique λ value was inferred for each locus, only the partition with four clusters was obtained.
in this partition, cluster c was absorbed in clusters d and e. thus, we created three data sets, each including clusters c and d, clusters c and e, or clusters d and e, and reanalyzed them with dpart.
λ was inferred for each locus.
because clusters in each data set differed from each other, we concluded that the bull data set included five clusters.
the dendrogram was generated on the basis of co-assignment probabilities calculated from  <dig>  mcmc samples .
the data set was also analyzed using fum, fcm, and structure with k =  <dig>  as shown in figure  <dig>  although the result obtained with fum was consistent with that obtained with dpart , the results obtained with fcm and structure were not consistent with those obtained with fum and dpart .
in the results of fcm and structure, the smallest cluster in the results of dpart and fum, i.e., cluster a , absorbed the members of clusters b , c , and e .
in addition, the moderate cluster, cluster e, also absorbed the members of cluster d .
on the other hand, the second smallest cluster in the results of dpart and fum, cluster b, absorbed very few members of the larger clusters in the results of fcm and structure.
the pairwise fst values between the clusters detected by dpart are shown in table  <dig>  we consider the interpretation of these inconsistent results in the discussion section.
discussion
the bayesian method based on the dp prior can infer the number of populations  and assign individuals, whereas the selection of the appropriate k value is often problematic when methods that run under a predefined k value are used.
we examined the properties of this method to provide information that will be useful in practice.
we showed that the sams sampler, which assigns multiple individuals simultaneously, was effective for the inference of population structures.
because sams sampling was faster than mcmcmc techniques, the implementation of this sampling technique may be especially useful for large data sets.
we also showed that a hyperparameter, named λ, which defines the prior distribution of allele frequencies, affected the performance of the method and its specification was important.
this problem could be resolved by considering λ a variable.
furthermore, we demonstrated that the dp prior method was suitable for data sets having unbalanced sample sizes among populations, whereas methods that implement structure-like algorithms were sensitive to unbalance.
in particular, we found that the allele frequencies correlated model was the most sensitive.
our results showed that both the sams sampler and mcmcmc were effective in improving the mixing of mcmc algorithms; however, the sams sampler was more effective than mcmcmc.
we implemented the sams sampler at a frequency in which four iterations with the sams sampler were followed by one iteration with the gibbs sampler .
although we examined other frequencies, such as two sams and five gibbs, one sams and one gibbs, and five sams and two gibbs, the results were almost the same .
we selected this frequency simply to shorten the run time, because one attempt of the sams sampler is faster than one scan of the gibbs sampler.
however, an accurate inference of the posterior distribution is hardly possible with the sams sampler alone and the gibbs sampler is necessary.
the acceptance rate of the proposed states was extremely low , suggesting the difficulty of proposing new partitions.
structurama implementing mcmcmc can probably increase its performance by increasing the number of chains or adjusting the temperature parameter to obtain the appropriate exchange rates among chains.
however, calculation time increases as the number of chains increases and multiple runs may have to be performed to find appropriate values for the temperature parameter.
therefore, we concluded that the sams sampler was more useful in practice and implemented this sampler in dpart to improve the mixing of mcmc algorithms.
we showed that the choice of λ values affected the behavior of the dp prior method in the simulated and real data sets.
since the preferable λ value varies depending on the uniformity of frequencies among alleles, it is desirable that the λ value at each locus is determined according to the allele frequencies.
because hwler set λ to  for any loci, its performance is probably inadequate for some data sets, although it implements the sams sampler.
structurama does not state how it specifies λ, but it appears to fix λ for any loci at a certain value.
inferring a unique λ value for each locus is a method of specifying the parameter for each locus.
however, this approach resulted in less accurate inferences than assuming a single value for all loci when the levels of uniformity of frequencies among alleles were relatively similar among loci.
we speculate that increasing the number of hyperparameters to be estimated may make inferences unstable.
in contrast, assuming a single value for all loci was less accurate in data sets in which the levels of uniformity differed significantly among loci.
therefore, it is difficult to state which assumption is more suitable for the data set of interest.
in the chicken data set, the unique value assumption  and the single value assumption  provided the same partition.
however, analysis in which a single value was assumed and inferred for all loci gave a slightly rougher partition.
in the bull data set, the single value assumption  gave a finer partition than the unique value assumption .
in general, if the loci included in the data have not been evaluated well with regard to polymorphism, some loci may be much less polymorphic than others, and thus, the allele frequencies at these loci will be skewed.
for such data sets, the unique value assumption will be suitable.
on the other hand, if the loci included in the data are selected from a large number of candidate loci, they will be highly polymorphic, and thus, the allele frequencies at these loci are expected to be close to uniformity.
therefore, the single value assumption will be suitable for such data sets.
we speculate that the chicken data may be closer to the former case and, in such a case, even if an appropriate single value that leads to an accurate inference exists, the inference of such a value may be difficult.
in contrast, the bull data set may be closer to the latter case because the loci included in the data set had been selected from a large number of candidates for parentage testing.
therefore, the single value assumption may be preferable.
our results showed that the behavior of the dp prior method depends on the selection of λ. we speculated that integration out of allele frequencies involves this dependency to some extent.
thus, although we have not examined this speculation, the dependency may decrease by inferring allele frequencies.
however, this will increase the calculation time and thus will not be suitable for large data sets.
we found that unbalanced sample sizes among populations affect the behavior of dpart, fum, and fcm.
fum and fcm were found to be sensitive to unbalanced sample sizes, and their performances were substantially affected, particularly by the presence of multiple minor subsets.
the reason why dpart is suitable for unbalanced sample sizes is probably its prior assumptions about the assignment of individuals.
as seen in eq.
, in the algorithm implementing the dp prior, clusters absorb individuals with higher probability as sample sizes increase, i.e., the "rich get richer" phenomenon occurs.
thus, this algorithm is suitable for data sets with unbalanced sample sizes among populations.
on the other hand, the algorithms implemented in fum, fcm, and structure assume that each population can contribute to the data set with equal probability.
thus, these algorithms are suitable for data sets in which sample sizes are uniform among populations.
when these methods are used for data sets with unbalanced sample sizes, they tend to cluster individuals such that the sizes become uniform among clusters, as observed in our experiments.
the extent of sensitivity varied depending on the number of loci and the migration rates when the level of unbalance was fixed.
thus, if genetic differences among populations are small and sample sizes are unbalanced, then the number of loci needed by these methods to correctly detect population structures is higher than that needed by dpart.
although we compared these programs using only microsatellite data, the differences in behavior among them will not vary if snp data are analyzed because the differences will be due to differences among the prior assumptions and will not depend on the number of alleles.
the sensitivity of fum and fcm to unbalance can probably be resolved by adding parameters that represent the mixing weights to the algorithms.
as described by neal  <cit> , when the prior distribution of the mixing weights is assumed to be a dirichlet distribution such as  where k is the assumed number of populations and these weights are integrated out, we can obtain a predictive distribution for the assignment of the ith individual,
which is similar to eq.
, but now in the parametric framework.
fum, fcm, and other structure-like algorithms will be able to deal with unbalanced sample sizes by implementing this prior distribution.
fcm, which implements the allele frequency correlated model, was shown to be more sensitive than the other programs to unbalanced sample sizes.
a refutation of this may be that comparison among methods was performed for simulated data sets generated using the isolated population models that do not accord with the assumption of fcm.
however, fcm was most sensitive to unbalance in the data sets that were generated by simulation method  <dig>  which was based on the assumption of fcm .
we speculated that correlated allele frequencies are involved in the sensitivity of fcm.
in the allele frequencies correlated model, the conditional posterior distribution for allele frequencies of population j is given as follows:
where yjlh is the number of copies of allele h at locus l in individuals assigned to population j, palh is the ancestral frequency of allele h at locus l, jl is the number of alleles at locus l, and fj is  where fj is the drift parameter for population j. thus, as the sample size for a population decreases, the effect of ancestral allele frequencies on the inference of allele frequencies for this population increases.
on the other hand, ancestral allele frequencies are inferred from information on allele frequencies of populations and drift parameters.
therefore, if unbalanced sample sizes are present, ancestral allele frequencies will be affected more strongly by the inferred allele frequencies for major subsets because those for minor subsets are substantially affected by ancestral allele frequencies themselves, and are thus, less informative for the inference of ancestral allele frequencies.
consequently, the inferred allele frequencies for minor subsets will be affected by those for major subsets and will approach them.
therefore, minor subsets may be prone to absorbing members of major subsets.
we interpreted the inconsistency found in the bull data set on the basis of the knowledge obtained from simulations.
if the results obtained with dpart and fum are correct, this inconsistency can be explained by the sensitivity to the unbalance of the allele frequencies correlated model.
fcm and structure probably failed to detect the smallest cluster  because of their sensitivity.
the failure of cluster a to absorb the members of cluster d in the results of fcm and structure was due to the high level of differentiation between these clusters .
although the unbalance between clusters d and e was moderate , fcm and structure also presumably failed to distinguish these clusters because of the relatively low pairwise fst between these clusters .
when the data set that included only clusters d and e was analyzed by fcm, this program also failed to distinguish the two clusters .
on the other hand, because cluster b was well differentiated from larger clusters , every program was able to detect the cluster.
the presence of multiple minor clusters was also considered to reduce the performance of fcm and structure.
in fact, when the two data sets that included only clusters a and b and clusters a and e were generated and analyzed by fcm, this program was able to distinguish the two clusters in each data set .
although we admit that we have not proved that the results of dpart and fum are correct, we believe that our interpretation is appropriate because it can clearly explain the inconsistency.
for the bull data set,  <dig> runs were performed with fcm and structure.
the results of the  <dig> runs of each program were almost similar, and we considered them to be incorrect.
however, in the simulated data sets with unbalanced sample sizes among populations, we occasionally observed that these programs or fum provided both correct and incorrect results when multiple runs were performed .
thus, when the ad hoc statistic proposed by evanno et al.
<cit>  is used to select the true k value, this phenomenon possibly confuses the selection because it increases the standard deviation of likelihood at the true k value.
the effect of sample sizes of populations on the performance of clustering programs was addressed in some studies .
however, the effect of unbalanced sample sizes has been overlooked, and simulation studies for comparison of clustering methods usually assume uniform sample sizes among populations  <cit> .
because our results showed that sensitivity to unbalance in sizes varied among the methods, we recommend that comparative studies consider the effect of unbalance during analyses.
through analyses of simulated data sets by dpart, we observed overestimation of k caused by small clusters that included only one or two individuals.
this phenomenon increased the average k values and slightly affected the average .
as discussed above, the dp prior method can efficiently detect minor subsets because of the "rich get richer" phenomenon.
thus, we speculated that the overestimation was due to the fact that dpart detected individuals that were slightly distanced from the other members because even in simulated data sets, individuals harboring rare genotypes can be generated with low probabilities.
this interpretation was supported by phylogenetic analysis based on genetic distances between individuals .
in addition, this may be supported by the fact that overestimation became prominent as the number of loci, i.e., the power to detect population structures, increased .
therefore, such small clusters are interpreted as overestimations in simulation studies, but will provide useful information in empirical studies because they indicate the presence of genetic discontinuity in the data sets.
CONCLUSIONS
this study characterized the bayesian method of implementing the dp prior and introduced a program, named dpart, in order to infer population structures more accurately than preceding programs based on the dp prior.
first, we showed that the sams sampler, which is a technique for improving the mixing of mcmc algorithms, was effective for population structure analysis.
implementation of this sampler was useful with regard to the accuracy of inference and computing time.
second, we showed that a hyperparameter for the prior distribution of allele frequencies affected the behavior of the dp prior method.
appropriate values can be specified by inferring this parameter.
third, the dp prior method was shown to be suitable for analysis of data sets with unbalanced sample sizes among populations.
in contrast, methods that implement structure-like algorithms were shown to be suitable for data sets with uniform sample sizes among populations, but not for data sets with unbalanced sample sizes.
because these differences can yield inconsistent results among methods, we recommend using these methods concurrently.
when the results obtained are inconsistent among methods, considering the effect of unbalanced sample sizes may be a key to interpreting the inconsistency.
competing interests
the authors declare that they have no competing interests.
authors' contributions
ao conceived the study, performed all statistical analyses, and drafted the manuscript.
mn and mm collaborated on the algorithm design and the manuscript.
all authors have read and approved the final manuscript.
endnote a - gibbs sampler
one scan of the gibbs sampler consists of the following steps.
step  <dig>  remove the ith individual from η.
step  <dig>  assign i to existing populations according to probabilities
where  <dig> ≤ j ≤ k, or a new population according to probability
step  <dig>  update η.
step  <dig>  repeat steps 1- <dig> for all i ∈ { <dig> ...,n}.
endnote b - sams sampler
the algorithm of the sams sampler is as follows.
step  <dig>  select two individuals, i and j, at random.
step  <dig>  if i and j belong to the same population s, remove i and j from s and form two singletons, si = {i} and sj = {j}.
if not, proceed to step  <dig>
step  <dig>  assign the individuals remaining in s to si or sj.
the order of the assignment is randomly determined.
the kth individual is assigned to si with probability
otherwise, add the individual to sj.
step  <dig>  the proposed partition η' is accepted with mh probability
where q and q are the transition probabilities.
using eq.
, the ratio of likelihoods in eq.
can be written as follows:
using eq.
, the ratio of prior probabilities in eq.
can be written as follows:
q is  <dig> and q results from the products of probabilities of the assignment given in eq.
.
step  <dig>  if i and j belong to different populations, say si and sj, propose a new population s by merging si and sj.
the proposed partition η' is accepted with the probability given in eq.
.
the ratio of likelihoods can be written as follows:
and the ratio of prior probabilities is given as follows:
q is  <dig> and q is computed by splitting s into si and sj in a randomly determined order.
endnote c - prior distribution of k
prior distribution for k can be inferred using the following monte carlo procedure.
step  <dig>  let the first individual belong to the first population and let k =  <dig>
step  <dig>  assign individual i = { <dig> ,...,n} to existing or new populations with the probabilities noted in eq.
.
step  <dig>  record k after the assignment of the last individual.
step  <dig>  repeat steps 1- <dig> for sufficient cycles to infer the distribution of k .
endnote d - mean partition
the algorithm for calculating the mean partition described by huelsenbeck and andolfatto  <cit>  is as follows.
step  <dig>  pick a sampled partition as the initial state of the mean partition and calculate d, which is the sum of the partition distances between the mean partition and every sampled partition.
step  <dig>  pick an individual i in the mean partition.
propose new mean partitions by moving i to other populations in the mean partition and to a new partition.
calculate the sum of partition distances between each proposed mean partition and each sampled partition, which is denoted as d'.
step  <dig>  let d'min denote the minimum value of d's.
if d'min <d, the corresponding proposed mean partition is accepted and d is replaced by d'min.
step  <dig>  repeat step  <dig> and  <dig> for i = { <dig> ,...,n}.
step  <dig>  repeat steps  <dig>   <dig>  and  <dig> until d stops decreasing.
supplementary material
additional file 1
an r function for agglomerative hierarchical clustering.
this performs agglomerative hierarchical clustering based on the results of dpart.
click here for file
additional file 2
executable file of dpart.
this file will work on a windows platform.
click here for file
additional file 3
parameters.
this file defines the parameters for dpart.
click here for file
additional file 4
input file names.
this file defines the input file names for dpart.
click here for file
additional file 5
source code.
this is the source code of dpart.
click here for file
additional file 6
header file.
this is the header file of dpart.
click here for file
additional file 7
manual.
this is the manual for dpart.
click here for file
additional file 8
bull data set.
this file contains the genotypes of the bulls used in this study.
this file is also an example of input data for dpart.
click here for file
acknowledgements
we would like to thank two anonymous referees for their helpful comments.
we also would like to thank the technical staff at maebashi institute of animal science for genotyping bulls.

 in population genetics, inference of population structures is important for various purposes such as assessment of genetic diversity, detection of genetic discontinuities in natural wildlife habitats, and correction for stratification in association studies. to infer population structures without prior knowledge about the population, various statistical approaches using neutral molecular markers have been proposed. bayesian approaches using markov chain monte carlo methods have been widely used to infer population structures since pritchard et al. proposed the bayesian clustering algorithms implemented in the well-known program structure. this program can infer the assignment of individuals to populations or the admixture proportions of individuals for a given number of populations. researchers have extended bayesian algorithms for various purposes such as to take advantage of spatial information, estimate inbreeding coefficients, allow for allele mutations, and infer k values. pella and masuda used a dirichlet process to infer k values. dp is a stochastic process that was proposed by ferguson to treat nonparametric problems in bayesian frameworks. the merit of using dp to infer k is that k can take any value between and the number of individuals, and thus, few assumptions about k are needed for inference. pella and masuda considered k and the assignment of individuals to populations as random variables using dp as a prior distribution for k and allele frequencies unique to populations. huelsenbeck and andolfatto also used the dp prior for the inference of population structures, and reich and bondell proposed a clustering algorithm using the dp prior, which can incorporate spatial information. besides the inference of population structures, dp priors have been used to infer the number of ancestral haplotype blocks, to model nonsynonymous/synonymous rate ratios, and to model the selfing rates of individuals. to date, two clustering programs that implement the dp prior have been provided, hwler and structurama. both programs implement the gibbs sampling procedure to infer the posterior distribution. these programs differ in their approach to improve the mixing of mcmc algorithms. hwler implements the sequentially-allocated merge-split sampler, which moves multiple observations simultaneously, and structurama implements the metropolis-coupled mcmc technique, which runs multiple chains, some of which are closer to a uniform distribution than the target distribution, and attempts to swap states among chains. although hwler and structurama are useful and have been used in some recent studies, their application to real data sets has been less common compared with that of structure. this may be because the properties of these methods have not been investigated in detail. when results obtained with different methods are inconsistent, it is difficult to interpret the results. for example, hwler and structurama may provide inconsistent results. in addition, although hwler detects three populations, the result obtained with structure under the assumption that k is may be different from that obtained with hwler. the purpose of this study was to characterize the bayesian method based on the dp prior and to provide information that will be useful in practice. first, we evaluated the sams sampler because its effectiveness has not been examined. second, we investigated the effect of a hyperparameter, named λ, which defines the prior distribution of allele frequencies, on the performance of this method. as described by pella and masuda, hwler set λ to where jl is the number of alleles at locus l for any loci. however, this specification resulted in inaccurate inferences in some cases. third, we compared the dp prior method with other bayesian methods that infer the assignment of individuals for a given k value. we focused primarily on the effect of unbalanced sample sizes among populations on the behavior of these methods because unbalanced sample sizes were found to affect the behavior of these methods in our preliminary study. methods assumption and goal we assumed that a sampled data set consists of diploid individuals derived from an unknown number of populations, which are defined by unique values of allele frequencies. although our novel program, named dpart, can also analyze haploid populations, we assumed only diploid individuals in this study for simplicity. our goal was to infer the number of populations and assign individuals to their populations based on their genotypes. the estimation of admixture proportions was not of interest. dp prior the mathematical details of dp can be found in studies conducted by ferguson, antoniak, and neal. here we provided a brief description of the dp prior. approximately and intuitively, dp can be seen as a stochastic process that converts a continuous distribution to a discrete distribution. we assumed that observations are generated from the following model. let g denote the prior distribution of allele frequencies at the locus, which is continuous. dp is defined by g and the concentration parameter α. dp divides g into a number of classes, each of which is represented by a single point and yields a discrete probability distribution, g. the number of classes, which is determined by the number of observations and α, can be infinite. a vector of allele frequencies for each genotype is drawn from g instead of g and genotypes are drawn from the corresponding allele frequency vectors. this model can be written as follows: where xi is the genotype of individual i, ϕi is the allele frequency vector for genotype i, and dp is dp. this model is known as the dp mixture model. because g is discrete, allele frequency vectors for some genotypes may have values in common, i.e., these genotypes can be seen as members of the same population, which is characterized by the shared allele frequency vectors. parameters by integrating out g, we can obtain a simpler representation of the model. when g is integrated out, the predictive distribution for the allele frequency vector ϕi conditional on {ϕ..., ϕi-1} can be written as follows: where nj is the number of allele frequency vectors that share values with ϕj. this sequence of predictive distributions is known as a polya urn scheme. to represent the clustering property of the model more explicitly, we let {ϕ...,ϕk} denote allele frequency vectors unique to populations {...,k}. in addition, we introduced parameters that represent the partition of individuals {...,n}, according to the parameterization of dahl. we let η = {s...,sk} define a partition for {...,n} such that, si ∩ sj = ∅ for all i ≠ j, and si ≠ ∅ for all i. eq. includes the following two types of prior information: where η<i is η before assigning the ith individual and |sj| is the number of individuals included in sj, which is the prior for the number of populations and assignment of individuals, and which is the prior for allele frequencies for each population. now, we can rewrite eq. as follows: where x = {x...,xn} is a vector of the genotypes of the individuals {...,n}, i{i ∈ sj} = if i ∈ sj, and otherwise and where Γ is the gamma function. note that eq. results from the products of eq.. integration of allele frequencies although some clustering methods infer allele frequencies, allele frequencies can be analytically integrated out as in hwler and structurama. thus, our goal is to infer not f but f. by bayes' theorem, this can be written as follows: where where denotes the genotypes of individuals included in sj and is the posterior distribution of ϕj updated from g on the basis of the genotypes of individuals assigned to sj preceding the ith individual. g is assumed to be a flat dirichlet distribution. the prior probability density of ϕ is given as follows: where l is the number of loci, jl is the number of alleles at locus l, λl is the hyperparameter for locus l, and ϕjlh is the frequency of allele h at locus l of population j. rannala and mountain provided the following equation: where yjlh is the number of copies of allele h at locus l in individuals assigned to sj preceding the ith individual and yjl = ∑yjlh. l in eq. can be rewritten as follows: where δ = if the genotype of individual i at locus l is heterozygous and if not is and nilh is the number of copies of allele h at locus l in individual i. now the integral in eq. can be solved analytically as follows: gibbs sampler neal and maceachern introduced the gibbs sampling procedure, which integrates out model parameters. let η-i denote η when only the ith individual is removed. the conditional probability of the assignment of individual i is given as follows: where b is the normalizing constant. because individuals are exchangeable, one can treat every individual as the last observation. the integrals in eq. can be solved analytically. the details of gibbs sampling is described in endnote a. in our experiments, one mcmc iteration consisted of one scan of the gibbs sampler. after burn-in, η is sampled in every predefined number of mcmc iterations. sams sampler because the gibbs sampler moves only one observation, a large change hardly ever occurs, and tends to stay at a local mode. a remedy is to move multiple observations simultaneously. jain and neal proposed a split-merge sampler, which proposes a new partition by splitting or merging components and accepts it with metropolis-hastings probability. when it splits a component, the proposed assignment of observations in the component is generated by repeating gibbs sampling only for these observations. this method can produce more acceptable proposals than the random assignment of observations. dahl improved the method in terms of efficiency and proposed an alternative sampler, the sams sampler, which generates proposals by one cycle of allocation of observations instead of repeating gibbs sampling. the algorithm of the sams sampling is described in endnote b. in this study, the sams sampling was performed once in an mcmc iteration. four iterations with the sams sampler were followed by one iteration with the gibbs sampler. this cycle was repeated until the end of the iterations. inference of the hyperparameter λl, which is the hyperparameter for the prior distribution of allele frequencies at locus l, can be inferred by the mh method. the joint posterior probability of η and λ = {λ...,λl} can be written as follows: when the prior distribution for λl is assumed to be uniform, equations in the gibbs and sams samplers can be used without modification. in this study, the prior for λl was assumed to be uniform, u. the proposal of λl, λl', was generated from the normal distribution n, where δ is an arbitrary value, e.g.,. the proposal was accepted with the probability: where λ-l is λ when λl is removed. we can also assume that each locus shares a single value, λa. in this case the probability becomes the following: prior distribution of the number of populations as seen in eq. or, the prior distribution of k depends on α and the number of individuals. we can infer this prior distribution by the monte carlo method. for example, when n = and α =, the expected k value is approximately. when α increases to, k increases to approximately. although it is possible to infer α as well as λ, in this study, we fixed α through the mcmc iterations. the effect of α has been thoroughly examined by huelsenbeck and andolfatto. the authors reported that the misspecification of α could affect the results, especially when the number of loci was small. summarizing the sampled partitions two methods were used to summarize the partitions that were sampled from the posterior distribution of η. for simulated data sets, we used the mean partition proposed by huelsenbeck and andolfatto, which is defined as follows: where v is the number of sampled partitions and d is the partition distance between the sampled partition ηi and mean partition. the partition distance between ηi and u is defined as the minimum number of individuals that must be removed from both ηi and u such that the two partitions are the same. the partition distance was calculated using eq. of konovalov et al. to obtain the cost matrix, and we solved this using the hungarian algorithm. the algorithm for calculating the mean partition is described in endnote d. this statistic is useful for evaluating the posterior distribution of η automatically. we divided the partition distances by the number of individuals. thus, the partition distances ranged from to for real data sets, we performed agglomerative hierarchical clustering on the basis of co-assignment probabilities. this method was introduced by dawson and belkhir. briefly, after all samples were obtained, the co-assignment probabilities for all individual pairs were calculated from the sampled partitions. next, complete linkage clustering was performed on the basis of the probabilities. we used the hclust function of the stats package of r for complete linkage clustering. an example r code is provided in additional file note that neither the mean partition nor hierarchical clustering are affected by the label-switching problem that often emerges in analyses using methods implementing the dp prior. programs and mcmc parameters we wrote the clustering programs in c. the program implementing the dp prior is referred to as dpart. in addition to dpart, we also produced two clustering programs that infer only the assignment of individuals under a given k value. these programs are based on algorithms proposed by pritchard et al. and falush et al.. in these algorithms, the prior probability that an individual belongs to each population is equal among the populations, and individuals are assigned to populations according to likelihood. these algorithms infer the allele frequencies of populations that are integrated out in dpart in addition to the assignment of individuals. one program assumes that there is no correlation among the allele frequencies of populations, and the other assumes a correlation. thus, the former is equivalent to the "no admixture and no f model" setting of structure and the latter is equivalent to the "no admixture and f model" setting. these programs are referred to as fixed k and uncorrelated model and fixed k and correlated model, respectively. structure infers admixture proportions of individuals until burn-in is half completed even when no admixture option is selected. because this setting helps the program to avoid the generation of an empty cluster, which has no individuals, we followed the setting. in fcm, the prior distribution of the drift parameter was assumed to be u. the programs provided in this study are summarized in table dpart is provided in additional files and fum and fcm can be reproduced by structure. structurama was used for simulated data sets. the expected k value was set equal to the true k value. four chains were run simultaneously and the temperature was set to. we also used structure ver. with its default setting, "admixture and f model," for the bull data set because it is a widely used method for the inference of population structures, although the estimation of admixture proportions was not of interest. for comparison between dpart and structurama and examination of the effect of λ, the number of mcmc iterations was, and the first half of the iterations was discarded as burn-in. sampling was performed every iterations. for comparison among dpart, fum, and fcm, the number of iterations was, and the first half of the iterations was discarded. sampling was performed every iterations. for the real data sets, the number of iterations was and the first iterations were discarded. sampling was performed every iterations. simulation method 1 this simulation method was used to compare among clustering methods. easypop was used to simulate populations under the wright-fisher model. the island model was assumed. the number of populations was and the number of individuals per population was. the number of generations, number of loci, mutation rate, and possible number of alleles were,, and respectively, for microsatellites, and, × 10- and respectively, for single nucleotide polymorphisms. new alleles occurred in every possible allelic state with equal probability. diploidy, random mating, one sex, and linkage equilibrium were also assumed. the genotypes of the first generation were generated by randomly assigning alleles from the possible allelic states. migration rates were,, and for microsatellites and for snps. easypop was run five times for each scenario. the resulting data sets are referred to as base population sets. sampling was performed times from each base population set to obtain data sets per scenario. when two, three or four populations were sampled from the base population sets, base populations to be sampled were selected randomly. similarly, when or loci were sampled, the loci to be sampled were also selected randomly. the pairwise fst values between the base populations for each migration rate are summarized in table typically, the observed heterozygosities were approximately for microsatellites and for snps. the pairwise fst and heterozygosity values were calculated using genetix. m indicates the migration rate. simulation method 2 this method was used to investigate the effect of λ. the variance of a dirichlet distribution, dirichlet, decreases as λ increases. when allele frequencies are drawn from the dirichlet distribution, it means that the frequencies closely approach uniformity among alleles as λ increases. thus, it was expected that λ would affect clustering behavior depending on the uniformity of frequencies among alleles. to examine this hypothesis, we devised a simulation method on the basis of the allele frequency correlated model. in this method, allele frequencies were generated from the dirichlet distribution where pa is the allele frequency of a virtual ancestral population, f is the drift parameter, which is equivalent to wright's fst, and ja is the number of ancestral alleles. we can determine any level of uniformity of frequencies among alleles by varying ancestral allele frequencies. we assumed two marker types, microsatellites and snps. for microsatellites, the ancestral allele frequencies were {,,..., } or {,,..., } and f was, and for snps, the ancestral frequencies were {, } or {, } and f was. k was assumed to be and the number of individuals was per population. genotypes were generated from the allele frequencies assuming random mating. if only one allele was observed at a locus, then that locus was excluded. one hundred data sets were generated for each scenario. comparison among methods for simulated data sets to compare methods for simulated data sets, we used the partition distance between the true and inferred partitions, which was denoted as. we calculated average over the simulated data sets and counted the number of data sets in which was or less. for dpart and structurama, we used the mean partition as the inferred partition and calculated the average k values in the mean partitions. when fum and fcm were used, the inferred partition was determined on the basis of the probabilities of assigning individuals to populations, which were calculated from the sampled partitions. although fum and fcm may be evaluated with the mean partition, we used the partition that was based on the probabilities of assignment because this is computationally more feasible than the mean partition and label-switching occurs rarely in structure-like algorithms as indicated by pritchard et al.. both approaches provided almost the same partitions in our preliminary study. note that can not be used for evaluation without modification when unbalanced sample sizes are present among populations. suppose that two populations are included in a data set. if the sample sizes are uniform between the populations and an analysis fails to detect any population structures, i.e., all individuals are assigned to one cluster, is. however, if the sample sizes are and and an analysis fails similarly, decreases to. thus, in such cases, we calculated as if an individual in the smaller subset was equivalent to individuals in the larger subset. for example, if the sample sizes are and and only one individual in the smaller subset is incorrectly assigned to the larger subset, is but not. chicken data set this data set represents chickens of european breeds, which were genotyped for microsatellites by rosenberg et al.. this data set was previously analyzed by pella and masuda using hwler. with only one run of hwler, pella and masuda were able to obtain a result similar to that obtained by rosenberg et al. with multiple runs of structure. here we used this data set to demonstrate the extent to which the choice of λ affects the behavior of the dp prior method. bull data set this real data set consists of bulls maintained in japan. these animals were born between and and had been genotyped for parentage testing. they included some half sibs but excluded full sibs. the number of microsatellites was and the mean observed heterozygosity was. this data set was used to demonstrate how unbalanced sample sizes among populations affect the behavior of clustering methods. this data set is provided in additional file  evaluation of the sams sampler the sams sampler was implemented in hwler to improve the mixing of mcmc algorithms. however, the effectiveness of this sampler in population structure analysis was unknown. therefore, we evaluated the sams sampler with our program, dpart, using simulated data sets generated by simulation method and compared it with structurama. the number of populations was and and the number of individuals per population was for dpart, λ was for all loci and α was set such that the expected k value was equal to the true k value. furthermore, we analyzed the data sets using only the gibbs sampler by dpart for comparison. although both dpart with the gibbs and sams samplers and structurama were superior to gibbs sampling alone, the former provided more accurate results than the latter, and the difference became prominent as k increased. these results showed that the sams sampler was effective in population structure analysis. note that the sams sampler has an additional advantage with regard to calculation time because one attempt of this sampler is faster than one scan with the gibbs sampler, and unlike mcmcmc, the sams sampler does not need multiple chains. average, which is the partition distance between the true and inferred partition, the number of data sets in which was or less in simulated data sets, and the average k values in the inferred partitions are shown. the number of individuals in each population was the number of loci was for microsatellites and and for snps. the migration rate was for microsatellites. mc indicates mcmcmc. effect of λ in the dp prior method, allele frequencies of populations are assumed to be drawn from the dirichlet distribution dirichlet. because the variance of the distribution decreases as λ increases, the frequencies approach uniformity among alleles as λ increases. thus, it was expected that λ would affect clustering behavior depending on the uniformity of frequencies among alleles, i.e., the preferable values of λ would vary depending on uniformity. we examined this hypothesis by analyzing data sets that were generated with simulation method in this simulation method, the level of uniformity of frequencies among alleles could be determined by varying ancestral allele frequencies. three scenarios were used for both microsatellites and snps. in one scenario, frequencies among alleles were relatively uniform at all loci. in another scenario, frequencies were skewed at all loci, and in the last scenario, frequencies were relatively uniform for half of the loci and skewed for the other half. we analyzed these data sets with dpart using the gibbs and sams samplers under different settings of λ. as expected, the results showed that the preferable values of λ varied depending on the uniformity of frequencies among alleles. when the frequencies were closer to uniformity, higher values were preferable. in addition, when a data set included loci that differed significantly in the uniformity of frequencies among alleles, analysis with a single λ value was less accurate than that with inferring λ for each locus. these results suggest that to maximize the performance of the dp prior method, the λ value should to be determined properly for each locus according to the allele frequencies. although inferring λ for each locus is a solution for this problem, the specification of λ in this manner tended to make inferences less accurate than assuming single values for all loci when the uniformities of allele frequencies were relatively equal among loci. thus, we recommend that analyses be repeated under two different assumptions, a single λ value for all loci and a unique value for each locus. in each assumption, inferring λ would be useful. average, the number of data sets in which was or less, and the average k values are shown. the number of populations was and the number of individuals per population was vectors in parentheses indicate ancestral allele frequencies. "mean major allele frequency" indicates the mean values of major allele frequencies in the data sets. jl is the number of observed alleles. "inferred " indicates that a unique λ value was inferred for each locus, and "inferred " indicates that a single value was inferred for all loci. the number of populations was and the number of individuals per population was analysis of the chicken data set this data set, representing chickens of european breeds, was analyzed previously using structure and hwler. rosenberg et al. examined k values ranging from to using structure and selected as the proposed k value according to likelihood. then, runs of structure were performed assuming that k was the authors reported that most breeds could be distinguished from each other, but breeds and shared a cluster in all runs. pella and masuda analyzed the data set with hwler assuming that k was and detected clusters. similar to that in structure analysis, breeds and were not distinguished. hwler divided breed into two clusters that included and individuals and detected three additional clusters that included and individuals in breed the three individuals in breed were sampled from a flock of zoo animals, which were reported to be frequently assigned incorrectly in structure analyses. if the λ value is determined according to the number of alleles at each locus, similar to that in hwler, λ ranges from to and the mean value is because the number of alleles at each locus ranges from to the average major allele frequency across all loci and breeds was, indicating that the allele frequencies were relatively skewed. we analyzed this data set with dpart using the gibbs and sams samplers, varying the λ value to demonstrate the extent to which λ affects the behavior of the dp prior method. α was set to, resulting in the expected k value of. the results are summarized in table when a unique λ value was inferred for each locus and λ was, the result was the same as that obtained with hwler. the λ value inferred for each locus ranged from to and the average across loci was. on the other hand, when a single value was inferred for all loci, breed was not divided into two clusters. the λ value inferred for all loci was. as λ was increased from to the number of detected clusters decreased. thus, the finest partition was obtained when λ was, inferred for each locus, and set to. this indicates that for this data set, the specification of λ according to the number of alleles was appropriate. this is probably because of the fact that the allele frequencies at any loci were relatively skewed. these results also indicate that inferring λ was actually useful in empirical data sets. comparison among dpart, fum, and fcm we compared dpart, fum, and fcm, focusing on cases with unbalanced sample sizes among populations. data sets were generated by simulation method first, we assumed situations in which k = the size of the smaller subset was fixed at the number of microsatellites was and the migration rate was. the size of the larger subset was and hereafter, the sample size is denoted as n, n, and so on. dpart was used with the gibbs and sams samplers. the true k value was used for fum and fcm. the results showed that dpart was insensitive to unbalanced sample sizes. in contrast, fcm was the most sensitive of the three programs to unbalance. fum was less sensitive than fcm, but was inferior to dpart when the sizes were n. the difference among methods was most prominent at n, but it decreased or disappeared when the migration rate decreased or the number of loci increased. average, the number of data sets in which was or less, and the average k values are shown. nl and m indicate the number of loci and the migration rate, respectively. n denotes sample sizes. next, we increased the size of the smaller subset to to create a moderate unbalance. the number of loci was decreased to and the migration rate was increased to in order to compare the differences more clearly. again, fcm was found to be most sensitive to unbalance. although the performance of dpart decreased slightly as the sample size became increasingly unbalanced, dpart provided the highest number of data sets in which was or less at n. the difference among methods decreased or disappeared when the migration rate decreased or the number of loci increased. furthermore, we examined whether the number of minor subsets affected performance of these methods. we compared the methods in situations in which the sample sizes were n, n, n, and n. the migration rate was in each situation. when multiple minor subsets were included in the data sets, i.e., at n and n, dpart outperformed fum and fcm, suggesting that fum and fcm were severely affected by multiple minor subsets. in contrast, when only one minor subset was included in the data set, i.e., at n and n, the effect of the minor subset was relatively small. the migration rate was in these analyses, fum and fcm often assigned individuals to clusters such that the sizes of the clusters were uniform, resulting in the failure of analysis. for example, at k = or k = with single minor subsets, the smaller subsets tended to absorb members in the larger subsets. at k = with multiple minor subsets, fum and fcm often failed to distinguish the two smaller subsets and divided the larger subset into three clusters such that their sizes were uniform. when the sizes were n and the number of loci was fum frequently generated an empty cluster. in such cases, fum detected only two clusters, consisting of two smaller subsets and the larger subset. empty clusters were not observed when the number of loci was because the larger subset had no population structure, fum probably detected the larger subset correctly as one cluster because of the increase in the number of loci. however, since fum failed to divide the smaller subsets, only two clusters were detected by fum. we did not observe this phenomenon in analyses with fcm; the difference in performance between fum and fcm may be relevant in this situation. analysis of the bull data set the bull data set, representing bulls genotyped with microsatellites, was analyzed by dpart, fum, fcm, and structure to demonstrate how unbalanced sample sizes among populations affect the results of these methods. in dpart, the α value was, resulting in the expected k value of a single value of λ was inferred for all loci. as a result, we obtained two partitions, a partition with five clusters denoted as clusters a to e and a partition with four clusters. when a unique λ value was inferred for each locus, only the partition with four clusters was obtained. in this partition, cluster c was absorbed in clusters d and e. thus, we created three data sets, each including clusters c and d, clusters c and e, or clusters d and e, and reanalyzed them with dpart. λ was inferred for each locus. because clusters in each data set differed from each other, we concluded that the bull data set included five clusters. the dendrogram was generated on the basis of co-assignment probabilities calculated from mcmc samples. the data set was also analyzed using fum, fcm, and structure with k = as shown in figure although the result obtained with fum was consistent with that obtained with dpart, the results obtained with fcm and structure were not consistent with those obtained with fum and dpart. in the results of fcm and structure, the smallest cluster in the results of dpart and fum, i.e., cluster a, absorbed the members of clusters b, c, and e. in addition, the moderate cluster, cluster e, also absorbed the members of cluster d. on the other hand, the second smallest cluster in the results of dpart and fum, cluster b, absorbed very few members of the larger clusters in the results of fcm and structure. the pairwise fst values between the clusters detected by dpart are shown in table we consider the interpretation of these inconsistent results in the discussion section. discussion the bayesian method based on the dp prior can infer the number of populations and assign individuals, whereas the selection of the appropriate k value is often problematic when methods that run under a predefined k value are used. we examined the properties of this method to provide information that will be useful in practice. we showed that the sams sampler, which assigns multiple individuals simultaneously, was effective for the inference of population structures. because sams sampling was faster than mcmcmc techniques, the implementation of this sampling technique may be especially useful for large data sets. we also showed that a hyperparameter, named λ, which defines the prior distribution of allele frequencies, affected the performance of the method and its specification was important. this problem could be resolved by considering λ a variable. furthermore, we demonstrated that the dp prior method was suitable for data sets having unbalanced sample sizes among populations, whereas methods that implement structure-like algorithms were sensitive to unbalance. in particular, we found that the allele frequencies correlated model was the most sensitive. our results showed that both the sams sampler and mcmcmc were effective in improving the mixing of mcmc algorithms; however, the sams sampler was more effective than mcmcmc. we implemented the sams sampler at a frequency in which four iterations with the sams sampler were followed by one iteration with the gibbs sampler. although we examined other frequencies, such as two sams and five gibbs, one sams and one gibbs, and five sams and two gibbs, the results were almost the same. we selected this frequency simply to shorten the run time, because one attempt of the sams sampler is faster than one scan of the gibbs sampler. however, an accurate inference of the posterior distribution is hardly possible with the sams sampler alone and the gibbs sampler is necessary. the acceptance rate of the proposed states was extremely low, suggesting the difficulty of proposing new partitions. structurama implementing mcmcmc can probably increase its performance by increasing the number of chains or adjusting the temperature parameter to obtain the appropriate exchange rates among chains. however, calculation time increases as the number of chains increases and multiple runs may have to be performed to find appropriate values for the temperature parameter. therefore, we concluded that the sams sampler was more useful in practice and implemented this sampler in dpart to improve the mixing of mcmc algorithms. we showed that the choice of λ values affected the behavior of the dp prior method in the simulated and real data sets. since the preferable λ value varies depending on the uniformity of frequencies among alleles, it is desirable that the λ value at each locus is determined according to the allele frequencies. because hwler set λ to for any loci, its performance is probably inadequate for some data sets, although it implements the sams sampler. structurama does not state how it specifies λ, but it appears to fix λ for any loci at a certain value. inferring a unique λ value for each locus is a method of specifying the parameter for each locus. however, this approach resulted in less accurate inferences than assuming a single value for all loci when the levels of uniformity of frequencies among alleles were relatively similar among loci. we speculate that increasing the number of hyperparameters to be estimated may make inferences unstable. in contrast, assuming a single value for all loci was less accurate in data sets in which the levels of uniformity differed significantly among loci. therefore, it is difficult to state which assumption is more suitable for the data set of interest. in the chicken data set, the unique value assumption and the single value assumption provided the same partition. however, analysis in which a single value was assumed and inferred for all loci gave a slightly rougher partition. in the bull data set, the single value assumption gave a finer partition than the unique value assumption. in general, if the loci included in the data have not been evaluated well with regard to polymorphism, some loci may be much less polymorphic than others, and thus, the allele frequencies at these loci will be skewed. for such data sets, the unique value assumption will be suitable. on the other hand, if the loci included in the data are selected from a large number of candidate loci, they will be highly polymorphic, and thus, the allele frequencies at these loci are expected to be close to uniformity. therefore, the single value assumption will be suitable for such data sets. we speculate that the chicken data may be closer to the former case and, in such a case, even if an appropriate single value that leads to an accurate inference exists, the inference of such a value may be difficult. in contrast, the bull data set may be closer to the latter case because the loci included in the data set had been selected from a large number of candidates for parentage testing. therefore, the single value assumption may be preferable. our results showed that the behavior of the dp prior method depends on the selection of λ. we speculated that integration out of allele frequencies involves this dependency to some extent. thus, although we have not examined this speculation, the dependency may decrease by inferring allele frequencies. however, this will increase the calculation time and thus will not be suitable for large data sets. we found that unbalanced sample sizes among populations affect the behavior of dpart, fum, and fcm. fum and fcm were found to be sensitive to unbalanced sample sizes, and their performances were substantially affected, particularly by the presence of multiple minor subsets. the reason why dpart is suitable for unbalanced sample sizes is probably its prior assumptions about the assignment of individuals. as seen in eq., in the algorithm implementing the dp prior, clusters absorb individuals with higher probability as sample sizes increase, i.e., the "rich get richer" phenomenon occurs. thus, this algorithm is suitable for data sets with unbalanced sample sizes among populations. on the other hand, the algorithms implemented in fum, fcm, and structure assume that each population can contribute to the data set with equal probability. thus, these algorithms are suitable for data sets in which sample sizes are uniform among populations. when these methods are used for data sets with unbalanced sample sizes, they tend to cluster individuals such that the sizes become uniform among clusters, as observed in our experiments. the extent of sensitivity varied depending on the number of loci and the migration rates when the level of unbalance was fixed. thus, if genetic differences among populations are small and sample sizes are unbalanced, then the number of loci needed by these methods to correctly detect population structures is higher than that needed by dpart. although we compared these programs using only microsatellite data, the differences in behavior among them will not vary if snp data are analyzed because the differences will be due to differences among the prior assumptions and will not depend on the number of alleles. the sensitivity of fum and fcm to unbalance can probably be resolved by adding parameters that represent the mixing weights to the algorithms. as described by neal, when the prior distribution of the mixing weights is assumed to be a dirichlet distribution such as where k is the assumed number of populations and these weights are integrated out, we can obtain a predictive distribution for the assignment of the ith individual, which is similar to eq., but now in the parametric framework. fum, fcm, and other structure-like algorithms will be able to deal with unbalanced sample sizes by implementing this prior distribution. fcm, which implements the allele frequency correlated model, was shown to be more sensitive than the other programs to unbalanced sample sizes. a refutation of this may be that comparison among methods was performed for simulated data sets generated using the isolated population models that do not accord with the assumption of fcm. however, fcm was most sensitive to unbalance in the data sets that were generated by simulation method which was based on the assumption of fcm. we speculated that correlated allele frequencies are involved in the sensitivity of fcm. in the allele frequencies correlated model, the conditional posterior distribution for allele frequencies of population j is given as follows: where yjlh is the number of copies of allele h at locus l in individuals assigned to population j, palh is the ancestral frequency of allele h at locus l, jl is the number of alleles at locus l, and fj is where fj is the drift parameter for population j. thus, as the sample size for a population decreases, the effect of ancestral allele frequencies on the inference of allele frequencies for this population increases. on the other hand, ancestral allele frequencies are inferred from information on allele frequencies of populations and drift parameters. therefore, if unbalanced sample sizes are present, ancestral allele frequencies will be affected more strongly by the inferred allele frequencies for major subsets because those for minor subsets are substantially affected by ancestral allele frequencies themselves, and are thus, less informative for the inference of ancestral allele frequencies. consequently, the inferred allele frequencies for minor subsets will be affected by those for major subsets and will approach them. therefore, minor subsets may be prone to absorbing members of major subsets. we interpreted the inconsistency found in the bull data set on the basis of the knowledge obtained from simulations. if the results obtained with dpart and fum are correct, this inconsistency can be explained by the sensitivity to the unbalance of the allele frequencies correlated model. fcm and structure probably failed to detect the smallest cluster because of their sensitivity. the failure of cluster a to absorb the members of cluster d in the results of fcm and structure was due to the high level of differentiation between these clusters. although the unbalance between clusters d and e was moderate, fcm and structure also presumably failed to distinguish these clusters because of the relatively low pairwise fst between these clusters. when the data set that included only clusters d and e was analyzed by fcm, this program also failed to distinguish the two clusters. on the other hand, because cluster b was well differentiated from larger clusters, every program was able to detect the cluster. the presence of multiple minor clusters was also considered to reduce the performance of fcm and structure. in fact, when the two data sets that included only clusters a and b and clusters a and e were generated and analyzed by fcm, this program was able to distinguish the two clusters in each data set. although we admit that we have not proved that the results of dpart and fum are correct, we believe that our interpretation is appropriate because it can clearly explain the inconsistency. for the bull data set, runs were performed with fcm and structure. the results of the runs of each program were almost similar, and we considered them to be incorrect. however, in the simulated data sets with unbalanced sample sizes among populations, we occasionally observed that these programs or fum provided both correct and incorrect results when multiple runs were performed. thus, when the ad hoc statistic proposed by evanno et al. is used to select the true k value, this phenomenon possibly confuses the selection because it increases the standard deviation of likelihood at the true k value. the effect of sample sizes of populations on the performance of clustering programs was addressed in some studies. however, the effect of unbalanced sample sizes has been overlooked, and simulation studies for comparison of clustering methods usually assume uniform sample sizes among populations. because our results showed that sensitivity to unbalance in sizes varied among the methods, we recommend that comparative studies consider the effect of unbalance during analyses. through analyses of simulated data sets by dpart, we observed overestimation of k caused by small clusters that included only one or two individuals. this phenomenon increased the average k values and slightly affected the average. as discussed above, the dp prior method can efficiently detect minor subsets because of the "rich get richer" phenomenon. thus, we speculated that the overestimation was due to the fact that dpart detected individuals that were slightly distanced from the other members because even in simulated data sets, individuals harboring rare genotypes can be generated with low probabilities. this interpretation was supported by phylogenetic analysis based on genetic distances between individuals. in addition, this may be supported by the fact that overestimation became prominent as the number of loci, i.e., the power to detect population structures, increased. therefore, such small clusters are interpreted as overestimations in simulation studies, but will provide useful information in empirical studies because they indicate the presence of genetic discontinuity in the data sets.  this study characterized the bayesian method of implementing the dp prior and introduced a program, named dpart, in order to infer population structures more accurately than preceding programs based on the dp prior. first, we showed that the sams sampler, which is a technique for improving the mixing of mcmc algorithms, was effective for population structure analysis. implementation of this sampler was useful with regard to the accuracy of inference and computing time. second, we showed that a hyperparameter for the prior distribution of allele frequencies affected the behavior of the dp prior method. appropriate values can be specified by inferring this parameter. third, the dp prior method was shown to be suitable for analysis of data sets with unbalanced sample sizes among populations. in contrast, methods that implement structure-like algorithms were shown to be suitable for data sets with uniform sample sizes among populations, but not for data sets with unbalanced sample sizes. because these differences can yield inconsistent results among methods, we recommend using these methods concurrently. when the results obtained are inconsistent among methods, considering the effect of unbalanced sample sizes may be a key to interpreting the inconsistency. competing interests the authors declare that they have no competing interests. authors' contributions ao conceived the study, performed all statistical analyses, and drafted the manuscript. mn and mm collaborated on the algorithm design and the manuscript. all authors have read and approved the final manuscript. endnote a - gibbs sampler one scan of the gibbs sampler consists of the following steps. step remove the ith individual from η. step assign i to existing populations according to probabilities where ≤ j ≤ k, or a new population according to probability step update η. step repeat steps 1- for all i ∈ {...,n}. endnote b - sams sampler the algorithm of the sams sampler is as follows. step select two individuals, i and j, at random. step if i and j belong to the same population s, remove i and j from s and form two singletons, si = {i} and sj = {j}. if not, proceed to step step assign the individuals remaining in s to si or sj. the order of the assignment is randomly determined. the kth individual is assigned to si with probability otherwise, add the individual to sj. step the proposed partition η' is accepted with mh probability where q and q are the transition probabilities. using eq., the ratio of likelihoods in eq. can be written as follows: using eq., the ratio of prior probabilities in eq. can be written as follows: q is and q results from the products of probabilities of the assignment given in eq.. step if i and j belong to different populations, say si and sj, propose a new population s by merging si and sj. the proposed partition η' is accepted with the probability given in eq.. the ratio of likelihoods can be written as follows: and the ratio of prior probabilities is given as follows: q is and q is computed by splitting s into si and sj in a randomly determined order. endnote c - prior distribution of k prior distribution for k can be inferred using the following monte carlo procedure. step let the first individual belong to the first population and let k = step assign individual i = {,...,n} to existing or new populations with the probabilities noted in eq.. step record k after the assignment of the last individual. step repeat steps 1- for sufficient cycles to infer the distribution of k. endnote d - mean partition the algorithm for calculating the mean partition described by huelsenbeck and andolfatto is as follows. step pick a sampled partition as the initial state of the mean partition and calculate d, which is the sum of the partition distances between the mean partition and every sampled partition. step pick an individual i in the mean partition. propose new mean partitions by moving i to other populations in the mean partition and to a new partition. calculate the sum of partition distances between each proposed mean partition and each sampled partition, which is denoted as d'. step let d'min denote the minimum value of d's. if d'min <d, the corresponding proposed mean partition is accepted and d is replaced by d'min. step repeat step and for i = {,...,n}. step repeat steps and until d stops decreasing. supplementary material additional file 1 an r function for agglomerative hierarchical clustering. this performs agglomerative hierarchical clustering based on the results of dpart. click here for file additional file 2 executable file of dpart. this file will work on a windows platform. click here for file additional file 3 parameters. this file defines the parameters for dpart. click here for file additional file 4 input file names. this file defines the input file names for dpart. click here for file additional file 5 source code. this is the source code of dpart. click here for file additional file 6 header file. this is the header file of dpart. click here for file additional file 7 manual. this is the manual for dpart. click here for file additional file 8 bull data set. this file contains the genotypes of the bulls used in this study. this file is also an example of input data for dpart. click here for file acknowledgements we would like to thank two anonymous referees for their helpful comments. we also would like to thank the technical staff at maebashi institute of animal science for genotyping bulls.
third, we compared the dp prior method with other bayesian clustering methods and showed that the dp prior method was suitable for data sets with unbalanced sample sizes among populations. in contrast, although current popular algorithms for population structure analysis, such as those implemented in structure, were suitable for data sets with uniform sample sizes, inferences with these algorithms for unbalanced sample sizes tended to be less accurate than those with the dp prior method. the clustering method based on the dp prior was found to be useful because it can infer the number of populations and simultaneously assign individuals into populations, and it is suitable for data sets with unbalanced sample sizes among populations. a bayesian approach based on a dirichlet process prior is useful for inferring genetic population structures because it can infer the number of populations and the assignment of individuals simultaneously. here we presented a novel program, dpart, that implements the sams sampler and can consider the hyperparameter for the prior distribution of allele frequencies to be a variable. we showed that this sampler was effective for population structure analysis. second, we examined the effect of a hyperparameter for the prior distribution of allele frequencies and showed that the specification of this parameter was important and could be resolved by considering the parameter as a variable.
 a bayesian approach based on a dirichlet process prior is useful for inferring genetic population structures because it can infer the number of populations and the assignment of individuals simultaneously. however, the properties of the dp prior method are not well understood, and therefore, the use of this method is relatively uncommon. we characterized the dp prior method to increase its practical use.  first, we evaluated the usefulness of the sequentially-allocated merge-split sampler, which is a technique for improving the mixing of markov chain monte carlo algorithms. although this sampler has been implemented in a preceding program, hwler, its effectiveness has not been investigated. we showed that this sampler was effective for population structure analysis. implementation of this sampler was useful with regard to the accuracy of inference and computational time. second, we examined the effect of a hyperparameter for the prior distribution of allele frequencies and showed that the specification of this parameter was important and could be resolved by considering the parameter as a variable. third, we compared the dp prior method with other bayesian clustering methods and showed that the dp prior method was suitable for data sets with unbalanced sample sizes among populations. in contrast, although current popular algorithms for population structure analysis, such as those implemented in structure, were suitable for data sets with uniform sample sizes, inferences with these algorithms for unbalanced sample sizes tended to be less accurate than those with the dp prior method.  the clustering method based on the dp prior was found to be useful because it can infer the number of populations and simultaneously assign individuals into populations, and it is suitable for data sets with unbalanced sample sizes among populations. here we presented a novel program, dpart, that implements the sams sampler and can consider the hyperparameter for the prior distribution of allele frequencies to be a variable. 