# Genomic Prediction Using Individual-Level Data and Summary Statistics from Multiple Populations

^{*}Wageningen University and Research, Animal Breeding and Genomics, 6700 AH, The Netherlands^{†}The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Easter Bush Research Centre, Midlothian EH25 9RG, UK

- 1Corresponding author: Wageningen University and Research, Animal Breeding and Genomics, P.O. Box 338, 6700 AH Wageningen, The Netherlands. E-mail: jeremie.vandenplas{at}wur.nl

## Abstract

This study presents a method for genomic prediction that uses individual-level data and summary statistics from multiple populations. Genome-wide markers are nowadays widely used to predict complex traits, and genomic prediction using multi-population data are an appealing approach to achieve higher prediction accuracies. However, sharing of individual-level data across populations is not always possible. We present a method that enables integration of summary statistics from separate analyses with the available individual-level data. The data can either consist of individuals with single or multiple (weighted) phenotype records per individual. We developed a method based on a hypothetical joint analysis model and absorption of population-specific information. We show that population-specific information is fully captured by estimated allele substitution effects and the accuracy of those estimates, *i.e.*, the summary statistics. The method gives identical result as the joint analysis of all individual-level data when complete summary statistics are available. We provide a series of easy-to-use approximations that can be used when complete summary statistics are not available or impractical to share. Simulations show that approximations enable integration of different sources of information across a wide range of settings, yielding accurate predictions. The method can be readily extended to multiple-traits. In summary, the developed method enables integration of genome-wide data in the individual-level or summary statistics from multiple populations to obtain more accurate estimates of allele substitution effects and genomic predictions.

- meta-analysis
- quantitative trait
- statistical method
- Genomic Prediction
- GenPred
- Shared Data Resources

GENOME-WIDE markers are nowadays widely used in animal and plant breeding to predict complex traits. This prediction is based on a linear model that partitions for each individual the observed complex phenotype value into systematic effects, comprising at least a population mean, an individual genetic value, and an environmental deviation (Fisher 1918). With genome-wide markers, individual genetic values can be computed from allele substitution effects estimated from individual-level phenotype and genotype data (Meuwissen *et al.* 2001). Subsequently, genetic values can be also computed for individuals of interest that are genotyped, but not phenotyped. This process is commonly called genomic prediction. In animal and plant breeding, genetic values are used to identify genetically superior individuals and use them as parents of the next generation to improve complex traits like milk yield (Meuwissen *et al.* 2001; VanRaden 2008) or grain yield (Schulthess *et al.* 2016). In human genetics, genetic values can be used to predict individual genetic risk for complex diseases to inform preventive and personalized medicine (de los Campos *et al.* 2010; Wray *et al.* 2013; Pasaniuc and Price 2017).

Accuracy of estimated allele substitution effects and of resulting genetic values for complex traits are foremost a function of the number of individuals with available phenotypes and genotypes (Daetwyler *et al.* 2008). To maximize the prediction accuracy, use of all available data are recommended (Henderson 1984; Wray *et al.* 2013; Vilhjálmsson *et al.* 2015). In some small populations, collecting large amounts of data are not possible, and a joint analysis across multiple populations is needed to achieve high accuracy (Hozé *et al.* 2014; Wientjes *et al.* 2016). However, such joint analysis is often impossible, because of logistic or privacy considerations (Powell and Norman 1998; Maier *et al.* 2018). Therefore, several methods were proposed to enable analysis of data from multiple populations when individual-level data are not available (Pasaniuc and Price 2017; Liu and Goddard 2018; Maier *et al.* 2018). These methods, often called meta-analyses (Pasaniuc and Price 2017), approximate a joint analysis by first obtaining summary statistics from separate analyses of individual-level data for each population, and then combining these summary statistics to estimate genetic values. In human genetics, summary statistics usually consist of publically available allele substitution effects, *i.e.*, genome-wide associations, together with their SE, estimated independently for each marker (Yang *et al.* 2012; Vilhjálmsson *et al.* 2015; Maier *et al.* 2018). In livestock, summary statistics more likely consist of allele substitution effects estimated jointly for all markers, together with prediction error (co)variances (Liu and Goddard 2018). While these methods may increase prediction accuracy in comparison to separate analyses, a loss in prediction accuracy is expected relative to an analysis using all individual-level data due to approximations (Maier *et al.* 2018). Further, these methods are based on some assumptions that make them difficult to apply outside their context of development. For example, Maier *et al.* (2018) implicitly assumed that only a single phenotype record per trait was associated with an individual. While this is usually the case in human genetics, it is not in breeding populations where individuals may have repeated phenotype records for the same trait, *e.g.*, repeated longitudinal production or reproduction records in livestock or replicated field trials in crops, or when phenotype records are measured on a group of individuals and linked to a genotyped relative, *e.g.*, progeny tested bulls for dairy production. Also, these developed methods do not allow combining individual-level data from some and summary statistics from other populations in one analysis (Liu and Goddard 2018; Maier *et al.* 2018).

