For a complex trait, heritability () gives the genetic determination of its variation. Given the emergence of biobank-scale data, a more powerful method is needed to estimate . Based on the framework of Haseman-Elston regression (RHE-reg), we integrate a fast randomization algorithm to estimate , and RHE-reg can tackle biobank-scale data, such as UK Biobank (UKB), very efficiently. Furthermore, we present an analytical solution that balances computational cost and precision of the estimation, a property that is important in dealing with biobank-scale data. We investigated the performance of the RHE-reg in simulated data and also applied it for 81 UKB quantitative traits; as tested in UKB data of nearly 300,000 unrelated individuals, it took on average about 4.5 hours to complete an estimation when used 10 CPUs. We extended the application of RHE-reg into distributed datasets when privacy is not compromised. As shown in UKB and simulated data the performance of RHE-reg was accurate in estimating . The software for estimating SNP-heritability for biobank-scale data is released.

Method description

Wu and Sankararaman proposed a randomized implementation for the Haseman-Elston regression (RHE-reg), which dramatically reduced the computational time from to in dealing with [5]; is the genetic relationship matrix for individuals on markers, and see its detailed definition in the section below. It is clear that a large , indicating more iteration of the presented algorithm, is helpful in improving precision, but it is unsolved how to get an estimate for and its role in determining the boundaries of key statistics, upon the standard errors of the randomized estimator [8]. This work is in general consistent with Wu and Sankararaman’s work, but we present the analytical sampling variance of the estimated and its corresponding test statistics after correction of some technical errors in their original work. We can consequently evaluate how influences the estimation of heritability and its corresponding score, and, as data can be very large, the control of is of theoretical as well as practical importance. An analytical resolution crystallizes a computational procedure, and we further extend the method to another two new scenarios, called vertical-RHE-reg, which is a global implementation for LD score regression [15], and horizontal-RHE-reg, which enables Federated Learning but we estimate heritability in distributed data without compromise of privacy [14].

A framework for Randomized Haseman-Elston regression (RHE-reg)

In essence, Haseman-Elston regression is a kind of method of moments (MoM) estimation for heritability, and can provide equivalent estimates of heritability for complex traits after IBD is replaced with IBS [4,11]. As we extend the work of Wu and Sankararaman [5], we similarly assume that

in which is the standardized phenotype of the traits of interest, is the standardized genotype matrix of individuals, is the number of double allelic markers, is the cumulative effect related to each of the markers, is the residual effect, is an identity matrix, is the SNP heritability, and is an identity matrix, is the residual variance. Under the general assumption for a polygenic trait, it is easy to see that

is the genetic relationship matrix (GRM); the moment estimator, or randomized Haseman-Elston regression, is to minimize . Of , by differentiating and , respectively, we have the following normal equations:

The preliminary estimators for and are given as

For ease of discussion, we now only focus on the expression without adjustment of covariates. The denominator involves , a high-order function for GRM. Alternatively, according to the trace property of a matrix, it can be calculated that , a summation of the square of each element in . We proved that the expectation of , where is the effective number of markers that depicts the average squared Pearson’s correlations among all genomic markers as often used for measuring linkage disequilibrium [12,16,17]; a brief sketch of how can be transferred into is also presented in the section "Estimation for the effective number of markers" below. Therefore, the expectation for the preliminary estimator of is for a typical polygenic trait as established [4,18]; is the averaged LD between a marker and a causal variant, and is the averaged LD between any pair of markers – including the LD of a marker with itself. At first glance at Eq 2, it seems inevitable to compute , the computational cost of which is , a substantial cost given a large sample size, such as for UKB of about 500,000 samples. We obtain the estimate of according to the properties of matrix algebra, and is the exponential index and takes the value of 1, 2, 3, or 4 upon the application in this study.

Where is a linear estimator for , is a vector of length and each element of is sampled from the standard normal distribution, and is the round of iterations. Of note, the sampling variance of . As will be shown below, will be a plugin parameter in the analysis below, and we suggest a robust estimation of from rather than . Eq 3 is the most innovative part in the work of Wu and Sankararaman, and it is known as Girard-Hutchinson estimation for stochastic trace estimation [19, 20]. Of note, , which was incorrectly derived as in Wu and Sankararaman’s work [5], and to fix their problem directly led to the present work.

