Genome position specific priors for genomic prediction

Data

Genotypic data was a mixture of true and imputed SNP markers from the Illumina 777K SNP chip. For HOL-AUS there were 843 Holstein heifers genotyped on the 777K SNP chip as well as 93 key ancestor bulls. For JER-AUS 93 key ancestor bulls were genotyped for the 777K SNP chip. Quality control steps included removing SNPs with very low minor allele frequencies, ambiguous or undefined map positions, and no heterozygote genotypes. For full details see [11]. These animals were used as reference to impute the high density genotypes for the remaining 2204 Holstein and 447 Jersey bulls which were genotyped with the 50k chip.

For HOL-NOR 557 bulls from the EuroGenomics project [4] were genotyped using the 777K chip and these bulls were used as reference to impute the 777K markers for the bulls genotyped with the Illumina 50K SNP chip. After imputation LD of each marker with the previous one in the assembly was inspected. If two adjacent markers were in complete LD one of the markers was deleted, so that r² of any pair of adjacent markers was less than one. The marker data was further edited by deleting markers with a minor allele frequency less than 0.01.

Imputation was done using Beagle [13] in all three populations. Since the purpose was to compare segments across populations and use this information for genomic prediction, the SNP datasets were further edited to only keep the markers that were in common across the populations. After data editing 465,542 markers remained for analysis.

Each of the datasets was split into a reference and validation set (Table 1) to allow for cross validation of the accuracy of DGV. In HOL-NOR the bulls were separated by birth date before or after 2001-10-01, and in JER-AUS and HOL-AUS the bulls were split by onset of progeny test before or after 2007. In both cases the younger animals were assigned to the validation set. This cross validation strategy was chosen as the resulting accuracy is the most meaningful in the context in which the genomic predictions will be used: in the dairy industry. Here reference sets of older bulls are used to predict the DGV for young bulls which are then selected for use based on these DGV. In the Australian data all of the available data was used to estimate segment variances to maximize the data. In the Nordic dataset only the reference set was used to estimate segment variances.

All genotypic and phenotypic data was obtained from pre-existing routine genetic evaluation data for the dairy cattle populations and required no ethical approval.

Estimation of genetic variances explained by different segments

Where A is the additive relationship matrix, σ_a² is the variance of residual polygenic effects, and r _y ² is the reliability of DRP/DTD. The four-distribution mixture chosen for the SNP effects, does not reflect any biological hypothesis, but was chosen to allow for easier mixing between SNPs with no effect and SNPs with effects of different sizes. The Dirichlet prior on the proportions of different SNP variances with all parameters set to one, is actually a uniform prior, but specifying it in this manner reflects the fact that the posterior distribution on the proportions follows a Dirichlet distribution with a pseudo count of 1 from each of the four distributions. The prior is not uninformative in any statistical sense since it states that all distributions have equal probabilities, but it adds very little information compared to the posterior, as the data gives information on almost half a million counts, and the prior only adds 4, see [11] for more detail.

Where W _s is the sub-matrix of W corresponding to the SNPs in segment s, and g _s is the vector of estimated SNP marker effects for the same segment, i.e. the segment variance is the variance across individuals of the partial direct genetic values (DGVs, marker only estimated breeding values) belonging to segment s. Variances of the partial DGVs for all segments were calculated at each iteration in the Gibbs sampler, and the estimated segment variances were obtained as the posterior means. Segment variances were estimated for segment sizes of 10, 25, 50, 100, 250, 500, 1000, 2000 or 3000 SNPs and for entire chromosomes. The approach is similar to [14] where sliding windows of five consecutive SNPs are used to estimate the genetic variance of chromosomal regions. In our approach the windows are however not overlapping.