The objective of this study was to develop a method that jointly analyses individual-level data and summary statistics from multiple populations with no, or a limited amount of, approximation. The method assumes that individual-level data are composed of marker genotypes and phenotype records that potentially have a variable number of replicates per individual. Further, summary statistics are assumed to be composed of estimated allele substitution effects with an associated measure of accuracy. Different measures of accuracy can be used, which controls the amount of approximation. The developed method is validated with simulated data. The results show that the method enables accurate integration of different sources of information across a wide range of settings.

## Materials and Methods

The first part of this section describes the theory of (1) separate and joint analyses of two individual-level datasets, (2) an exact integration of estimated allele substitution effects from one population into the analysis of another, (3) approximate integrations, and (4) generalization for multiple populations. The second part describes simulations used for validation of the developed method.

### Theory

Assume we have two populations with independent individual-level datasets of phenotyped and genotyped individuals. The two populations and their corresponding datasets are hereafter referred to as 1 and 2. Further assume that both datasets contain the same markers. From this data we want to obtain accurate estimates of allele substitution effects and genetic values for complex traits. We can achieve this by a joint analysis of the two datasets. When one of the datasets is not available, we can achieve this by integrating the results of a separate analysis of the unavailable data into the separate analysis of the available dataset. We show how to perform this integration exactly or approximately.

#### Separate and joint analyses:

A standard marker model, using random regression on marker genotypes, for the separate analysis of dataset i (i = 1, 2) is:(1)where is a vector of phenotypes, is a vector of fixed effects that are linked to by a incidence matrix , is a vector of allele substitution effects that are linked to by a incidence matrix and a matrix of genotypes and is the vector of residuals. In this work we consider single-nucleotide polymorphism markers, which we code in as 0 for homozygous aa, 1 for heterozygous aA or Aa, and 2 for homozygous AA. Other genotype coding and centering, that is of the form with 1 being a vector of ones and being a vector, can be used with no difference in obtained estimates of allele substitution effects (Strandén and Christensen 2011). We assume a prior multivariate normal (MVN) distribution for allele substitution effects for the separate analysis of the dataset *i*, , with mean zero and covariance where is a diagonal matrix (*e.g.*, an identity matrix I), and is the variance of allele substitution effects. We also assume that residuals are multivariate normally distributed with mean zero and covariance where is a diagonal matrix (*e.g.*, an identity matrix I), and is the residual variance. For simplicity, and without loss of generality, it is assumed in the following that residual variances are the same for all separate and joint analyses. Variance components and are assumed known, as they will have been estimated from the data previously. This marker model is the ridge regression model (Hoerl and Kennard 1976; Whittaker *et al.* 2000; Meuwissen *et al.* 2001; de los Campos *et al.* 2012) with optional different weights in (to differentially shrink different loci) and (to account for heterogeneous residual variance due to variable number of repeated phenotype records per individual).

Separate estimates of allele substitution effects are obtained by solving the following system of equations:(2)Separate estimates of genetic values for individuals in a dataset i (i = 1, 2) are obtained by

A marker model for the joint analysis of two datasets 1 and 2 is:(3)where phenotypes from the two populations are modeled with population-specific fixed effects but a joint set of allele substitution effects We assume a MVN prior distribution for allele substitution effects with mean zero and covariance where is a diagonal matrix, and is the variance of allele substitution effects in the joint analysis. We also assume that residuals are multivariate normally distributed, specifically where is a diagonal matrix.

Joint estimates of allele substitution effects (4) are obtained by solving the following system of equations:Joint estimates of genetic values for individuals in a dataset i (i = 1, 2) are obtained by

#### Exact integration:

The integration of estimates of allele substitution effects from one dataset into the analysis of another can be performed by means of absorbing corresponding equations in the joint system of equations. We choose to integrate estimates from the dataset 1 into the analysis of dataset 2. Derivations in Appendix A1 lead to the following system of equations that performs such integration and gives equivalent estimates of allele substitution effects to the joint analysis (Eq. 4):(5)where are estimates of allele substitution effects from the separate analysis of dataset 1 using (Eq. 2), and is the inverse of the corresponding prediction error covariance (PEC) matrix. The latter can be obtained as with . Note that only the individual-level dataset 2 and summary statistics from the dataset 1 (*i.e.*, the estimated allele substitution effects and their PEC) are required. Individual-level dataset 1 is therefore not required.

It is worth noting that the integration of estimates of allele substitution effects from the dataset 1 into the analysis of dataset 2 can also be obtained from a Bayesian context. Bayes estimators for linear mixed models were discussed by several authors (Lindley and Smith 1972; Dempfle 1977; Gianola and Fernando 1986). In a Bayesian context, we can assume the following prior multivariate normal distributions for the marker model (Eq. 1) applied to dataset 2:,where is a mean vector and is a (co)variance matrix,andAssuming a noninformative prior for the system of equations (2) for dataset 2 can be obtained by differentiating the joint posterior distribution of and with respect to and and setting the derivatives equal to 0 (Gianola and Fernando 1986). Integration of estimates of allele substitution effects from dataset 1 into the analysis of dataset 2 can be therefore obtained by defining a MVN prior distribution for allele substitution effects in the analysis of dataset 2 using the posterior distribution for allele substitution effects from a separate analysis of dataset 1:(6)The matrix Q can be considered as the PEC matrix of a hypothetical separate analysis of dataset 1 using the MVN prior distribution for allele substitution effects of the joint analysis, that is and and the vector can be considered as the estimated allele substitution effects of this hypothetical separate analysis. In animal breeding, a similar approach was used to integrate estimated genetic values and associated accuracies from one genetic evaluation into another genetic evaluation (Quaas and Zhang 2006; Legarra *et al.* 2007; Vandenplas and Gengler 2012).

