Kernel density estimation of allele frequency including undetected alleles

KDE of allele frequency

Allele frequency ranges from 0 to 1, and the sum of all frequencies is always 1. Here, we consider the nucleotide sequence space as discrete; the distance between alleles is equivalent to the number of substitutions between them, and we do not consider mixed bases. We assume that all the alleles has the same length, and the nucleotide sequence space has dimensions equal to the length of the alleles. In this case, the allele frequency can be considered as a probability mass function defined on the nucleotide sequence space. Although KDE is usually applied to a continuous probability density function, the allele frequency lies on the discrete nucleotide sequence space. Furthermore, whereas ordinal real-number components of space change bidirectionally (positive and negative), a base as a component of the sequence space can change (mutate) in three directions (bases). Due to this discreteness and particularity of the sequence space, KDE of allele frequency is abnormal but possible.

Choice of the kernel function in KDE does not have a large impact on the result of KDE, but because the KDE of allele frequency is on a unique sequence space as described above, a special kernel function is required. To define the kernel function, we first assume the probability $P\left(g,\mu \right)$ that one base will be another after $g$ generations (e.g., A = >T). This is calculated from the probability recurrence relation as:

$$P\left(g,\mu \right)=1/4+\left(1/3\mu -1/4\right){\left(1-4/3\mu \right)}^{g-1},$$where $\mu$ is the mutation rate per generation per base locus, which is concordant with the Jukes-Cantor substitution model (Jukes & Cantor, 1969). This equation assumes no recombination and no population size changes, but the content of the equation is not important. Then, we consider the probability $Q_P\left(a,x,P\right)$ that one sequence $a$ of length $l$ will be another sequence with $x$ new bases after $g$ generations (e.g., AAA = >ATT when $l=3,x=2$). Assuming the independence of the mutations, it will be:

$$Q_P\left(a,x,P\right)={\left\{1-3P\left(g,\mu \right)\right\}}^{l-x}P{\left(g,\mu \right)}^x.$$

Now we can consider that $P\left(g,\mu \right)$ and $Q_P\left(a,x,P\right)$ correspond to the bandwidth (smoothing parameter) and kernel function, respectively. This is because summing up $Q_P\left(a,x,P\right)$ for all $x$ equals to 1, and $P\left(g,\mu \right)$ can be considered as one variable which defines the shape of the kernel function $Q_P\left(a,x,P\right)$. The mutation rate $\mu$ and the number of generations $g$ in $P\left(g,\mu \right)$ are not necessary to conduct KDE because the value of $P\left(g,\mu \right)$ is selected as one variable through a bandwidth optimization method (described in later section). In this scheme, we assume the sequence length $l$ is fixed. Insertions and deletions must be excluded since they are completely different from substitutions, especially in the point that multiple cites are often inserted or deleted by a single mutation. Using this kernel function, KDE of allele frequency can be expressed by the following equation:

(1) $$\begin{array}{c}\hat{f}\left(a\right)={\displaystyle \frac{1}{n}{\sum }_{i=1}^n{\sum }_{x=0}^l{Q}_P\left(a_i,x,P\right),\mathrm{\#}\left(1\right)}\end{array}$$where $\hat{f}(\cdot )$ is the estimated frequency mass function of the allele frequency, $a$ is all alleles, $a_i$ is an allele in the sample, $n$ is the sample size and $x$ indicates the number of base changes from $a_i$. Please note that the probabilities given by $Q_P\left(a_i,x,P\right)$ are recorded on coordinates or alleles on the nucleotide sequence space. We will discuss the space later.

There are many kinds of methods to optimize the bandwidth (Heidenreich, Schindler & Sperlich, 2013), but we could only implement the least square cross validation (LSCV) method and the likelihood cross validation (LCV) method (Zambom & Dias, 2013). This is because many other methods, such as the plug-in methods, depend on the Taylor expansion for approximation (Zambom & Dias, 2013), but differential of the nucleotide sequence space is difficult to define, and only the first order differential can be considered for the nucleotide sequence space since the maximum distance between bases is one.

Nucleotide sequence space and its compression

KDE calculation is conducted on the nucleotide sequence space. However, the sequence space is generally too large for ordinary computers. For example, the space of 1,000-bp sequences has $4^{1000}\approx {10}^{602}$ coordinates. Only when sequences have a few bases, KDE can be conducted directly on the sequence space. Otherwise, the sequence space must be compressed.

One way to compress the nucleotide sequence space is to use “distance space”. The distance space assigns one coordinate to each detected allele. Undetected alleles which are not obtained by sampling but estimated by KDE are recorded by their distances from the detected alleles. For example, when the detected alleles are “AA” and “AT”, the total distance space will be (0, 1), (1, 0), (1, 1), (1, 2), (2, 1) and (2, 2), where the first and second components of the coordinates correspond to the distance from the sequence “AA” and “AT”, respectively. Figs. 1 and 2 show the results of KDE in the sequence and distance spaces, respectively, when three “AA” and three “AT” alleles are detected. The size of the sequence space for the detected alleles “AA” and “AT” is $4^2$ whereas the corresponding distance space has only six coordinates as shown above. For example, the coordinate (1, 1) in the distance space corresponds to “AG” and “AC” in the nucleotide sequence space. Thus, the distance space can reduce the space size.

However, further compression of the sequence space is practically necessary. For this purpose, maximum range of the kernel function $Q_P\left(a,x,P\right)$ is limited to a given positive integer by changing $l$ in the Eq. (1) to another integer. We refer to this integer as the “mutation number”. Owing to this range limitation, the probability density of undetected alleles farther than the mutation number from a detected allele becomes zero, and such undetected alleles no longer need to be recorded on the distance space. Therefore, all components of the coordinates in the distance space larger than the mutation number can be recorded as infinity and serve for space compression. For example, when the detected alleles are “AAA”, “AAT” and “GGG”, the size of the distance space is 29. If mutation number 1 is applied to this distance space, the size of the space will be seven composed of (0, 1, 3), (1, 0, 3), (3, 3, 0), (1, 1, Inf), (1, Inf, Inf), (Inf, 1, Inf) and (Inf, Inf, 1). Setting a larger mutation number requires more memory and computation time but provides more precise results. However, this improvement in precision decreases as the mutation number increases because the estimated frequency of the undetected alleles decreases as the undetected alleles are farther from a detected allele. On the other hand, setting the mutation number 0 is equivalent to not applying KDE and calculating the allele frequency by counting (cf. gray bars in Figs. 1 and 2). Practically, we suppose that using mutation number of 1 or 2 is sufficient in most cases, but larger mutation number may be desirable when using data with high mutation rates, such as virus sequences.

Approximate nucleotide diversity

Calculation of nucleotide diversity from the result of KDE is not possible when using distance space. Therefore, we devised approximate nucleotide diversity as follows. Nucleotide diversity $\pi$ is sum of the product of all pairs of allele frequencies and substitution rate between the pair, namely,

$$\pi ={\sum }_{i=1}^N{\sum }_{j=1}^N{p}_i{p}_j{\pi }_{ij},$$where N is the total number of allele types, p_i and p_j are allele frequency of the sequence i and j, and π_ij is the substitution rate (Nei & Li, 1979; Nei & Tajima, 1981). The substitution rate is the distance between alleles divided by the sequence length. The frequencies and distances between the detected alleles can be directly obtained from the frequencies estimated by KDE and the coordinates of the distance space, respectively. The distance between the undetected alleles cannot be obtained. The frequencies of the undetected alleles are generally much lower than those of the detected alleles. Therefore, we ignore the terms for two undetected alleles in the calculation of the approximate nucleotide diversity. For the terms between a detected allele and an undetected allele, frequencies and distance are obtained by extra calculation even when the distance is recorded as “Inf” in the coordinate of the undetected allele. Thus, approximate nucleotide diversity is calculated by ignoring the terms between undetected alleles. Additional calculations to obtain the frequency and distance between detected and undetected alleles are shown in the Supplemental Article.

Implementation

The KDE of allele frequency was implemented using R (version 4.3.2; R Core Team, 2023). This implementation used ape package (version 5.7.1; Paradis & Schliep, 2019) to import Fasta files and depended on matrixStats package (version 1.0.0; https://cran.rstudio.com/web/packages/matrixStats/index.html) to accurately calculate the logarithm of the sum of the exponentials of allele frequencies. Although we prepared the parallelization of the calculation using snow package (version 0.4-4; Tierney, Rossini & Li, 2008), the KDE calculation requires large memory usage for creating space, and the time to copy the memory for parallelization overwhelmed the merits of acceleration by parallelization. Thus, parallelization through multi processes using snow package was not useful.

Simulation experiment

The population nucleotide diversities under mutation rates 0.1 and 0.01 were 0.1611368 and 0.0175182, respectively. The squared accuracy of nucleotide diversity showed generally the same tendency as the accuracy. Therefore, we hereafter explain the accuracy of the nucleotide diversity. All the raw data, their graphs and their comparing graphs are shown in Supporting Information. Because all the data are too abundant to show concisely, we representatively show four graphs of the concordance rate and accuracy (Figs. 3–6).

Comparison of mutation number

A mutation number of zero showed the best accuracy for nucleotide diversity (no-KDE in Figs. 3 and 4). Mutation number 2 was the worst in most cases, but mutation number 1 was the worst under a mutation rate of 0.1 and bandwidth optimization by LCV. The bias of nucleotide diversity showed the same trend as nucleotide diversity accuracy. Whereas the bias under mutation number 0 was around 0, those under mutation numbers 1 and 2 were mostly negative.

The concordance rate of allele frequency between the sample and population increased mostly proportional to the sample size (Figs. 5 and 6). Mutation number zero was the highest in the almost all cases. Exceptionally, the parameter set of bandwidth optimization under LSCV and LCV, sample size 20, mutation rate 0.1 and mutation number 1 showed a better concordance rate than that of mutation number 0.

Comparison of bandwidth optimization method

The accuracy of nucleotide diversity under LCV was almost always better than that under LSCV (Figs. 3 and 4). Only when the sample size 20, mutation rate 0.1 and mutation number 1, LSCV showed slightly better accuracy of nucleotide diversity than LCV. The bias of nucleotide diversity was mostly negative, and the bias under LCV was always closer to 0 than that under LSCV.

LCV showed higher concordance rate than LSCV in most cases (Figs. 5 and 6). LSCV showed a higher concordance rate than LCV in two cases: The case of the sample size 20, mutation rate 0.1 and mutation number 1, and the case of the sample size 20, mutation rate 0.1 and mutation number 2.

Comparison of mutation rate

The accuracy of nucleotide diversity under mutation rate 0.01 was always better than that under mutation rate 0.1 when using LSCV. When using LCV, the accuracy of nucleotide diversity under mutation rate 0.1 was better than that under mutation rate 0.01 in the following cases: The case of the sample size 60 under mutation number 1 and case of the sample size 40 and 60 under mutation number 2. The bias of nucleotide diversity was mostly negative. Whereas the bias of mutation rate 0.01 under LSCV was clearly closer to zero, that of under LCV did not show such large difference between the mutation rate 0.1 and 0.01.

Concordance rate under mutation rate 0.01 was almost always higher than that under mutation rate 0.1. Only in the case of the sample size 80 using LCV, the concordance rate under mutation rate 0.1 was higher than that under mutation rate 0.01.

Real data experiment

The population nucleotide diversity was 0.019857. Again, the squared accuracy of nucleotide diversity showed generally the same tendency as the accuracy. Therefore, we hereafter explain the accuracy of nucleotide diversity. All the raw data, their graphs and their comparing graphs are shown in Supporting Information. We representatively show two graphs of the concordance rate and accuracy (Figs. 7 and 8).

Except when the sampling ratio was 20%, LCV selected almost zero bandwidths, resulting in the almost same values of statistics between LCV and noKDE. Accuracy of the nucleotide diversity under LCV was almost always better than that under LSCV (Fig. 7). LSCV showed better nucleotide diversity accuracy than LCV only when the sampling ratio was 20%. However, even when the sampling ratio was 20%, the value with noKDE was better than that with LSCV. The bias of nucleotide diversity under LSCV was always negative. The bias under LCV was always closer to zero than that under LSCV. LCV and noKDE showed higher concordance rates than LSCV in most cases, but LSCV was higher when the sampling ratio was 20% (Fig. 8).