Posterior means of the parameters were obtained using a Gibbs sampler run for 20,000 iterations with a burn-in of 10,000 in the Holstein populations. For the Jersey population results were not consistent with only 20,000 iterations, so a chain length of 100,000 with a burn-in of 50,000 was used instead. The relatively poor mixing properties of the Gibbs sampler for the Jersey data could be due to the small size of the reference population. Lengths of the chains were based on preliminary runs and comparisons of the obtained segment variances. With 20.000 iterations the Holsteins showed a mean pairwise correlation between segment variances from 10 consecutive runs of 0.95, whereas the Jerseys showed a mean correlation between segment variances from 10 consecutive runs of 0.80. Increasing the number of iterations for the Jerseys to 100.0000 increased the mean correlation of segment variances between consecutive runs to .96.

Prediction using location specific prior information

Here π _s is the probability vector for the four SNP effect distributions in segment s, and α_s is the vector of prior parameters for the Dirichlet distribution in segment s. The model is similar to the original BayesR model, with the modification that the probability to sample SNPs from the four different distributions now depends on the segment. By setting the location specific information via the Dirichlet prior, instead of using constant proportions, the model estimates the proportions using both the data and the prior information. As this is a BayesR by segment approach, the model will be referred to as BayesRS.

HOL-AUS (random) is a random subset of 500 animals from the HOL-AUS reference population, which was generated to test the hypothesis that the advantage of the BayesRS model would be greater in smaller populations. The second and third setups were tested using the same validation animals.

Validation of DGV accuracy

Where w '  _k is the row of W belonging to animal k. Accuracies of the DGV were calculated as r(DGV, DTD) and validated in HOL-AUS and JER-AUS. Differences in accuracy between BayesR and BayesRS were tested for significance using a Hotelling-Williams t-test, which takes account of the number of individuals in the validation set [15].

Segment variances

Correlations of segment variances between populations are large if the DGAT1 segment is included since the proportion of variance explained for this QTL is very large compared to all others, as seen from Table 3 and Figure 1. Hence SNPs associated with the DGAT1 mutation have been removed from all plots in Figure 2 for clarity. For fat yield the QTL on chromosome 5 with large effect has been removed as well. Figure 2 shows a plot of segment size versus correlations of segment variances between the three populations for the three traits. With the large effect segments removed, the patterns of correlations generally follow expectations based on knowledge of genetic relationship between populations. Highest correlations are found between HOL-AUS and HOL-NOR, and second and third place are respectively HOL-AUS versus JER-AUS and HOL-NOR versus JER-AUS. In Australia some crossbreeding has taken place between the Holstein and Jersey populations and, accordingly, correlations of variance explained by the segments are higher for these two populations than between HOL-NOR and JER-AUS. Genotype by environment interactions could also influence these correlations.

Correlations for small segment sizes are close to zero. A possible explanation for this is that differences in LD patterns and SNP allele frequencies across breeds cause the SNPs with the highest associations to actual QTL to be placed in different segments when these are very small. The rapid increase in correlation of segment variances with segment size for HOL-AUS vs. HOL-NOR suggests that these population share QTL in similar locations, as would be expected given they are genetically closely related. For this pair of populations, even reasonably small segments would convey information between the populations. When taking an entire chromosome as a segment, correlations of segment variances ranging from 0.6 between JER-AUS and HOL-NOR up to 0.8 between the HOL-AUS and HOL-NOR were found. A connection between chromosome size and variance explained has previously been reported by e.g. [20]. As seen from Figure 3 a similar pattern can be found in our populations, provided the effect of DGAT1 is removed from the analysis.

BayesRS

The accuracies of DGV, measured as r(DGV,DTD) in the validation population, are shown in Figure 4 and 5. In the figures, a horizontal line depicts the accuracy using BayesR. For JER-AUS a second (higher) horizontal line shows the accuracies obtained from BayesR when simply pooling JER-AUS with HOL-AUS. For the HOL-AUS with HOL-NOR information scenario this accuracy was not available due to data sharing policies.