Finally, it is worth noting that the term can be interpreted as a vector of hypothetical or pseudophenotype records associated with allele substitution effects, and, as such, summarize available information in dataset 1. In this sense, the system (Eq. 5) is similar to approaches that compute pseudorecords associated with individuals, from available estimated genetic values where individual-level phenotypic information is not readily available, or is not measured on the individuals themselves but on close relatives. In animal breeding, these approaches are commonly known as deregression of estimated genetic values (Jairath *et al.* 1998).

#### Approximate integration:

Exact integration requires the inverse of PEC matrix from the separate analysis, which could be approximated when unavailable. Genomic analyses of complex traits that combine different datasets commonly have access to estimated allele substitution effects and associated prediction error variances (in different forms), but not the whole PEC matrix required in (5). We propose several ways to accommodate this situation. We assume that we know, at least, the prediction error variances (PEV) of estimated allele substitution effects the number of individuals and variance components used in the separate analysis of dataset 1 ( and ).

When only the PEV of the estimated allele substitution effects are known, while PEC are not, then we can approximate with This approximation would be accurate if the matrix product has (close to) zero off-diagonal elements, which is dependent on the characteristics of genotypes in dataset 1 (*e.g.*, allele frequencies, linkage disequilibrium (LD), and population/family structure). If this is not the case, the approximation will bias the analysis by ignoring off-diagonal elements.

When allele frequencies and LD correlations in dataset 1 are known, we can obtain a good approximation of under some conditions (one phenotype record per individual, homogenous residual variance, overall mean is the only fixed effect, and Hardy-Weinberg equilibrium). Derivations in Appendix A2 show that, under these conditions, we can approximate with with the unknown matrix approximated from commonly available population parameters (*i.e.*, allele frequencies and LD correlation) as where p is a vector of allele frequencies, V is a diagonal matrix of expected genotype sum of squares with the i-th diagonal element equal to , and C is a matrix of pairwise genotype correlations between markers. In practice, the matrix C for dataset 1 could be unknown, but we can approximate it by using a reference panel that includes, for example, available genotypes of nonphenotyped individuals originating from this population (Yang *et al.* 2012; Vilhjálmsson *et al.* 2015; Maier *et al.* 2018).

Finally, we relax the assumption of having a single phenotype record per individual in the preceding approximations. This is relevant when individuals have repeated phenotype records, *e.g.*, repeated longitudinal production or reproduction records in livestock or replicated field trials in crops. A related issue is the violation of assumption of homogenous residual variance when phenotype records are first preprocessed and then used in genomic analyses, *e.g.*, deregressed progeny proofs in livestock (*e.g.*, Garrick *et al.* 2009) or adjusted field trial means in crops (*e.g.*, Schulz-Streeck *et al.* 2013; Oakey *et al.* 2016; Damesa *et al.* 2017). For these situations, we show in Appendix A3 that we can approximate with where Ψ is a diagonal matrix with the j-th diagonal element equal to , and is a diagonal matrix with the *j*-th diagonal element representing the square root of effective number of records for the *j*-th marker. The matrix can be obtained by solving the nonlinear system of equationsthrough a fixed-point iteration algorithm (Burden and Faires 2010) detailed in Appendix A3. It is worth noting that the proposed algorithm requires the inversion of a dense matrix at each iteration. This computational cost can be reduced by performing the algorithm for each chromosome separately.

#### Integration with multiple populations:

When more than two populations or datasets are available, the developed methods can be easily extended. With n datasets, the prior distribution for allele substitution effects in the separate analysis of the n-th dataset is defined using the posterior distributions for allele substitution effects from the separate analyses of datasets:

### Simulations

We tested developed methods with simulated data that either had low or high genetic diversity. The data were simulated in five replicates with the AlphaSim program, which uses the coalescent method for simulation of base population chromosomes and the gene drop method for simulation of chromosome inheritance within a pedigree (Hickey and Gorjanc 2012; Faux *et al.* 2016).

A diploid genome was simulated with 30 chromosomes, each 10^{8} bp long. Coalescent mutation and recombination rate per base pair were set to 10^{−8}, while effective population size was modeled over time to mimic population history of a livestock population in line with the values reported by MacLeod *et al.* (2013). Specifically, for the low diversity scenario, the effective population size of the base population was set to 100 and increased to 120, 250, 350, 1000, 1500, 2000, 2500, 3500, 7000, 10,000, 17,000, and 62,000 at, respectively, 6, 12, 18, 24, 154, 454, 654, 1754, 2354, 3354, 33,154, and 933,154 generations ago. For the high diversity scenario, effective population size of the base population was set to 10,000 and increased above this value in the same way as in the low diversity scenario; to 17,000 and 62,000 at 33,154, and 933,154 generations ago. For each chromosome, 10,000 whole chromosome haplotypes were sampled, which, on average, hosted ∼700,000 markers (21 million per genome) for the low diversity scenario and 1,400,000 markers (42 million per genome) for the high diversity scenario. Out of these loci, 100 per chromosome (3000 per genome) were sampled as causal loci affecting a complex trait. The allele substitution effect of causal loci was sampled from a normal distribution with mean zero and variance 1/3000. The effects were used to simulate a complex trait with additive genetic architecture. In addition, 2000 loci per chromosome (60,000 per genome) were selected as markers with the restriction of having minor allele frequency above 0.05.