Randomized estimation for via RHE-Reg

When there is random mating, , and substituting the expressions given as Eq 3 into Eq 2, a randomized estimator of heritability is

The component in the denominator is no other than a shuffling nature of the estimation with rounds of resampling.

Sampling variance of RHE-reg

Of Eq 4, we have and , and their respective mean and variance are

The randomized estimator of can be seen as a ratio of , in which both and are variables, and, according to Delta method, its sampling variance can be expressed as , in which the covariance term can be zeroed out in this scenario [21]. So, we can obtain

For the definition of please refer to the section "Estimation for key parameters". As is a random variable, using Taylor approximation can be obtained by . .

in which the second term is the bias of the RHE-reg estimator. At the same time, we can also find the mean squared error (MSE), the summation of the sampling variance and squared bias, for as below

In this polynomial expression, as will be shown in the simulation and real data analysis, is largely upon the sampling variance, which can be further reduced with sufficient iterations (). As will be shown for UK Biobank examples, dynamically ranges from 10 to 200, even greater upon many factors.

Constructing test statistics

Given the estimation of heritability, we can construct the z-score statistic below:

in which , a quantity that will be zeroed out after sufficient iterations, and can be estimated from Eq 5. Obviously, when is large enough, the optimal z score is the following:

There is an obvious relationship between two z scores in Eq 8 and Eq 9 (practically ). Given we can predict optimal test statistic as below:

It means that after iteration the expectation of the test statistic is predictable in certain degree.

In summary, given iterations the test statistic observed is , which is subject to the realized values of , , and . is the expected optimized test statistic when is very large and zeroed out all uncertainty due to iteration. is a reconstruction of , and , indicating how a larger seem to bring out advantage in such as a more significant p-value.

Estimation for key parameters

There are several key quantities/parameters involved in the above equations for RHE-reg, and we present how to estimate them. These parameters are – effective number of markers, the trace of fourth-order GRM, and .

Estimation for the effective number of markers ()

can be decomposed into four terms upon (or ) and (or ), and according to Isserlis’s Theory [22], having integrated these four terms, we have

in which the effective number of markers and the squared Pearson’s correlation of LD between a pair of SNPs [23]. Often , and if all markers are in linkage equilibrium (see the note of Table 1 ). Here, is a population parameter, a summary statistic that encompasses allelic frequencies and linkage disequilibrium of makers. According to Eq 3, , we consequently propose a randomization algorithm, which estimates as below

A more detailed estimation procedure for can be found in our recent work [16]. See Note I in S1 Text for more details.

Estimation for

The benchmark estimation for is , a fourth-order summation of the eigenvalues of . However, it is computationally expensive when is large. There are two alternative choices to estimate . Method I: , and would affect its precision. Method II: , which uses the fourth-order randomized estimation in Eq 3. Both Method I and Method II can be realized via Eq 3. As will be shown below, Method II provides more stable estimates than Method I.

About —high-dimension structure of genetic architecture

For , , in which if is replaced by because is too expensive to constructed as aforementioned, we consequently have

can be estimated as in Eq 3 if is replaced by the phenotype itself; it is similarly for and . They reflect high-dimensional structure between and . So the sampling variance of is not only related to itself, but is eventually upon the high-order structure between and . See Note II in S1 Text for more discussion about .

About — the term determines the iteration B

We define the ratio as below

in which can be estimated as as above. However, it should be noticed that is a heavy penalty for higher heritability; for example, comparing with , leads to a 100-fold penalty for the latter in the numerator of Eq 13. Easily, we can estimate if we want to know how many iterations are needed to reach the preset ratio of

In practice, can take the value of 0.1 or 0.05 as in our simulation and real data analysis below.

Extended utilities for distributed GWAS datasets

Because datasets are often distributed across institutes, we consequently consider two scenarios for the application of RHE-reg in distributed datasets. As the estimation for (Eq 4) can be split into the numerator and the denominator, the numerator and the denominator are estimated from two different sources. In the other scenario, the whole dataset has been distributed into small slices at different institutes. We call the first scenario the vertical RHE-reg and the latter horizontal RHE-reg.

Vertical RHE-reg

