A Faster and More Accurate Algorithm for Calculating Population Genetics Statistics Requiring Sums of Stirling Numbers of the First Kind
- *Infectious Diseases Programme, Department of Medicine, Yong Loo Lin School of Medicine, National University of Singapore, Singapore 119228, Singapore & Infectious Diseases Group, Genome Institute of Singapore, Singapore 138672
- †Centrum Wiskunde & Informatica (CWI), Science Park 123, 1098 XG Amsterdam, The Netherlands
- 2Corresponding author: Division of Infectious Diseases, Department of Medicine, Yong Loo Lin School of Medicine, National University of Singapore, Singapore 119228 & Laboratory of Bacterial Genomics, Genome Institute of Singapore, Singapore 138672. E-mail: slchen{at}gis.a-star.edu.sg
Abstract
Ewen’s sampling formula is a foundational theoretical result that connects probability and number theory with molecular genetics and molecular evolution; it was the analytical result required for testing the neutral theory of evolution, and has since been directly or indirectly utilized in a number of population genetics statistics. Ewen’s sampling formula, in turn, is deeply connected to Stirling numbers of the first kind. Here, we explore the cumulative distribution function of these Stirling numbers, which enables a single direct estimate of the sum, using representations in terms of the incomplete beta function. This estimator enables an improved method for calculating an asymptotic estimate for one useful statistic, Fu’s 
. By reducing the calculation from a sum of terms involving Stirling numbers to a single estimate, we simultaneously improve accuracy and dramatically increase speed.
- Population genetics statistics
- Evolutionary inference
- Stirling numbers of the first kind
- Asymptotic analysis
- Numerical algorithms
- Cumulative distribution functions
The dominant paradigm in population genetics is based on a comparison of observed data with parameters derived from a theoretical model (Casillas and Barbadilla 2017; Nielsen 2001). Specifically for DNA sequences, many techniques have been developed to test for extreme relationships between average sequence diversity (number of DNA differences between individuals) and the number of alleles (distinct DNA sequences in the population). In particular, such methods are widely used to predict selective pressures, where certain mutations confer increased or decreased survival to the next generation (Nielsen 2001). Such selective pressures are relevant for understanding and modeling practical problems such as influenza evolution over time (Grenfell et al. 2004) and during vaccine production (Chen et al. 2019); adaptations in human populations, which may impact disease risk (Wollstein and Stephan 2015; Quintana-Murci 2016); and the emergence of new infectious diseases and outbreaks (Wu et al. 2016).
Many population genetics tests are therefore formulated as unidimensional test statistics, where the pattern of DNA mutations in a sample of individuals is reduced to a single number (Nielsen 2001; Casillas and Barbadilla 2017; Fu 1997). Such statistics are heavily informed by combinatorial sampling and probability distribution theories, many of which are built upon the foundational Ewens’s sampling formula (Ewens 1972), which describes the expected distribution of the number of alleles in a sample of individuals, given the nucleotide diversity.
Ewens’s sampling formula not only was a seminal result for population genetics, but also established connections with combinatorial stochastic processes, algebra, and number theory (Crane 2016). For population genetics, in particular, Ewens’s sampling formula provided a key analytical result that finally enabled mathematical tests of the neutral theory of evolution (Crane 2016; Nielsen 2001). It has given rise to several classical population genetics tests for neutrality, including the Ewens-Watterson test, Slatkin’s exact test, Strobeck’s S, and Fu’s 
(Fu 1997; Strobeck 1987). Calculation of subsets of this distribution are useful for testing deviations of observed data from a null model; such subsets often require the calculation of Stirling numbers of the first kind (hereafter referred to simply as Stirling numbers). In particular, Fu’s 
has recently been shown to be potentially useful for detecting genetic loci under selection during population expansions (such as an infectious outbreak) both in theory and in practice (Wu et al. 2016). However, Stirling numbers rapidly grow large and thus explicit calculation can easily overwhelm the standard floating point range of modern computers.
In previous work, an asymptotic estimator for individual Stirling numbers was used to solve the problem of computing Fu’s 
for large datasets, which are now becoming common due to rapid progress in DNA sequencing technology (Chen 2019). Without such improved numerical methods, Fu’s 
calculations for data sets as small as 170 sequences can cause overflow, preventing the use of these statistics for genome-wide screens of selection. This algorithm based on estimating individual Stirling numbers solved problems of numerical overflow and underflow, maintained good accuracy, and substantially increased speed compared with other existing software packages (Chen 2019). However, there was still a need to estimate multiple Stirling numbers (up to half the total number of sequences). Here, we explore the potential for further increasing both accuracy and speed in calculating Fu’s 
by using a single estimator for the entire sum, which involves multiple Stirling numbers.
Methods
General definitions and theory
We take a population of n individuals, each of which carries a particular DNA sequence 
(referred to as the allele of individual i). We define a metric, 
to be the number of positions at which sequence 
differs from 
. Then, we denote the average pairwise nucleotide difference as 
(hereafter referred to simply as θ), defined as:
(1)We also define a set of unique alleles 
which have the property of 
. The ordinality of 
is denoted m, i.e., the number of distinct alleles in the data set.
Building upon on Ewens’s sampling formula (Fu 1997; Ewens 1972), it has been shown that the probability that, for given n and θ, at least m alleles would be found, is
(2)where 
is the Pochhammer symbol, defined by
(3)
is a Stirling number and is defined by:
(4)Fu’s 
is then defined as:
(5)Fu’s 
thus measures the probability of finding a more extreme (equal or higher) number of alleles than actually observed. It requires computing a sum of terms containing Stirling numbers, which rapidly become large and therefore impractical to calculate explicitly even with modern computers (Chen 2019).
Because of the relation in (4), the statistics quantity 
satisfies 
. Also, this relation and (3) show that 
are non-negative. We have the special values
(6)There is a recurrence relation
(7)which easily follows from (4). For a concise overview of properties, with a summary of the uniform approximations, see (Gil et al. 2007, §11.3).
We introduce a complementary relation
(8)leading to an alternate calculation for Fu’s 
of
(9)The recent algorithm considered in (Chen 2019) is based on asymptotic estimates of 
derived in (Temme 1993), which are valid for large values of n, with unrestricted values of 
. It avoids the use of the recursion relation given in (7).
In the present paper we derive an integral representation of 
and of the complementary function 
, for which we can use the same asymptotic approach as for the Stirling numbers without calculating the Stirling numbers themselves. From the integral representation we also obtain a representation in which the incomplete beta function occurs as the main approximant. In this way we have a convenient representation, which is available as well for many classical cumulative distribution functions. We show numerical tests based on a first-order asymptotic approximation, which includes the incomplete beta function. In a future paper we give more details on the complete asymptotic expansion of 
, and, in addition, we will consider an inversion problem for large n and m: to find θ either from the equation 
, when 
is given, or from the equation 
, when 
is given.
Remarks on computing 
When computing the quantity 
defined in (5), numerical instability may happen when 
is close to 1. In that case, the computation of 
suffers from cancellation of digits. For example, take 
, 
, 
. Then 
, and 
becomes about 6.6561 when using the first relation in (9). However, when we calculate 
and use the second relation, then we obtain the more reliable result 
.
We conclude that, when 
, it is better to switch and obtain 
from the sum in (8) and 
using the second relation in (9). A simple criterion to decide about this can be based on using the saddle point 
(see Remark 1 below).
A second point is numerical overflow when n is large, because 
rapidly becomes large when m is small with respect to n. For example, when 
, 
we have
(10)Therefore, it is convenient to scale the Stirling number in the form 
. In addition, the Pochhammer term 
in front of the sum in (2) will also be large with n; we have 
.
We can write the sum in (2) in the form
(11)Leading to a corresponding modification in the recurrence relation in (7) for the scaled Stirling numbers:
(12)To control overflow, we can consider the ratio
(13)This function satisfies 
if 
. For small values of n we can use recursion in the form
(14)For large values of n and all 
, we can use a representation based on asymptotic forms of the gamma function.
It should be observed that using the recursions in (7) and (12) is a rather tedious process when n is large. For example, when we use it to obtain 
for all 
, we need all previous 
with 
for all 
. A table look-up for 
in floating point form may be a solution. When n is large enough, the algorithm mentioned in (Chen 2019) evaluates each needed Stirling number by using the asymptotic approximation derived in (Temme 1993).
Data availability
Code implementing the new estimator for Fu’s 
in R is available at https://github.com/swainechen/hfufs.
Results And Discussion
Analytical results
The new algorithm is based on the following results, which we describe in two theorems.
Theorem 1. The statistics quantity 
introduced in (2) has the representation as an integral in the complex z-plane
(15)where n and m are positive integers, 
, θ is a real positive number, and 
is a circle at the origin with radius 
. The symbol 
denotes the Pochhammer symbol introduced in (3).
Observe that we have raised in 
the parameters n and m with unity; this is convenient in the mathematical analysis. The proof of this theorem will be given in the Appendix (Proof of Theorem 1).
Corollary 1. The complementary quantity 
introduced in (8) has the representation
(16)The main asymptotic result is given in the second theorem.
Theorem 2. 
has the representation
(17)where 
is the incomplete beta function defined by
(18)with
(19)The term 
is a function of which we give a one-term approximation in (32).
Corollary 2. The complementary quantity 
has the representation
(20)This follows from Theorem 1 and the complementary relation of the incomplete beta function
(21)Note also that the incomplete beta function in (17) has the representation (see (Paris, 2010, §8.17(i)))
(22)and from the complementary relation in (21) it follows that the function in (20) has the expansion
(23)The representation in this theorem in terms of the probability function 
shows the characteristic role of 
as a cumulative distribution function of the Stirling numbers. The representation can also be viewed as an asymptotic representation in which the incomplete beta function is the main approximant.
The proof of Theorem 2 can be found in the Appendix (Proof of Theorem 2), but we give here some preliminary information about functions used in the proof to explain the definition of the parameter τ in (17). It is a function of θ and arises in certain transformations of the integral given in Theorem 1. For this we need the function
(24)and its derivative
(25)With the function 
we can write (15) in the form
(26)Then the saddle point of the integral in (26) follows from the equation
(27)There is a positive saddle point 
when 
.
Next to these functions we introduce a function for complex values of a variable t:
(28)where 
. These functions are related by
(29)with condition 
. In this way, using this relation as a transformation of the variable z to t, we can write (26) as
(30)The parameter τ in Theorem 2 is defined as the positive solution of the equation
(31)In Figure 1 we show the graphs of 
(Figure 1A) and 
(Figure 1B) for 
, 
. For these values the saddle points are 
and 
. The sign condition 
for the relation in (29) means that the left branches of the convex curves correspond with functions values for 
and 
, and the right branches with values for 
and 
. Clearly, we have a one-to-one relation between the positive z and t-variables.
Graphs of 
(A) and 
(B) for 
, 
, with 
and 
.
A first-order approximation of the function 
in (17) and (20) reads
(32)where
(33)and the function 
is defined in (30). The value 
follows from evaluating 
(see (30)) at 
, by observing that both functions 
and 
vanish when 
(hence, 
). Then, l’Hôpital’s rule can be used to obtain 
.
In Figure 2 we show the error curves 
in (34) for Fu’s 
(9) for 
. We show examples for 
, 
(Figure 2A) and 
, 
(Figure 2B). The solid curves are for 
when using 
, the dashed curves when using 
with the asymptotic estimate given in (32). For ease of visualization, the error 
has been multiplied by a factor 10 or 100 in Figure 2. We have used the following mollified error function
(34)where 
is the approximation of 
. This mollified error is exactly the relative error unless 
is small. Because 
will vanish when 
(which also means that θ is near the transition value 
(in Figure 2A) and 
(in Figure 2B) (see Remark 1)), we cannot use relative error for all 
. This explains the non-smooth curves in Figure 2.
Mollified error in estimating Fu’s 
for 
, 
and 
(A) and for 
and 
(B). The data for the dashed curves are multiplied by a factor of 10 (A) and 100 (B), to make the error curves visible in the figures. Refer to the text for further details.
The final estimator is based on the representations in (17) and (20) and the first order approximation in (32), which are used to calculate Fu’s 
with one of the two relations in (9) depending on whether 
, decided as described above.
Implementation and numerical results
We first summarize the steps to compute Fu’s 
by using (9) and the first-order approximations (see (32) and (17) or (20))
(35)or
(36)for large n, 
and 
.
1 Compute the saddle point
, the positive zero of ; see (27).2 With
, the positive zero of (see (28)), compute τ, the solution of the equation (see (31))(37)with defined in (24) and defined in (28). When there is one solution . When there are two positive solutions, and we take the one that satisfies the condition .3 When
, hence , compute the approximation of by using (35), and from the first relation in (9).4 When
, hence, , compute the approximation of by using (36), and from the second relation in (9).
In Table 1, we show the relative errors in the computation of 
defined in (5). The values of n, m, and θ correspond with those in Table 1 of (Chen 2019). The asymptotic result is from (35). Computations were done with Maple, with Digits = 16. The “exact” values were obtained by using Maple’s code for 
, which computes the Stirling numbers of the first kind.
defined in (5) using the asymptotic estimator in (35)We additionally performed a comparison with the recently published algorithm in (Chen 2019). We performed 10,000 calculations with each algorithm and compared the results with an exact calculator. As expected, since the previous algorithm required estimating a Stirling number for each term of the sum, while the current asymptotic estimate directly calculates the sum, both error and compute speed were improved. Relative error for the single term estimate in (35) was well controlled at 
for nearly 99% of the calculations; for 411 calculations where the previous hybrid estimator had an error 
, the estimate in (35) was more accurate in all but one case (
; 3.08e-3 relative accuracy using (Chen 2019); 3.32e-3 relative accuracy using (35)) (Figure 3). Further analysis of the relative error demonstrated that it peaks at intermediate values of 
, depending on θ. These correspond to parameter choices near the transition values 
, where t approaches 
and z approaches 
in the calculation; notably, they remain well controlled (all values 
mollified error) regardless of θ. The asymptotic behavior (lower relative error) can also be seen as both n and m increase in Figure 3B.
(A) Comparison of relative error of the estimator from (Chen 2019) and the single term asymptotic estimator in (35). Relative error for each is calculated against the arbitrary precision implementation described in (Chen 2019). In total, 10,000 calculations were performed with n randomly sampled from a uniform distribution between 50 and 500; m between 2 and n; and θ between 1 and 50. A solid diagonal line is drawn at 
. Dotted lines are drawn at a relative error of 0.001. Numbers within each quadrant defined by the dotted lines indicate the number of points in each quadrant. The red dot indicates the one case where the relative error was 
and the error of (35) was greater than the estimator from (Chen 2019). (B) Comparison of mollified error (34) as a function of m. For this plot, we fixed 
(solid lines) or 500 (dotted lines) and 
(as indicated by different line colors).
The fewer calculations led to a clear improvement in calculation speed (median 54.6x faster; Figure 4). The speedup also depends on the parameter choices; in general, the speed advantage is greater when the hybrid calculator requires many calculations (namely, when m is small relative to n, as the hybrid calculator performs the sum in (2)) (Figure 4B).
(A) Comparison of run times between the hybrid algorithm from (Chen 2019) and the single term asymptotic estimator in (35). 100 iterations were run, each with 10,000 calculations; the time elapsed for each set of 10,000 calculations was recorded and plotted here. The same set of parameters were used for each algorithm. The order of running the algorithms was alternated with each iteration. The dark horizontal line indicates the median, the box indicates the first and third quartiles, the whiskers are drawn at 1.5x the interquartile range, and outliers are represented by open circles. The median for the hybrid algorithm is 62.64 s; the median for the asymptotic algorithm is 1.17 s. (B) Detailed benchmarking for 
(open violins) or 500 (gray violins), 
, and 
. Fold speedup (ratio of the time taken for the hybrid calculator to that taken for the aysmptotic estimator) is plotted on the y-axis. Each dot represents one set of parameters; the violin plots summarize the density of points on the y-axis. Times were calculated for 100 iterations of each estimator for the same parameter values.
Conclusion
The rapid growth of sequencing data has been an enormous boon to population genetics and the study of evolution. Traditional population genetics statistics are still in common use today. The statistics Fu’s 
and Strobeck’s S have been difficult to calculate on modern, large data sets using previous methods; we now further improve both accuracy and speed for the calculation of Fu’s 
for such data sets, using the main estimator in (35). Our plan for a paper about the ability to invert the calculation provides additional future directions in understanding the performance of these statistics. Therefore, the methods used herein may be useful for the development of new statistics that more effectively capture different types of selection.
Acknowledgments
The authors thank the reviewers for their helpful comments and suggestions. SLC acknowledges Shyam Prabhakar and members of the Chen lab for fruitful discussions. NMT acknowledges CWI, Amsterdam, for scientific support; he thanks Edgardo Cheb–Terrab (Maplesoft) for his help with the Maple code. SLC was supported by the National Medical Research Council, Ministry of Health, Singapore (grant numbers NMRC/OFIRG/0009/2016 and NMRC/CIRG/1467/2017). NMT was supported by the Ministerio de Ciencia e Innovación, Spain, projects MTM2015-67142-P (MINECO/FEDER, UE) and PGC2018-098279-B-I00 (MCIU/AEI/FEDER, UE).
APPENDIX
Proof of Theorem 1
We use the integral representation of the Stirling numbers that follows from the definition given in (4). That is, by using Cauchy’s formula,
(38)where 
is a circle around the origin with radius R. We can take R as large as we like. As in (Temme, 1993, §3), it is convenient to proceed with
(39)Using the definition of 
in (2) we have
(40)and using (39) we obtain
(41)We can take 
to have 
on the circle 
, and we can perform the summation to ∞, because all terms with 
do not give contributions. In this way we obtain the requested integral representation
(42)This concludes the proof of Theorem 1.
Corollary 1 now follows by using the theory of integrals of analytical functions on complex contours. We have assumed that 
, but we can take 
while picking up the residue at 
. The result is
(43)This gives the relation in Corollary 1.
Remark 1. When θ crosses the value 
, 
becomes (almost) 
. Especially when the parameters m and n are large, 
starts with very small values for small θ, becomes close to 
when 
, and quickly becomes 1 as θ increases. We call 
the transition value for θ.
For fixed values of n and θ, there is also a transition value for m; refer to this transition value as 
. When n is large, 
starts at values near 1 for small m, it becomes about 
when m nears 
, and it becomes quickly small as 
.
Proof of Theorem 2
The relation in (29) between the functions 
(see (24)) and 
(see (28)) can be used as a transformation of the variable z to t, as in (Temme, 1993, §3). The result is the integral representation in (30). In Figure 1 we have shown the relationship between z and t.
The function 
in (30) has a pole in the t-domain; refer to this pole as 
. This then corresponds with the pole at 
in the z-domain. The relation between τ and θ follows from the transformation given in (29). In other words, τ is defined by the equation
(44)where the sign-convention follows from the one used for (29). We can express the existence of the pole of the function 
defined in (30) by writing
(45)In asymptotic analysis, the presence of such a pole is of great interest, especial when (in the t-domain) the saddle point (here 
) is close to a pole (here τ), or even when these points coalesce. See, for example, (Temme, 2015, Chapter 21). Usually, the error function is introduced to handle the asymptotic analysis; in the present case, we use an incomplete beta function. We split off the pole from 
and write
(46)where we assume that 
is well defined at 
. To find A we use the analytical relation in (29) between t and z, in particular at 
(or 
). Applying l’Hôpital’s rule, we find that 
as 
, which gives 
1. Hence, substituting this form of 
in (30), we find
(47)The radius of the circle 
in the first integral is larger than τ. For the second integral, we take a circle C around the origin such that the singularities of 
are outside the circle.
In Proof of the incomplete beta relation below, we prove that the first integral in (47) can be evaluated in terms of the incomplete beta function as shown in Theorem 2. We can then write (47) as
(48)where
(49)A first-order approximation of (49) follows from replacing 
by its value at the saddle point 
. This gives
(50)where
(51)This expression of 
follows from (33). In a future publication we will give details about the complete asymptotic expansion of the term 
.
Proof of the incomplete beta relation
We give a proof of the claim that the incomplete beta function in (48) equals the first integral in (47). That is,
(52)where 
is a circle at the origin with radius larger than τ. We have, using the definition of 
in (28),
(53)which is the relation in (22). In the second line we have used a finite number of terms of the infinite expansion of 
because terms with 
do not give a contribution.
Footnotes
Communicating editor: A. Kern
- Received March 9, 2020.
- Accepted August 23, 2020.
- Copyright © 2020 Chen, Temme
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.