For JER-AUS no gain in accuracy was observed for milk yield when using prior information from HOL-AUS, for protein yield a small gain of around 1% is seen for the smallest segment size, and for fat yield gains in accuracy of up to 3.5% are seen when using the genomic location specific prior information compared to using BayesR. Compared with accuracies obtained with a simple pooling of reference data, the BayesRS approach leads to an extra gain of up to 1.5% for fat yield, but not for the other two traits. Although differences in accuracy were seen, none of the differences were significant at a 5% level, reflecting the small size of validation population.

For HOL-AUS the largest gain in accuracy is found for protein yield with gains of up to 2%. For milk and fat yield smaller gains are seen, and these are not significantly different for milk yield. Using the prior information derived from HOL-NOR, however, seems consistently better than the model without location specific priors. Results from the HOL-AUS (random) setup are shown in Figure 6. The gain in accuracy is here slightly higher for protein and fat yield, whereas the gain for milk yield seems unchanged. Results were significant for protein and fat but not for milk yield. Looking at the accuracies of DGV from the regular BayesR model (no prior information from other populations) the mean accuracy is now on the same level as in JER-AUS due to the same size of the reference population (though the traits vary considerably).

In all three scenarios the highest gains in accuracy are found for a segment size of 100 markers, implying that using smaller segments gives a stronger advantage from the location specific priors. Furthermore, significant results are only found in two cases: the largest and smallest segments. For the largest segment size of 3000 markers, it is surprising that the increase in accuracy is significant although larger gains in accuracy are seen for smaller segments. However, this could be an artifact of the test chosen for the significance. With a large segment size the added information becomes very unspecific which could lead to results more similar to those obtained from the regular BayesR method. With a high correlation between DGVs from the two methods, the Hotelling Williams t-test would cause even small differences in accuracy to be significant.

The different scaling factors (weights) applied to the parameters in the Dirichlet priors, seems to make little or no difference on the accuracy of the BayesRS model, which suggests that the accuracies obtained from BayesRS could be random fluctuations. This is in many, but not all cases, supported by the lack of significance of the results.

To summarize, BayesRS gave accuracies comparable to, but not always higher than or significantly different from, a simple pooling of the data. For closely related populations pooling is expected to be superior. So a simple pooled multi-breed or multi-population reference could be a better approach in some cases, but not necessarily for all traits. For example, the BayesRS approach gave higher accuracies than a pooled reference for fat yield in JER-AUS. Further studies are needed to confirm the validity of the results in a larger validation population.

One advantage of the method presented here is a large reduction in computational demand. Since the BayesRS model only uses very condensed information from the other population, the increase in memory demand is negligible, and the extra complexity of the model only slightly increases the CPU run time. For JER-AUS running the BayesR model for 100.000 iterations required 33 hours, whereas the BayesRS model could be run for the same number of iterations in 39 hours. When using BayesR with the combined JER-AUS HOL-AUS reference data, 100,000 iterations takes about 150 hours, and more than quadruples the memory requirements.

Although the accuracies obtained using BayesRS in most cases cannot compete with pooling of the data, the results seem consistently better than when using only data from the target population and a non-informative prior, for example only the JER-AUS data. In some cases where the extra data itself is not available, the BayesRS model or a similar approach could improve the accuracy of genomic predictions using only summary statistics. This might be in cases when intellectual property issues prevents sharing of the raw data, but allows use of summary statistics as in this study. The approach could also be useful for meta-analysis of many data sets from different sources.

The model presented here would also allow the use of other prior information such as known QTL or expression pathways, by assigning a higher prior probability to sample large effects in the involved genomic regions. In this study segments were chosen arbitrarily with a fixed length, but another approach could be to define coding and non-coding regions of the genome as different segments and set different Dirichlet priors. A challenge here would, however, be how to choose the counts in the Dirichlet prior without sampling them from a different population. Previous results show that SNPs near genes found in both human and bovine genomes are significantly associated with stature [12]. By considering evolutionary conserved regions as segments the method using external information sources presented in this study could be used for genomic prediction across species for traits of common interest such as growth in meat-production animals or production traits in dairy species.