A novel method for an unbiased estimate of cross-ancestry genetic correlation using individual-level data

Cross-ancestry genetic correlation is an important parameter to understand the genetic relationship between two ancestry groups for a complex trait. However, existing methods cannot properly account for ancestry-specific genetic architecture, which is diverse across ancestries, producing biased estimates of cross-ancestry genetic correlation. Here, we present a method to construct a genomic relationship matrix (GRM) that can correctly account for the relationship between ancestry-specific allele frequencies and ancestry-specific causal effects. Through comprehensive simulations, we show that the proposed method outperforms existing methods in the estimations of SNP-based heritability and cross-ancestry genetic correlation. The proposed method is further applied to six anthropometric traits from the UK Biobank data across 5 ancestry groups. One of our findings is that for obesity, the estimated genetic correlation between African and European ancestry cohorts is significantly different from unity, suggesting that obesity is genetically heterogenous between these two ancestry groups.


Introduction
Complex traits are polygenic and influenced by environmental factors 1,2 , which can be distinguished from Mendelian traits that are regulated by single or few major genes and minimal environmental influences. In humans, causal loci and their effects on complex traits are dynamically distributed across populations such that the same trait can be genetically heterogenous between two ancestry groups. Studies of cross-ancestry genetic correlation are therefore important to understand the genetic relationship between two ancestry groups when investigating a complex trait 2 . Understanding genetic correlation between two ancestry groups can inform us on how well we can predict a complex disease in one ancestry group based on the information on another.
Genome wide association studies (GWAS) have been successful in identifying variants associated with causal loci, and can also provide a useful resource to investigate crossancestry genetic correlations. Using GWAS datasets, genomic relationships between samples across multiple ancestry groups can be constructed based on genome-wide single nucleotide polymorphisms (SNPs), which can be fitted with phenotypes across ancestries in a statistical model 3,4 . This paradigm-shift approach to estimate cross-ancestry genetic correlations have greatly increased the potential to dissect the latent genetic relationship between ancestry groups for complex traits.
Most GWAS have focused on European ancestry samples (> 80%) [5][6][7] although Europeans represent only 16% of the global population [8][9][10] . Because of large samples, the estimated SNP associations in Europeans are far more accurate than those in other ancestries. As a matter of fact, the performance of polygenic risk prediction depends on the accuracy of estimated SNP associations, causing a disparity in genetic prediction across populations 5 . Therefore, crossancestry genetic studies are urgently required to bridge the disparity, e.g. estimated crossancestry genetic correlations may be able to inform if SNP effects estimated in European ancestry samples can be useful to predict the phenotypes of samples in other ancestry groups.
To date, several studies have been undertaken to estimate cross-ancestry genetic correlations between diverse ancestry groups for a range of complex traits 2,3,11-13 . For example, Yang et al. 4 estimated the cross-ancestry genetic correlation between East Asians and Europeans for ADHD, where ancestry-specific allele frequencies were used to standardise samples' genotype coefficients in estimating their genomic relationships. This method of cross-P a g e | 3 ancestry genomic relationships has been widely used for cross-ancestry genetic studies 3,14 .
However, the method cannot account for the trait genetic architecture specific to each ancestry group, which is a function of the relationship between the genetic variance and allele frequency (one important aspect of a heritability model 15,16 ). With an incorrect heritability model, estimated genetic variances within and covariance between ancestry groups are biased 15,16 , and hence the cross-ancestry genetic correlations cannot be correctly estimated.
Another method based on GWAS summary statistics has been introduced (Popcorn) 2 .
However, this method also cannot correctly account for diverse genetic architectures across ancestry groups.
Here, we develop a novel method that can properly account for ancestry-specific genetic architecture and ancestry-specific allele frequency in estimating a genomic relationship matrix (GRM). In addition, we revisit the existing theory to correct mean bias in genomic relationships. In simulations, the SNP-based heritability and cross-ancestry genetic correlation estimated from our proposed method are shown to be unbiased, whereas other existing methods generate biased estimates. We apply the proposed method to six anthropometric traits from the UK Biobank data that include standing height, body mass index (BMI), waist circumference, hip circumference, waist-hip ratio, and weight. For each trait, we estimate SNP-based heritabilities and cross-ancestry genetic correlations across 5 ancestry groups, i.e., white British, other Europeans, Asian, African, and mixed ancestry cohorts.