From the base population, founder genomes for four populations (A, B, C, and D) were obtained by random sampling of chromosomes with recombination. The populations were ancestrally related through the common base population, but otherwise maintained independently, *i.e.*, there was no migration between the four populations. Each population was initiated with 10,000 founders (half males and half females) and maintained for seven generations with constant size. In the low diversity scenario, with the effective population size of 100, 25 males and 5000 females were selected as parents of each generation, while in the high diversity scenario, with the effective population size of 10,000, all 5000 males and 5000 females were used. The 25 males were selected on true genetic value, assuming accurate progeny test was available.

For every individual in the population we simulated two types of phenotypes. First, an own single phenotype was simulated as the sum of the true genetic value and a residual sampled from a normal distribution with mean zero and residual variance scaled relative to the variance of true genetic value in the base population such that heritability was 0.3. These simulated single phenotype records mimic records measured on the individual. Second, a weighted phenotype was simulated as the sum of the true genetic value and the mean of residuals. Each residual was sampled from a normal distribution with mean zero and residual variance scaled relative to the variance of true genetic value in the base population such that heritability was 0.3. The weight was equal to where the real value was sampled from a geometric distribution with a probability p of 0.15 and a probability mass function of with . The average was 6.6. These weighted phenotypes mimic either repeated records of an individual or records on multiple progeny of an individual. To satisfy the assumption of identical residual variance across all analyses, phenotype records were divided by the residual SD specific for each population, such that . For every individual in each population we stored the true genetic value, own single and weighted phenotype records, associated weight, and 60,000 marker genotypes.

### Analysis

The data were analyzed in several ways to evaluate the developed methods. In each case the aim was to obtain accurate genetic values utilizing all the available information. Specifically, we integrated results from separate analysis of populations B, C, and D, into the analysis of population A. We assumed throughout that variance components were known and equal to the rescaled variances. We analyzed three scenarios in total. The first and second scenario used population specific training data of randomly sampled 30,000 individuals with single phenotype record from generations 1–6 under low and high diversity settings. The third scenario used population specific training data of randomly sampled 10,000 individuals with weighted phenotype record from generations 1–6 under low diversity setting. In all scenarios all of the 10,000 individuals from generation 7 of each population were considered as validation individuals. The following analyses were performed:

A joint analysis of four populations. This was the reference that the other analyses were compared against;

A separate analysis for each of the four populations;

An exact integration of separate analyses of populations B, C, and D, into the analysis of population A;

The same as 3, but approximating the PEC matrix with a partial PEC matrix for each chromosome,

*i.e.*, PEC between markers on different chromosomes were set to zero;The same as 3, but approximating the PEC matrix with a diagonal PEV matrix,

*i.e.*, PEC between all markers were set to zero;The same as 3, but approximating the PEC matrix with PEV, allele frequencies, and LD correlations between markers obtained from the training sets. For the scenario with weighted phenotype records, the algorithm for estimating the effective number of records per marker was performed for each marker separately and for each chromosome separately.

The same as 6, but with LD correlations between markers computed from validation individuals instead of the training data.

For each analysis we calculated genomic prediction accuracy as the Pearson correlation between the true and estimated genetic value in validation individuals. Further, we evaluated the different integrations by comparing estimated genetic values of validation individuals against the estimated genetic values obtained from the joint analysis, which was considered as the reference because it used information from all populations. If integration was fully accurate, there should be no difference between the joint analysis and the analysis with integration. We assessed this by (a) accuracy of integration as a Pearson correlation between estimated genetic values from the joint analysis and the analysis with integration (desired value equals 1), (b) calibration of integration as a regression of estimated genetic values from the joint analysis on estimated genetic values from analysis with integration, and (c) magnitude of error in integration as a mean square error (MSE) between estimated genetic values from the joint analysis and from the analysis with integration (desired value equals 0). By calibration, we mean the slope of relationship of the estimates from the integration analysis onto the estimated genetic values from the joint analysis. The desired slope value is 1, which indicates a well calibrated model. Values above or below 1 indicate an uncalibrated model.

### Data availability

Supplemental figures are available in Supplemental Material, File S1. A description of the simulated genotype and phenotype datasets for each scenario is provided in File S2. Simulated genotype and phenotype datasets for the five replicates of each scenario are provided in Files S3–S5. Data simulation scripts and Fortran codes developed to perform the different analyses, as well as a short description of each of them, are provided in File S6. Supplemental material available at Figshare: https://doi.org/10.25386/genetics.6216533.

## Results

### Genomic prediction accuracy of separate and joint analyses

