Genotype Value Decomposition: Simple Methods for the Computation of Kernel Statistics

Recent advances in sequencing technologies enable genome‐wide analyses for thousands of individuals. The sequential kernel association test (SKAT) is a widely used method to test for associations between a phenotype and a set of rare variants. As the sample size of human genetics studies increases, the computational time required to calculate a kernel is becoming more and more problematic. In this study, a new method to obtain kernel statistics without calculating a kernel matrix is proposed. A simple method for the computation of two kernel statistics, namely, a kernel statistic based on a genetic relationship matrix (GRM) and one based on an identity by state (IBS) matrix, are proposed. By using this method, calculation of the kernel statistics can be conducted using vector calculation without matrix calculation. The proposed method enables one to conduct SKAT for large samples of human genetics.


Introduction
A very large number of human genome sequences are now available for the study of human genetics, because of recent advances in genome sequencers. A recent study has shown that rare variants substantially contribute to phenotype variation. [1] Because each linkage disequilibrium block can be analyzed independently, increases in the number of sites can be tackled with parallel computation. [2] However, the statistical power of classical single-marker association analysis for rare variants is quite limited.
To address this challenge, rare and low-frequency variants are often grouped into gene or pathway levels, and the effects of multiple variants evaluated are based on collapsing methods. [3,4] The sequential kernel association test (SKAT) [5,6] is one such popular method. SKAT applies a test statistic S, which is defined by the quadratic form, S = y T Ky, where y is column vector of the phenotype defined by Equation (1).

DOI: 10.1002/ggn2.202100066
where y(i) is the phenotype value of the i-th individual and n is the sample size. In the following section, we assume the average of the elements of y is 0.
Evaluating the probability density of the null distribution of S is important for conducting SKAT, but it requires computing a matrix related to the genotype covariance between markers, which requires a very long computational time. When the length of y is n and the size of the matrix K is n 2 . For example, when the number of people is 10 000, the size of the kernel matrix is 100 000 000. A genetic relationship matrix (GRM) among individuals is used in genome-wide complex trait analysis [7] and in principal component analysis. [8] Identity by state (IBS) defines similarity between individuals as the number of shared alleles. The IBS kernel is used in linear regression [9,10] and SKAT. [5] The aim of the present study is to develop simple methods for the computation of these two kernel statistics without calculating a GRM and an IBS matrix explicitly.

Genotype Value Vectors
In the present study, all sites are assumed to be biallelic, namely, each site has a reference allele and an alternative allele. Let us define genotype value vectors at site k. Let a k (i) be 1 when the individual i is a homozygote of a reference allele at site k, otherwise a k (i) = 0. Let b k (i) be 1 when the individual i is the heterozygote at site k, otherwise b k (i) = 0. Let c k (i) be 1 when the individual i is a homozygote of an alternative allele at site k, otherwise c k (i) = 0. The vectors a k , b k , c k are defined by Equation (2).
Let us denote a k , b k , and c k as the genotype value vectors. Because a k (i) + b k (i) + c k (i) = 1, it is worth noting that where 1 is defined by 1 = ( 1 1 ⋯ 1 ) T The alternative allele frequency, p, is obtained by p = (b k T 1 + 2c k T 1)∕(2n). In the follow- Genotype ing section, the Hardy-Weinberg equilibrium is assumed for this site. Namely, the frequencies of heterozygotes and homozygotes of the alternative alleles are 2p(1 − p) and p 2 , respectively.

The GRM Kernel
The allele values are 0 for the reference allele, and 1 for the alternative allele. The separator between the alleles is "/" as used in the variant call format. [11] Let g k (i) be the genotype value for individual i at site k. In the present study, g k (i) is the number of alternative alleles. The relationship between the genotype of the individual i at site k and g k (i) is shown in Table 1. The vector g k is defined by Equation (4).
Let us denote the GRM at site k as X k . We subtract the mean }∕n to obtain a matrix with row sums equal to 0. The ij-th element of X k at site k is obtained using Equation (5).
Subsequently, the matrix X k can be obtained using the genotype value vectors.
Let us define a new matrix G k . As shown in Table 1, the ij-th element of G k at site k can be calculated using Equation (7).
Subsequently, the matrix G k would be obtained using the genotype value vectors.
It is worth noting that y T X k y = y T G k y + y T 1g k T y + y T g k 1 T y + 2 y T 11 T y Because y T 1 = 1 T y = 0, we obtain By using the distributivity and associativity of matrix production, we obtain Q k is a scalar value of site k defined by Equation (13): because the transpose of a product of matrices is the product, in the reverse order, of the transposes of the factors. Note y T b k + 2y T c k is a scalar that can be obtained as