Estimation for can be implemented in summary statistics that the numerator and the denominator can be from different components [18]. We denote the correspondingly heritability for this subtle difference, as well as all tilded symbols from a reference panel that is related to genotypes. Alternatively, Eq 4 can be rewritten as

the denominator . , in which has the dimension of ; is the genotype matrix of the reference sample that is employed to estimate . So, is (assuming )

The bias is , which will be zeroed out when increases (see Note III in S1 Text for more general situations). The corresponding test statistic is

in which . For the population of similar ancestry, the ratio is cancelled out after sufficient iteration, and leads to

Horizontal RHE-reg

For this application, it is assumed that the entire dataset is divided into institutes ( is subscript for and ). Consequently, , the whole data and are distributed in institutes, and the length of upon how the proportion of data has in the institute; , similarly the dimension of is in which is the number of individuals in the institute. One only needs to receive the mean and summation of square for each , and similarly for receiving the allele frequencies of the reference alleles of . So after scaling for and ,

, a matrix, can be generated from , by each institute, and consequently independently generate and without compromise of privacy; the subscript indicates Frobenius norm of a matrix that . Upon the precision requirement, after rounds of iterations, can be calculated so as to evaluate whether further iterations are needed. Unlike the vertical RHE-reg, the horizontal RHE-reg is identical to the RHE-reg under this simple scenario. An R script is attached for its detailed implementation (S1 Data).

Summary for RHE-reg

Now we discuss some computational issues about RHE-reg. So, eventually will creep into the RHE-reg. The focus here is to investigate how would affect the RHE-reg, in particular the stability of and z scores. All the above analyses are based on three computational units, , , and – if covariates are taken into account, and the operation between them lead to the whole computational procedure, of which their elementary operations can be implemented hierarchically (Table 2). We give an atlas for the computational route. Furthermore, we have covariates, and the covariate matrix is of dimensions. After inclusion of the covariates, the equations for stopping rules can be updated accordingly (see Note IV in S1 Text). We finish the description of the statistical approaches and go to their applications now.

Software

We have developed computer software that can handle biobank-scale algorithm presented in this study. The software reads genotype in binary format as defined in such as PLINK. For fast vector-matrix multiplication, Mailman algorithm is employed here [24]. We adapt the implementation of the Mailman algorithm from Agrawal’s fast PCA project [25]. It is known that using the Mailman algorithm the vector-matrix multiplication in is reduced from to . There is no conceptual obstacle to applying the method for genotype data in dosage format, but the Mailman algorithm cannot proceed in such a scenario. There are many matrix multiplication included, and in programming some suggested tips to take them out is as shown (S1 Fig).

Simulation results

We conducted simulations to evaluate the aforementioned theoretical results under various parameters. The reference allele frequency was evenly sampled from 0.1 ~ 0.5, and was set three values of 0, 0.1, and 0.25, and all SNPs were considered causal after a typical polygenic model, which follows Normal distribution. 1) The linkage disequilibrium (Lewontin’s ) for each pair of consecutive SNPs were 0, 0.2, 0.4, 0.6, and 0.8 for consecutive SNPs. 2) We set three levels of unrelated samples 1,000, 5,000, and 10,000, respectively. 3) Three levels of SNP numbers 10,000, 50,000, and 100,000. These five parameters could totally carry out 45 simulation scenarios for each by our in-house simulation code, and its detailed implementation can be found in Zhang et al [16]. For each simulation scenario, we set the value of 10, 20, and 50 in order to find proper . , (as well as their allele frequencies), , and were considered to investigate how to determine . Although neither nor reaches real biobank-scale data, we investigate and summarize certain properties of RHE-reg under these 135 scenarios in the results below. The biobank-scale test is to be investigated in UK Biobank examples.

Result 1: Randomized estimation for

As shown in the method section, was appeared as one of the key parameters in determining the performance of the sampling of RHE-reg. The direct estimation of from the eigenvalues of was the golden standard, and we consequently compared Method I, , and Method II, , with its direct estimation. As shown in Fig 1, the above 135 simulation scenarios were compared with the direct estimation for . For Method I, increasing from 10 to 50 could increase the precision of the estimation. In contrast, Method II showed very consistent and high precision for the estimation of regardless of the sample size, an increasing of from 10 to 50 did not help improve precision. The advantage of Method II was probably because estimated as its mean, whereas as its sampling variance. So, hereafter we used Method II to estimate . Of note, as is computational expensive when is large so that only limited sample size and SNP numbers were tested in these 135 simulations; however, the principal results should be retained for an even larger sample size, as well as , but with more expensive computational cost in solving eigenvalues. Biobank-scale performance of the proposed method will be illustrated in UK Biobank 81 traits in Result 4.