Joint analysis increased genomic prediction accuracy in comparison to separate analyses. This is shown in Table 1. Analyzing separately the four datasets gave accuracies of ∼0.71 (low diversity) and 0.53 (high diversity) with single phenotype records, and of ∼0.73 (low diversity) with weighted phenotype records. Analyzing jointly the four datasets increased accuracy by at least 0.09 absolute points with single phenotype records and by at least 0.12 absolute points with weighted phenotype records.

### Integration based on PEC, partial PEC, or PEV matrices

For all scenarios, the developed method enabled exact integration when complete PEC matrices were used. Integration of estimated allele substitution effects by means of the complete PEC matrix led to the same estimated genetic values as with the joint analysis, as shown by correlation and regression coefficients of 1, and MSE close to 0 (Figure 1, Figure 2, Figure 3, Figure 4, and Figures S1–S8). For comparison, correlations between estimated genetic values from separate analyses and joint estimated genetic values were ∼0.87 (low diversity) and 0.77 (high diversity) with single phenotype records, and 0.85 (low diversity) with weighted phenotype records.

Approximate integration by means of partial PEC matrices for each chromosome, that is ignoring PEC between markers on different chromosomes, gave almost as accurate and calibrated estimated genetic values as the exact integration. This is illustrated in Figure 1, Figure 2, Figure 3, Figure 4, and Figures S1–S8 with correlations higher than 0.96, regression coefficients close to 1, and MSE close to 0. Increasing the diversity slightly deteriorated accuracy and calibration of genomic predictions (Figure 1, Figure 2, and Figures S1–S4).

Approximate integrations by means of PEV matrices, that is ignoring PEC between all markers, gave quite accurate, but not calibrated estimated genetic values. This is shown in Figure 1, Figure 2, Figure 3, and Figure 4 and in Figures S1–S8. Correlations between joint estimated genetic values and estimated genetic values with integration by means of PEV were between 0.95 and 0.98 with single phenotype records and between 0.93 and 0.95 with weighted phenotype records. Despite these correlations close to 1, estimated genetic values were not well calibrated, as depicted by regression coefficients below 0.77 for the low diversity scenarios with single and weighted phenotype records, and below 0.86 for the high diversity scenario with single phenotype records (Figure 2, Figure 4, and Figures S2 and S6).

### Integration based on PEV, allele frequencies, and LD information

When LD information was derived from training data of other populations, approximate integrations by means of PEV, allele frequencies, and LD information, resulted in highly accurate and well calibrated estimated genetic values with single phenotype records. This is shown in Figure 1 and Figure 2 (Figures S1–S4). Correlation and regression coefficients were equal to 1 for the low diversity scenario. Slightly lower values, but still close to 1, were observed for the high diversity scenario. For both low and high diversity scenarios, MSE were close to 0. In contrast, when LD information was derived from validation data of other populations, approximate integrations gave less accurate and calibrated estimated genetic values. This is shown in Figure 1 and Figure 2 (Figures S1–S4). For these scenarios, correlations were equal to at least 0.94, and regression coefficients varied between 0.87 and 1.05.

For the scenario with weighted phenotype records, approximate integrations by means of LD information from training data of other populations resulted in highly accurate and well-calibrated estimated genetic values when sets of markers per chromosome were used to estimate the effective number of records for each marker. Correlations between joint estimated genetic values and estimated genetic values with integration were ∼0.99 (Figure 3 and Figure S5), regression coefficients were ∼0.95 (Figure 4 and Figure S6), and MSE were close to 0 (Figure S7 and Figure S8). Using LD information from the validation data of other populations, instead from the training data of other populations, gave slightly less accurate (correlations higher than 0.95), and moderately less calibrated estimated genetic values (regression coefficients between 0.87 and 1.04; Figure 3, Figure 4, and Figures S5–S8). For both cases, estimating the effective numbers of records per marker, instead of for all markers per chromosome simultaneously, reduced accuracy and calibration of estimated genetic values (Figure 3, Figure 4, and Figures S5 and S6).

### Comparison of estimated allele substitution effects

Correlation and regression coefficients between estimated allele substitution effects from the joint analysis and analysis with integration largely followed patterns of the corresponding values for estimated genetic values (Table 2 and Table 3). Correlation and regression coefficients were close to 1 when the integration of estimated allele substitution effects was by means of the complete PEC matrices. Ignoring PEC between markers on different chromosomes, or ignoring PEC between all markers, reduced correlations to between 0.92 and 0.99 (Table 2 and Table 3). Using LD information with PEV led to correlations between joint estimates of allele substitution effects and estimates with integration ranging from 0.71 to 0.83 for the scenario with weighted phenotype records (Table 2 and Table 3).

## Discussion