The IBS Kernel
IBS defines similarity between individuals as the number of shared alleles. The IBS kernel is used in linear regression [9,10] www.advancedsciencenews.com www.advgenet.com Table 2. Relationship between genotype values and identities by state (IBS). and the SKAT. [5] Let IBS k (i) be the ij-th element of the IBS matrix, IBS k , at site k, which denotes the number of shared alleles by subjects i and j at site k.
Thus, the IBS matrix at site k is obtained by R k is a scalar value of site k defined by Equation (17).
By using the distributivity and associativity of matrix production, we obtain y T a k and y T c k are scalars that can be obtained using the inner product of two vectors. By using Equation (3), we can obtain When multiple single-nucleotide polymorphisms (SNPs) are investigated, the entire GRM and IBS matrices are obtained using S = ∑ l k=1 w k Q k and S = ∑ l k=1 w k R k , respectively, where w k is weight of site k and l is the number of sites. SKAT allows the incorporation of flexible weight functions. [12] Weights can normalize each data column to have the same variance [13] and can increase the power of tests. [13]

Computer Simulations
To evaluate the new method, I performed computer simulations. The python scripts used in the computer simulation are shown in Material 1, Supporting Information. The usages are in Material 2, Supporting Information. This program is ready for data analysis.

Genotype Selection
SNPs on the SLC22A2 gene that are known to affect uric acid levels [1,[14][15][16] were selected. Then, genetic variation of these SNPs of 2504 individuals were downloaded from the 1000 Genomes Project. Monomorphic sites were excluded. As a result, the sites in Table 3 were used in the computer simulation.

Phenotype Generations
The heterozygous individuals and the homozygous individuals of alternative allele of the uric acid level were set to be 1.0 µg dL −1 lower than the homozygotes of the reference allele. A random variable that follows the normal distribution with mean 0.0 µg dL −1 and standard deviation 1.0 µg dL −1 was added to the uric acid level of each individual in the simulation as an environmental factor of uric acids level. These values are similar to the observed values. [1]

Calculation of Test Statistics and Permutation Tests
For each of these phenotypes, the test statistics of GRM and IBS were observed. Then the permutation tests were performed with 1 000 000 permutations to calculate the probability of exceeding the observed score. The significance level was set to be 5 × 10 -6 , because the number of tests of genome-wide SKAT will be 10 4 . Each permutation test was repeated ten times (n = 10).  Table 4 shows that there is no significant difference between the GRM and IBS in the statistical power (the chi-square test, n = 10, P > 5%). Table 4 also shows that permutation tests can be conducted in a short period of time by using the methods proposed in this study.

Discussion
We demonstrate that necessary variant/phenotype association test statistics can be obtained without obtaining eigenvalues and eigenvectors of GRM and IBS matrices, in the present study. The method is referred to as genotype value decomposition. The new methods proposed in this study are conducted with computational time of O(n), where n is the sample size. Notably, these new methods are applicable for common variants as well as rare variants, even though the methods were developed for the association tests for rare variants. Sparse matrix computation can be used when all of variants are rare. When the alternative allele frequency is very small, homozygotes of the alternative allele are very rare, so that c k is ignorable. In other words, Q k can be approximately obtained byQ k ≈ (y T b k ) 2 . Under the same condition, (y T a k ) 2 ≈ (y T b k ) 2 and (y T c k ) 2 ≈ 0, so that R k is approximately equal to 2Q k .
On one hand, when all sites are independent, the necessary probability density functions can be calculated using convolution of the probability density functions of all sites. On the other hand, it is difficult to obtain convolution of the probability density functions when the sites are linked and dependent on each other. In such cases, a permutation test is used. [6] Because the statistics calculated by the new method are not approximations but exact values, the null distributions of these statistics are exactly the same as the test statistics with calculating GRM and IBS matrices. Wu et al. [5] showed that the test statistics approximately follow the chi-square distribution. Furthermore, because the distribution is derived from an asymptotic distribution of its statistics, the p-values for datasets with an insufficient number of samples may be inaccurate, which could cause inflation or power loss. [] In a permutation test, the test statistic null distribution can be approximated by fully resampling the observed traits without replacement. The proposed method can be useful for reducing computational time to obtain p-values using resampling methods.

Conclusion
In the present paper, a genotype value decomposition method is proposed for handling the kernel matrices. The method can be re-ferred to as genotype value decomposition. By using this method, calculation of the null distribution of the kernel statistics can be conducted with time complexity O(n). The proposed method enables one to conduct SKAT for large samples of human genetics.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.