Overview
We used publicly available data from the UK Biobank. Participants of the UK Biobank were stratified into multiple ancestries (white British, other European, Asian, African, and mixed ancestry cohorts) according to their underlying genetic ancestry based on a principal component analysis (see Supplementary Figure S1) 17 . In each ancestry group, we assume that the relationship between genetic variance and allele frequencies (heritability model)

‫ݒ‬
, and α is the scale factor determining the genetic architecture of complex traits in each ancestry groups. Note that with α = 0, equation (1) becomes the classical model of Falconer and Mackay (1996) 18 . By assuming that the genetic variance of causal variant is constant across minor allele frequencies (MAF) spectrum, a heritability model with α = -0.5 has been widely used [19][20][21]Speed et al. 15,16,22 reported a different α value, e.g., α = -0.125 for anthropometric traits, using multiple European cohorts. While it is intuitive that α values can be varied across ancestry groups, it has not been well studied.
First, we determine an optimal α value for each ancestry group, comparing model fits (maximum log-likelihood) of various heritability models with different α values for the 6 anthropometric traits from the UK Biobank (see Methods). Second, we simulate phenotypes based on the UK Biobank genotypic data to assess if SNP-based heritability and crossancestry genetic correlation are unbiasedly estimated. In the simulation, various α values are used to generate causal effects of SNPs in various ancestry groups, and the correlation of SNP effects between ancestry groups varies between 0 and 1 (see Methods). For simulated phenotypes of multiple ancestry groups, we estimate SNP-based heritability and crossancestry genetic correlation, using bivariate GREML 20,21,23 with four existing methods to construct GRM as shown in Table 1. In addition to the existing methods, we use a novel method to estimate GRM in the cross-ancestry analysis in which we use ancestry-specific α value and ancestry-specific allele frequency, so that the estimation model matches with the ancestry-specific genetic architecture (see Methods). The equation for the proposed method can be written as  Table 2).
Finally, we analyse real data using the proposed method (equation 2) to estimate SNP-based heritability and cross-ancestry genetic correlation for 6 anthropometric traits across different ancestries using bivariate GREML. (10) GRM4 -0.5 common only ancestry-specific 3,4 (10) Proposed method ancestryspecific e common only ancestry-specific (2) a Using all SNPs from both ancestry groups. b Allele frequencies estimated from the combined population of both ancestry groups when scaling the genotypes. c Using only common SNPs presented for both ancestry groups after QC including ancestry specific MAF > 0.01. d Allele frequencies estimated from each population when scaling the genotypes. e Different α value specific to each ancestry group can be used when scaling the genotypes.

Determination of scale factor (α) across ancestries
We compared the Akaike information criteria (AIC) values of heritability models with varying α values to determine which α value provides the best model fit (see Methods), which is analogue to the approach of Speed et al. 16 Figure 2 and Supplementary Tables 3-7).
When comparing GCTA-α and LDAK-thin-α models, the AIC value of GCTA-α model was much smaller than that of LDAK-thin-α model for white British or other European ancestry cohort (Supplementary Table 8). For Asian ancestry cohort, the AIC of GCTA-α model was slightly lower than LDAK-thin-α model when using the best α value = -0.625. In contrast, the AIC of LDAK-thin-α model was generally lower than that of GCTA-α model for African or mixed ancestry cohort. Δ AIC values from GCTA-α models are plotted against scaling factors, α , for each ancestry group. The lowest AIC (i.e. Δ AIC=0) indicates the best model. The sample sizes are 30,000, 26,457, 6,199, 6,179 and 11,797 for white British, other European, Asian, African, and mixed ancestry groups, respectively.