The results show that the developed method enables accurate and well-calibrated estimated genetic values for complex traits using both individual-level data and summary statistics. As expected from theory, the analysis of individual-level data and estimated allele substitution effects from other analyses by means of PEC matrices, yielded the same estimates as the joint analysis of all individual-level data. To our knowledge, this is the first time that individual-level data and summary statistics were analyzed simultaneously for genomic predictions. As illustrated by simulations, the combined analysis of multiple datasets may increase genomic prediction accuracy over separate analyses of a single dataset. Unfortunately, combining individual-level data from several sources is generally not feasible for several reasons, *e.g.*, political roadblocks, data protections concerns, or data inconsistencies (Powell and Sieber 1992; Vilhjálmsson *et al.* 2015; Maier *et al.* 2018). However, summary statistics, such as estimates of allele substitution effects and associated measures of accuracy (*e.g.*, PEV), are usually available for exchange in human genetics, or are discussed to be shared, *e.g.*, at an international level for dairy cattle breeding (Liu and Goddard 2018). The developed method enables increase in genomic prediction accuracy of complex traits by means of jointly analyzing the available individual-level data and summary statistics.

Accurate integration of estimated allele substitution effects is possible also when the complete PEC matrix is not available. This is important because computing the exact PEC matrix and exchanging it between analyses might be challenging in some cases. For the vast majority of marker arrays used in animal and plant breeding, the calculations and data transfers should be doable. For example, most arrays have between 10,000 and 100,000 markers, for which we need between ∼1 and ∼80 GB of memory to store the PEC matrix and between a minute and a day to invert it on current computers. For a larger number of markers, commonly used in human genetics, the memory requirements and computing time become prohibitive. The results show that in such cases we can still obtain accurate genomic predictions when the integration is done by means of partial PEC matrices for each chromosome. This is expected since high LD between markers mostly occurs within chromosomes. High LD between markers on different chromosomes may especially occur in structured populations and populations under selection (Farnir *et al.* 2000; Flint-Garcia *et al.* 2003; Rostoks *et al.* 2006). Both of these conditions are present in breeding populations. However, the results suggest that LD between chromosomes can be ignored for the purpose of integration for populations with both low and high diversity. The results also show that we can successfully integrate estimated allele substitution effects when only PEV and allele frequencies from each population are available together with LD information of a reference genotype panel representative of each population. Assuming that such reference genotype panels are available, only estimated allele substitution effects, associated PEV, and allele frequencies need to be exchanged between populations for such analyses. Similar conclusions were drawn from studies combining only summary statistics obtained from genome-wide association studies to perform multi-trait genomic predictions (Maier *et al.* 2018).

Accurate integration of estimated allele substitution effects is possible irrespective of the diversity of the populations and characteristics of genotypes (*e.g.*, allele frequencies, LD). This is obvious, and confirmed by our results, when integration is performed by means of complete PEC matrices. When complete PEC matrices are unavailable, accurate integration is possible if the inverses of the PEC matrices can be approximated accurately from available population parameters (*i.e.*, LD and allele frequency information), whatever the level of diversity and characteristics of the populations, as shown by our results or a study combining summary statistics in human genetics (Maier *et al.* 2018). In our study, the population parameters obtained from the reference panels adequately reflected the characteristics of the training sets. We expect that this would be the case for populations with substantial migration, such as, for example, Holstein dairy cattle populations. Future studies should be conducted to assess the impact of suboptimal reference panels. Therefore, the developed method is expected to perform well on any type of data, from animal and plant breeding to human genetics, provided accurate information is available.

The developed method has some simplifying assumptions that can be readily relaxed. For example, we assumed that the same genotype coding was used in all populations. This assumption can be relaxed when centered genotype coding (*i.e.*, of the form of ) is used because variance component estimates, estimates of allele substitution effects and PEC are the same irrespective of the centering of the genotype coding, provided that the model has a fixed general mean, which is considered in the integration (Strandén and Christensen 2011). Also, centered and scaled (standardized) genotype coding is often used in human genetics, instead of only centered genotype coding (Yang *et al.* 2010; Speed *et al.* 2012; Maier *et al.* 2018). In practice, estimates of genetic values are only slightly influenced by scaling of centered genotype coding (Strandén and Christensen 2011; Bouwman *et al.* 2017). Therefore, assuming that the same estimated genetic values are obtained with different scaling, allele substitution effects estimated using one type of genotype scaling could be obtained from a postanalysis by converting estimated genetic values computed for a reference genotype panel into allele substitution effects for another genotype scaling. Converting estimated genetic values into allele substitution effects is often referred to as back-solving of allele substitution effects (Strandén and Garrick 2009; Strandén and Christensen 2011; Wang *et al.* 2012; Bouwman *et al.* 2017). PECs associated with the converted estimated allele substitution effects could be derived from the (prediction error) covariances of the estimated genetic values (see derivations in Appendix A4).

Allele substitution effects estimated from analyses using different sets of markers or different residual variances, can be used in the integration as well. The assumption that all individuals were genotyped at the same loci could be considered as fulfilled if small differences in the sets of markers are corrected by assuming zero allele substitution effect and zero accuracy for markers not used in an analysis. When large differences between sets of markers are observed, this assumption can be accommodated following two approaches. A first, postanalysis, approach consists of assuming that estimated genetic values are the same for two different sets of markers, allowing the conversion of estimated allele substitution effects from one set of markers to another set of markers (Liu and Goddard 2018). The conversion can be performed by back-solving estimated allele substitution effects from estimated genetic values, as proposed previously for different genotype codings, or by applying a marker model to the estimated genetic values with the reference set of markers (Liu and Goddard 2018). A second approach consists of harmonizing genotype data across populations. This approach must be performed before the analyses, and requires therefore coordination between populations. Harmonization of genotype data could be performed by identifying a subset of markers for which all populations are genotyped, or by genotype imputation (*e.g.*, Marchini and Howie 2010). Finally, the assumption that residual variances were the same in all populations, can be relaxed by noting that separate estimates of allele substitution effects obtained by the system of equations (2), can be also obtained by the following different formulations:where is the residual variance used for the i-th (focal) analysis, and