Result 2: MSE of RHE-reg

In Eq 7, , and we defined according to Eq 14. In Fig 2, we showed how MSE and could be reduced by for these 90 simulated scenarios (excluded 45 scenarios under ). We only illustrated the results for and 10,000, respectively, because was too small a sample size here for efficient convergence. The top row of Fig 2 illustrated how MSE were reduced by , and obviously a much larger reduced MSE because was turned down. Actually the bias term played little weight in MSE, which was dominantly determined by . was at least one or two order of magnitude compared with . In the second row of Fig 2 reflected how quickly vanished after iterations. Neither LD nor played an important role in determining MSE for the simulated scenarios, but the ratio between and mattered much as under the same sample size, more SNPs always inflated MSE.

Result 3: Randomized estimation for and z-score

In result 3, we studied how could influence and its z-score. As the sampling variance of was reciprocal to the sample size under the null hypothesis and , it was obviously to see in simulation that: greater , and greater would help to bring out a more stable estimation for (Fig 3A-C). If we employed from as the benchmark, when sample size , there was very high consistent estimation for even (Fig 3C). is the sampling variance of REML when [26].

The availability of the score of the estimated heritability was important for statistical inference. We evaluated the influence of in determining the performance of the randomized algorithm (Fig 3D-F). It was known from the above analysis , so when the estimation of became stable the test statistic was stable too. So, was relative stable when (Fig 3E) or (Fig 3F). When the sample size was sufficiently large, a few iteration could guarantee high accuracy of the estimation. In addition, we also tested the estimation by setting , 0.75, and 0.9, respectively, and the results, as promised by our theory, were consistent to what observed.

Result 4: Application of horizontal RHE-reg

This study was to estimate heritability for distributed data as exact as a single piece of data. Two cohorts with 4,000 and 6,000 individuals, respectively, were generated to verify h-RHE-reg. was set the value of 0, 0.1, and 0.25, respectively. The effects of all 10,000 SNPs were sampled from the distribution . Heritability and z scores were estimated using individual-level RHE-reg as well as h-RHE-reg. was set of 10, 20, and 50. The genotypes of the two simulated cohorts were standardized by and for the -th locus, where was the average allele frequency. The phenotypes of the two cohorts were standardized by and , where and . Each simulation scenario had 10 repeats (see Source1.R and Source2.R in S1 Data for its implementation).

The estimates of heritability and its z score were consistent using individual-level RHE-reg and h-RHE-reg in all scenarios when the random vectors were the same (Fig 4). We also split the data into 2,000 and 8,000 individuals, and, as expected, the results were nearly identical and unbiased.

Real data analysis for UK Biobank

We chose the unrelated 292,223 British white who have no kinship found, as indicated by the genetic kinship provided in the UK Biobank (field 22021) for real data test [2]. After quality control, the inclusion criteria were: MAF > 0.01, missing call rate < 0.05 and Hardy-Weinberg proportion test p-value > 1e-6, whose genotype call rate > 0.95, and 525,460 autosome SNPs were included for analysis. We estimated heritability of the 81 quantitative traits, and included the top two principal components and sex as covariates.

We used two strategies to estimate heritability. In strategy I, denoted as B+ strategy hereafter, we set as a warm-up step to evaluate and was set of 0.05. After the warm-up of iteration, we then increased iteration by a step of 10, We then estimate final realized , , , and three kinds of scores until the convergence ratio of ; however, we set a hard stop for even if was still greater than 0.05. In strategy II, we directly set , 20, or 50 without further considering additional iteration anymore, and consequently denoted as , , and strategies hereafter.

Comparison of the UKB results between strategy I and II