Method validation by simulation
We simulated phenotypes based on the real genotypic data of multiple ancestry groups where the estimated α value for each ancestry group was used to generate SNP effects that were correlated between ancestry groups (Methods). In this simulation, we do not consider associations between the causal effects and LD structure for any SNP, i.e. LDAK simulation model, because LDAK model was not particularly plausible for the genetic architecture of the traits especially for white British, other European and Asian ancestry cohorts (Supplementary Table 8) and LDAK simulation model was not feasible for a bivariate framework. When using simulation models with α = -0.5 for all ancestry groups, estimated SNP-based heritabilities were mostly unbiased for all the methods, GRM1-4 (Supplementary Table 9). Estimated genetic correlations from GRM3 and 4 were unbiased for all scenarios (Figure 2 and Supplementary Table 9). However, estimated genetic P a g e | 7 correlations from GRM1 and 2 were biased when the true genetic correlation was high α is fixed to -0.5 for both simulation and estimation. Height of each bar is the estimated crossancestry genetic correlation and error bar indicates 95% confidence interval (CI) for 500 replicates. In our simulations, we used three combinations for estimating cross-ancestry genetic correlation (white British vs. Asian, white British vs. African and Asian vs. African ancestry cohorts). In each P a g e | 8 ancestry group, we used 500,000 SNPs that were randomly selected from HapMap3 SNPs after QC.
To simulate phenotypes, we selected a random set of 1,000 SNPs as causal variants, which were presented for both ancestry groups. We used α = -0.5 when scaling the causal effects by ancestryspecific allele frequency in each ancestry group. Various values of genetic correlation were considered (0, 0.25, 0.50, 0.75 and 1.0). In the estimation, the four methods (GRM1 -4) used α = -0.5 (standard scale factor in GRM estimation). GRM1 and 3 used all available SNPs from both ancestry groups (791581, 812332 and 777894 for figure panels a, b and c) whereas GRM2 and 4 used only the set of SNPs common between two ancestry groups (208419, 187668 and 222106 for a, b, and c). When scaled with α = -0.5, GRM1 and 2 used allele frequency averaged between two ancestry groups whereas GRM3 and 4 used ancestry-specific allele frequency estimated from each ancestry group (see Table 1).
To mimic real data, we simulated phenotypes based on the real genotypes of multiple ancestry groups, using realistic α values (α = -0.25 for white British ancestry cohorts, α = -0.625 for Asian ancestry cohorts and α = -0.75 for African ancestry cohorts) instead of using a constant α value across ancestries (Methods). For these simulated phenotypes, estimated cross-ancestry genetic correlations using four existing methods were biased when the true genetic correlation was 0.5 or higher between white British and African ancestry cohorts or between Asian and African ancestry cohorts. Biased estimates were still observed even when ancestry-specific allele frequencies were considered in GRM3 and 4 (Figure 3b and 3c For the same simulated phenotypes, we applied the proposed method that accounts for ancestry-specific allele frequency and ancestry-specific α values (equation 2) and found that it provided unbiased estimates for both SNP-based heritability and cross-ancestry genetic correlation (Figure 3 and Supplementary  Height of each bar is the estimated cross-ancestry genetic correlation and error bar indicates 95% confidence interval (CI) for 500 replicates. In our simulation, we used three combinations for estimating cross-ancestry genetic correlation (white British vs. Asian, white British vs. African and Asian vs. African ancestry cohorts). In each ancestry group, we used 500,000 SNPs that were randomly selected from HapMap3 SNPs after QC. To simulate phenotypes, we selected a random set of 1,000 SNPs as causal variants, which were presented for both ancestry groups. We used various α values that were specific to ancestries (α = -0.25, -0.625 and -0.75 for P a g e | 1 0 white British, Asian and African ancestry cohorts, respectively) when scaling the causal effects by ancestry-specific allele frequency in each ancestry group. Various values of genetic correlation were considered (0, 0.25, 0.50, 0.75 and 1.0). In the estimation, we used existing methods (GRM1 -4) that used the standard scale factor α = -0.5 in GRM estimation. GRM1 and 3 used all available SNPs from both ancestry groups (791581, 812332 and 777894 for figure panels a, b and c) whereas GRM2 and 4 used only the set of SNPs common between two ancestry groups (208419, 187668 and 222106 for a, b, and c). When scaled with α = -0.5, GRM1 and 2 used allele frequency averaged between two ancestry groups whereas GRM3 and 4 used ancestry-specific allele frequency estimated from each ancestry group. In addition, we applied the proposed method that used ancestry-specific α value and ancestry-specific allele frequency in GRM estimation.
Popcorn software, a GWAS summary statistics-based method for cross-ancestry genetic analysis, also generated biased estimates of cross-ancestry genetic correlation when using realistic α values in the simulation (Supplementary Figure 4). This agrees with Brown et al. 2 , that popcorn estimates can be deviated from the true values when using alternative genetic architectures (or heritability models).

SNP-based Heritability (h 2 ) estimates for anthropometric traits using real data
The estimated SNP-based heritability for each of the anthropometric traits was presented in

Cross-ancestry genetic correlations for anthropometric traits
Estimated cross-ancestry genetic correlations (r g ) between ancestry groups for 6 anthropometric traits are shown in Figure   Estimated cross-ancestry genetic correlations for waist circumference and hip circumference showed a similar pattern with BMI, i.e., there was a significant evidence for genetic heterogeneity between white British and African ancestry cohorts and between other Europeans and African ancestry cohorts (Figure 5c and 5d; Supplementary Tables 15 and   16). The estimated genetic correlation between African and mixed ancestry cohorts was low, but not significantly different from 1. As the same as in BMI and height, there was no genetic heterogeneity between white British and other European ancestry cohorts for both waist circumference and hip circumference.
For weight, the estimated genetic correlation between African and other European ancestry cohorts was significantly different from 1 ( ‫ݎ‬ =0.624; SE=0.139; p-value = 6.83e-03), indicating a significant genetic heterogeneity between these two ancestry groups (Figure 5f and Supplementary Table 18). Although the estimations are not significant, we have estimated lower genetic correlations (far from 1) between white British and African ancestry cohorts, between African and mixed ancestry cohorts and between white British and Asian ancestry cohorts.
We did not observe any significant heterogeneity across ancestries (genetic correlation estimate was not significantly different from 1) for waist-hip ratio ( The colour and size of each pie chart indicates the magnitude of estimated cross-ancestry genetic correlations. The value in each pie chart is a p-value (*, ** and *** indicates p value < 0.05, < 0.01 and <0.001, respectively) testing the null hypothesis of ‫ݎ‬ =1. Coloured asterisk indicates significantly different from 1 after Bonferroni correction (0.05/54). Non-estimable parameter is shown as NA, which was due to one ancestry group is nested within the other ancestry groups or estimated SNP-h 2 is zero for one ancestry group.

Discussion
We propose a novel method that provides unbiased estimates of ancestry-specific SNP-based heritability and cross-ancestry genetic correlations. This is possible because the proposed method correctly account for ancestry-specific genetic architectures or ancestry-specific heritability models. Our method provides a tool to dissect the ancestry-specific genetic architecture of a complex trait and can inform how genetic variance and covariance change across populations and ancestries. By using a meta-analysis across multiple ancestry groups 25,26 based on unbiased estimates of ancestry-specific heritability and cross-ancestry genetic correlations, we hope the current ancestry disparity and study bias in GWAS 5,8 can be reduced. Cross-ancestry genetic correlations can provide crucial information in cross-ancestry GWAS and cross-ancestry polygenic risk score prediction. We showed a significant genetic heterogeneity for obesity traits (BMI, weight and waist and hip circumferences) between African and European ancestry samples. For height, there was a significant genetic heterogeneity between African and Asian ancestry samples. However, without considering ancestry-specific α values (i.e. GRM 4 from equation 5), the findings were changed and could be over-interpreted, e.g. there was an additionally significant genetic heterogeneity between There are several limitations in this study. Firstly, we used the GCTA-α model only in the real data analysis, assuming that all SNPs contributed equally to the heritability estimation 21,33 . We did not use the LDAK-thin-α model that required to prune SNPs within each ancestry group, which could substantially reduce the number of common SNPs between two ancestry groups in cross-ancestry genetic correlation analyses. Secondly, we did not consider MAF stratified, LDMS, or baseline model 34-36 in estimating cross-ancestry genetic correlations because it was developed to estimate SNP-based heritability, but not suitably designed to estimate genetic correlations. Furthermore, it is required to fit multiple random effects (i.e. multiple GRMs), which is not computationally efficient. Nevertheless, these advanced models can improve the estimations of ancestry-specific SNP-based heritability and cross-ancestry genetic correlation. Thirdly, we estimated optimal scale factors (α) with a moderate sample size especially for Asian or African ancestry cohorts, resulting in a P a g e | 1 6 relatively flat curve of Δ AIC values. A further study is required to estimate more reliable α for Asian or African ancestry cohorts with a larger sample size.
In conclusion, we present a method to construct a GRM that can correctly account for the relationship between ancestry-specific allele frequencies and ancestry-specific causal effects.
As the result, our method can provide unbiased estimates of ancestry-specific SNP-based heritability and cross-ancestry genetic correlation. By applying our proposed method to anthropometric traits, we found that obesity is a genetically heterogenous trait for African and European ancestry cohorts, while human height is a genetically heterogenous trait between African and Asian ancestry cohorts.

Ethics statement
We used data from the UK Biobank (https://www.ukbiobank.ac.uk), the scientific protocol of which has been reviewed and approved by the North West Multi-center Research Ethics Committee, National Information Governance Board for Health & Social Care, and Community Health Index Advisory Group. UK Biobank has obtained informed consent from all participants. Our access to the UK Biobank data was under the reference number 14575.
The research ethics approval of this study has been obtained from the University of South Australia Human Research Ethics Committee.

Statistical model
Our main aim is to unbiasedly estimate cross-ancestry genetic correlation for a complex trait, using common SNPs presented for both populations. We use a bivariate linear mixed model (LMM) to estimate SNP-based heritability and cross-ancestry genetic correlation, using GWAS data comprising multiple ancestry groups. In the model, a vector of phenotypic observations for each ancestry group can be decomposed into fixed effects, random genetic effects and residuals. The LMM can be written as, where y i is the vector of phenotypic observations, b i is the vector of fixed (environmental) effect with the incidence matrix X i , g i is the vector of random additive genetic effects with the incidence matrix Z i and e i is the vector of residual effects for the i th ancestry group (i = 1 and 2).
is the genetic covariances between the two ancestry groups. Note that it is not required to model residual correlation in V because there are no multiple phenotypic measures for any individual (no common residual effects shared between two ancestry groups).

The variance of random additive genetic effects
Assuming that causal variants are in linkage equilibrium and that the phenotypic variance is

Genomic Relationship Matrix (GRM)
GRM is a kernel matrix and can be normalized with a popular form of (2008) where A ij is the genomic relationship between the i th and j th individuals, L is the total number where α is a scale factor, d is the expected diagonals, When α = -0.5, equation. (5) is equivalent to equation (3), and when α = 0, equation (5) becomes equivalent to equation ( where the expectation of the first term can be rewritten as , the diagonals from equation (5) or (7) can be expressed as which is deviated from 1 by a factor 1/n. Without loss of generality, the biased factor with any α values can be written as In a similar manner, the off diagonals can be written as where the expectation of the first term can be rewritten as , the off diagonals from equation (5) where d is redefined as is noted that with a sufficient sample size (> 1,000), the difference between equation (7) and (8)

GRM for cross-ancestry genetic analyses
Yang et al. 4 proposed a GRM method that uses ancestry-specific allele frequency to be applied to cross-ancestry genetic analyses, which can be expressed as It is intuitive that α value can be dynamically changed across ancestry groups because of genetic drift, natural selection, and various selection pressures. Combined with a revised formula derived above (equation (9)), a novel GRM equation for cross-ancestry genetic analyses, which accounts for ancestry-specific α and ancestry-specific allele frequency, can be proposed as

Data source and quality control
We tested this proposed method in real genotypic and phenotypic data obtained from the  Figure1). We did not include individuals who do not know their ancestry and who prefer not to answer (UK Biobank codes are -1 and -6). Gender mismatch and sex chromosome aneuploidy were also excluded during the quality control (QC) process.

Phenotype Simulation
The phenotypes of each ancestry group were simulated using a bivariate linear mixed model SNPs were randomly chosen for each ancestry group. To simulate phenotypes, we randomly selected 1000 SNPs that were common between the two ancestry groups and assigned causal effects to them. According to equation (4)

Determining scale factor (α) across ancestries for LDAK-thin-α and GCTA-α model
We analysed six different anthropometric traits (BMI, standing height, waist circumference, hip circumference, waist hip ratio and weight) from the UK Biobank across different ancestries (white British, other European, Asian, African, and mixed ancestry cohorts). These traits were adjusted for demographic variables, UK biobank assessment centre, genotype measurement batch and population structure measured by the first 10 principal components traits simultaneously. Note that in the multivariate model, we treated the six traits independent without considering residual and genetic correlations between traits. This was because the optimal alpha value should be trait-specific as our main analysis was to estimate trait-specific cross-ancestry genetic correlation, i.e. equal weights from all six traits to obtain the optimal α value for each ancestry group. Finally, we identified optimal α that gives lowest AIC values, i.e.
is the logarithm of the maximum likelihood from the model and k is the number of parameters.

Genetic correlation estimation using existing methods
Bivariate GREML 11,23,48,49 analyses were used to estimate heritability and cross-ancestry genetic correlation. In the analyses, we used four existing GRM methods ( Table 1) to assess their performance, compared to the proposed method using simulated phenotypes, using PLINK 42 (-make-grm-gz for command GRM1 and GRM2) and GCTA 21 software (-sub-popu command for GRM3 and GRM4). The distribution of diagonal elements and off-diagonal elements across ancestries represented in Supplementary Figure 6-10 (using α =-0.5 and estimated in PLINK2). We also used Popcorn 2,11 , that could be based on GWAS summary statistics. For these methods, we calculated empirical SE and its 95% confidence interval (CI) over 500 replicates to assess the unbiased estimation of the methods for each combination of ancestry pairs.

Web resources
The genotype and phenotype data of the UK Biobank can be accessed through procedures described on its webpage (https://www.ukbiobank.ac.uk/). Simulated data used in this paper can be obtained from the authors upon request.

Code availability
The source code for MTG2 and example code along with related files for estimating GRM in combined ancestries using MTG2 can be accessed from https://sites.google.com/site/honglee0707/mtg2

Competing interests
The authors declare that they have no competing interests.  G  e  n  o  m  e  -w  i  d  e  g  e  n  e  t  i  c  h  o  m  o  g  e  n  e  i  t  y  b  e  t  w  e  e  n  s  e  x  e  s  a  n  d  p  o  p  u  l  a  t  i  o  n  s  f  o  r  h  u  m  a  n  h  e  i  g  h  t  a  n  d  b  o  d  y  m  a  s  s  i  n  d  e  x  .   H  u  m  a  n  M  o  l  e  c  u  l  a  r  G  e  n  e  t  i  c