For integration of must be approximated using the residual variance of the focal population () and the effective numbers of records per marker estimated using variance components of the *i*-th analysis. Another way to relax this assumption is to extend our univariate model to a bivariate model, similarly to methods developed to combine different genetic evaluations in animal breeding (Schaeffer 1994; Vandenplas *et al.* 2015). In a bivariate model, one trait would represent individual-level data, while the other trait would represent summary statistics. The genetic correlation between the two traits could be estimated based on a subset of individual-level data available for both datasets or based on summary statistics (Bulik-Sullivan *et al.* 2015). Such an approach would also allow the integration of summary statistics expressed on a different scale (*e.g.*, different measure units, trait definitions) than the scale of the focal population (Vandenplas *et al.* 2015).

The developed method can be readily generalized to multi-trait models and is therefore a generalization of previous works that were based on several (implicit) assumptions (Liu and Goddard 2018; Maier *et al.* 2018). For example, previous works assumed that no individual-level data were available. It was also (implicitly) assumed that only single phenotype records with homogeneous residual variance (Maier *et al.* 2018), or that the least-squares part of the separate analyses (Liu and Goddard 2018), were available for integrating estimated allele substitution effects. Both assumptions lead to simple and accurate approximations of PEC matrices as shown in our study. However, we relax all these assumptions, such that our method can jointly analyze individual-level data and summary statistics, with possibly multiple phenotype records per individual.

With all the proposed generalizations, the developed method could be used in different contexts. For example, in human genetics, allele substitution effects with associated SE are publicly available (Yang *et al.* 2012; Vilhjálmsson *et al.* 2015; Maier *et al.* 2018). In animal breeding, individuals′ genetic values with associated reliabilities are publicly available and in the case of dairy cattle extensively combined across multiple populations (Schaeffer 1994; VanRaden and Sullivan 2010; Jorjani *et al.* 2012; Vandenplas *et al.* 2017). The developed method can be used in both contexts, but, in the latter case, individuals′ genetic values must be first back-solved to allele substitution effects (Strandén and Garrick 2009; Strandén and Christensen 2011; Wang *et al.* 2012; Bouwman *et al.* 2017). It is worth noting that our method assumes that summary statistics from one population are free of information from other populations. This suggest that it can be used when there is no, or limited, sharing of information between populations, as is, for example, the case in beef cattle, but not in dairy cattle populations such as Holstein, where pseudophenotypes summarizing information from multiple populations are used extensively (VanRaden and Sullivan 2010; Jorjani *et al.* 2012). This assumption can be relaxed by performing separate analyzes free of information from other populations, or by correcting for double-counting of information, which has bee developed for the integration of estimated genetic values from different populations (Vandenplas *et al.* 2014, 2017; VanRaden *et al.* 2014). This correction for double-counting of information is not yet developed for the integration of summary statistics, and should be investigated in future studies.

## Conclusions

We developed a method for genomic prediction that accurately integrates summary statistics obtained from analyses of separate populations into an analysis of individual-level data. The method accommodates use of multiple phenotype (pseudo)records per individual, and further extensions have been presented to accommodate for differences in residual variances or genotype codings used in the populations. When complete summary statistics information is available the method gives identical genomic predictions as the joint analysis of individual-level data from all populations. When summary statistics information is not complete we can use a series of approximations that give very accurate and well-calibrated genomic predictions.

## Acknowledgments

This study was financially supported by the Dutch Ministry of Economic Affairs (TKI Agri & Food project 16022), the Breed4Food partners Cobb Europe, CRV, Hendrix Genetics and Topigs Norsvin, and United Kingdom Biotechnology and Biological Sciences Research Council (BBSRC) Institute Strategic Programme Grant (ISPG) to The Roslin Institute BBS/E/D/30002275. The use of the high performance computing (HPC) cluster has been made possible by CAT-AgroFood (Shared Research Facilities Wageningen University and Research).

Author contributions: J.V. derived the equations, wrote the programs to do the analyses, performed the analyses, and drafted the outline of the manuscript. G.G. performed the simulations. All authors discussed the design of the simulations. J.V. and G.G. wrote the first version of the manuscript. All authors provided valuable insights throughout the analysis and writing process.

## Appendix A1: Exact Integration

Here, we detail the derivation of exact integration by means of absorbing the set of equations that pertain to one dataset. We start with the system of equations for separate analysis of dataset 1:

and the system of equations for the joint analysis of datasets 1 and 2:(A1.2)From the first set of equations in (A1.2) it follows:(A1.3)From the third set of equations in (A1.2) it follows: (A1.4)Inserting (A1.3) into (A1.4) gives, after some algebra:with .