Even a couple of traits were set to take a hard stop because their were greater than 200, the estimated for the 81 traits had a mean of 0.0518, which was very close to the preset (Fig 5A1). It indicated that our theory worked to control the precision of the sampling variance of RHE-reg. The trait "age of diabetes diagnosed" had , extremely large standard error compared to other traits, because of its smallest sample size of (Fig 5B1). In Fig 5C1, we got three z scores, which are score directly calculated given iterations (green colored, via Eq 8), the optimal z score when was infinitive (blue colored, via Eq 9), and the predicted score (pink colored, via Eq 10).

For comparison, we examined the corresponding statistics that were estimated under , , and , respectively. In strategy II, A larger led a smaller as expected (Fig 5A2-4). Interestingly, regardless the change of in strategy II, were very consistent to those estimated from strategy I, as shown that the fitted regression lines were very close to 1 (Fig 5B2-B4). Three types of z scores were compared (Fig 5C2-C4), and the optimal z scores from both strategies were nearly perfect (blue points and blue dashed lines). Then, as shown in Fig 5C4, the three kinds of z scores were nearly completely matched.

In addition, the estimates were also consistent with our previous results using a less efficient method [27], and see S1 Table for more details. The heritability estimated by the randomization algorithm exhibited a relative high degree of correlation (Pearson’s correlation coefficient of 0.77) with the previous estimates for 81 traits. Compared to the previous results, the was nearly consistent with the GRM-based estimates, and is with averaged 1.38% deviation after 10 iterations and further decreased to 1.23% deviation after 50 iterations (S1 Table).

We also compared the computational efficiency of RHE-reg with GCTA [28] and BOLT-REML [29] in estimating the heritability on BMI. The comparison was conducted on a sub-dataset in UKB with randomly selected 10,000 individuals and 523,945 SNP markers after filtration. The results indicate significant efficiency improvement in estimating the heritability of complex traits in biobank-scale datasets for RHE-reg, with computation times reduced by 96.6% and 83.8% compared to GCTA and BOLT-REML, respectively (S2 Table). More benchmark comparison of the computational performance could be found in earlier studies [18,5]. Even using a complete dataset, RHE-reg could also complete heritability estimation within an acceptable time (S3 Table). In our tested 81 UKB traits, with 10 threads, it on average took 453 mins to finish the analysis of a trait and the average iteration of .

Application of vertical RHE-reg

Of Eq 15, indicates that and can be from two independent sources. Consequently, we split each UKB trait evenly into halves to test the v-RHE-reg, and Eq 15 had four possible combinations: 1) split 1/1: both and were estimated from split 1; 2) split 2/1: was estimated from split 2 and split 1; 3) split 1/2: was estimated from split 1 and split 2; 4) both and were estimated from split 2. So, we had four estimators as below

Of each trait, its heritability and z score tests could be constructed within each split and between each split by exchanging the estimation, and consequently brought out v-RHE-reg. As shown in Eq 15, we compared the result for =10, 20, and 50, respectively, and observed consistent results between split 1 and split 2, and between split 1/2 and split 2/1.

Fig 6 showed the results of these four estimators under different . It illustrated that pairwise estimates against , and against , and as observed the pairwise estimates were quite consistent with each other both within and between splits.

S1 Text. Technical details on some mathematical derivations.

"Effective number of markers" (Note I); "Discussion of " (Note II); "Sampling variance for vertical RHE-reg" (Note III); "Adjustment for covariates" (Note IV); "Coding scheme and LD" (Note V).

https://doi.org/10.1371/journal.pcbi.1013568.s001

(DOCX)

S1 Fig. This figure gives the computational tips on how to program the software.

Some sequential operation of the matrix is suggested to make the program easy to write.

https://doi.org/10.1371/journal.pcbi.1013568.s002

(TIFF)

S1 Data. Implementation for Horizontal RHE-reg (R code: Source1. R and Source2.R).

https://doi.org/10.1371/journal.pcbi.1013568.s003

(ZIP)

S1 Table. SNP-heritability estimation for 81 UKB traits (XLSX).

https://doi.org/10.1371/journal.pcbi.1013568.s004

(XLSX)

S2 Table. The time cost for heritability estimation on BMI for RHE-reg, BOLT-REML and GCTA.

https://doi.org/10.1371/journal.pcbi.1013568.s005

(XLSX)

S3 Table. Computational time for 81 continuous traits using RHE-reg.

https://doi.org/10.1371/journal.pcbi.1013568.s006

(XLSX)