Now the system of equations (A1.2) can be rewritten with the first set of equations absorbed as: (A1.4)Similarly, the absorption of the first set of equations in separate analysis of dataset 1 (A1.1) leads to:(A1.5)where(A1.6)is the inverse matrix of prediction error covariances of

Combining (A1.4) and (A1.5) with the use of (A1.6) enables the exact integration of estimates from the separate analysis of dataset 1 into the separate analysis of dataset 2 with the following system of equations:

(A1.7)## Appendix A2: Approximate Integration

Here, we detail the derivation of different approximate integrations by means of simplified assumptions and use of summary statistics. We start with the expression for prediction error covariance matrix of allele substitution effects from dataset 1:

If we assume that: (1) every individual has a single phenotype record, *i.e.*, (2) residual variance is homogeneous, *i.e.* and (3) only overall mean is fitted as a fixed effect, *i.e.*, then we can simplify (A2.1) as:(A2.2)because will tend to the identity matrix I with increasing The matrix also known as the centering matrix, is a symmetric and idempotent matrix with off-diagonal elements equal to and with diagonal elements equal to

When genotypes from the dataset 1 are not available, but variance components and are, we “only” need to approximate the unknown matrix of genotype sum of squares in (A2.2). This product can be approximated from linkage-disequilibrium and allele frequency information of the dataset 1, as shown in the following (similarly to Yang *et al.* 2012, Vilhjálmsson *et al.* 2015, and Maier *et al.* 2018). Assume that LD between two markers is represented by the correlation of their unphased genotypes (Rogers and Huff 2009). Then, a matrix of all pairwise correlations between markers is:(A2.3)where the matrix contains centered genotypes of dataset 1 (). The matrix product can be computed as:(A2.4)where are allele frequencies in dataset 1 (Strandén and Christensen 2011). Assuming Hardy-Weinberg equilibrium, the i-th diagonal element of the matrix product is equivalent to expected genotype sum of squares at the i-th marker, with being the allele frequency of the i-th marker in dataset 1.

Combining (A2.3) and (A2.4) we can approximate the unknown matrix of genotype sum of squares as:(A2.5)where V is diagonal matrix of expected genotype sum of squares with the i-th diagonal element equal to

## Appendix A3: Estimation of the Effective Number of Records Per Marker

Here, we detail the algorithm for computing the effective number of records per marker by use of available population parameters (*i.e.* LD, and allele frequency information) and PEVs of of the dataset 1. We start with the expression for the PEC matrix of allele substitution effects from dataset 1:If the number of individuals and the number of records per individual are unknown, we can assume that a diagonal matrix exists such that:where Ψ is a diagonal matrix with the *j*-th diagonal element equal to , and the squared *j*-th diagonal element of represents the effective number of records for the *j*-th marker. The term is similar to the approximation of the unknown matrix of genotype sum of squares (*i.e.*, ) in Appendix A.2. However, it does not involve the number of individuals because it is confounded with the effective number of records.

The diagonal matrix can be estimated by solving the nonlinear system of equationsthrough a fixed-point iteration algorithm (Burden and Faires 2010) as follows:

1)

where is a diagonal matrix with the *i*-th diagonal element equal to the PEV of the *i*-th marker and contains the diagonal elements of

2)

3)

4)

5)

6)

7) If trace of is not sufficiently small:

a.

b. If any diagonal element in is negative, set it to 0

c.

d.

e. Repeat from 4)

8)

It is worth noting that the proposed algorithm is similar to algorithms to estimate effective number of records per individual, where “effective” means that they are free of contributions from relatives (Misztal and Wiggans 1988; Vandenplas and Gengler 2012). The *j*-th diagonal element of can therefore equivalently be considered as the effective number of records for the *j*-th marker.

## Appendix A4: Conversion of Allele Substitution Effects

Here we detail a postanalysis to obtain allele substitution effects estimated using one type of genotype coding by converting estimated genetic values computed for a reference genotype panel with allele substitution effects for another genotype coding . We assume that allele substitution effects are available with the associated prediction error (co)variance matrix , as well as the (co)variance matrix of , and genotypes of a reference panel using a particular type of genotype coding (). Estimates of genetic values for the reference individuals are obtained as

Assuming that estimated genetic values are not influenced by scaling of centered genotype coding (Strandén and Christensen 2011; Bouwman *et al.* 2017), and that the (co)variances of genetic values are the same irrespective of the genotype coding, we can write that with being a matrix with reference genotypes using another type of genotype coding than and being a vector of estimated genetic values using this type of genotype coding. Therefore, can be computed by back-solving as follows (Strandén and Garrick 2009; Wang *et al.* 2012; Bouwman *et al.* 2017):where is a diagonal matrix (*e.g.*, an identity matrix I) with optional different weights to differentially shrink different loci.

Based on the properties of mixed models (Henderson 1984), the prediction error covariance matrix of can be obtained as follows:

## Footnotes

Supplemental material available at Figshare: https://doi.org/10.25386/genetics.6216533.

*Communicating editor: W. Valdar*

- Received May 3, 2018.
- Accepted July 16, 2018.

- Copyright © 2018 by the Genetics Society of America